Common Mistake

Browsing through code that uses ojAlgo it’s quite common to see it being used inefficiently or at least not the intended way. Here’s an example code line, seen in the wild, that demonstrates common problems. Primitive64Matrix selectedColumns = Primitive64Matrix.**_FACTORY_**.rows(initialMatrix.logical().limits(-1, 5).get().toRawCopy2D());
That single line, combined with how the selectedColumns matrix was later used, demonstrates several common mistakes. Let’s try to explain using a code example.

Example Code

CommonMistake.java
import org.ojalgo.OjAlgoUtils;
import org.ojalgo.matrix.MatrixR064;
import org.ojalgo.matrix.MatrixR064.DenseReceiver;
import org.ojalgo.matrix.decomposition.SingularValue;
import org.ojalgo.matrix.store.RawStore;
import org.ojalgo.netio.BasicLogger;
import org.ojalgo.random.Uniform;

/**
 * An example with common mistake(s) seen in code using ojAlgo.
 *
 * @see https://www.ojalgo.org/2021/08/common-mistake/
 */
public class CommonMistake {

    private static final Uniform RANDOM = Uniform.standard();

    public static void main(final String[] args) {

        BasicLogger.debug();
        BasicLogger.debug(CommonMistake.class);
        BasicLogger.debug(OjAlgoUtils.getTitle());
        BasicLogger.debug(OjAlgoUtils.getDate());
        BasicLogger.debug();

        /*
         * Just some initial matrix.
         */
        MatrixR064 initialMatrix = MatrixR064.FACTORY.makeFilled(7, 9, RANDOM);

        /*
         * This line is copied from a GitHub repository using ojAlgo. It has several problems and, I believe,
         * demonstrates some things that are not obvious to new users.
         */
        MatrixR064 selectedColumns = MatrixR064.FACTORY.rows(initialMatrix.logical().limits(-1, 5).get().toRawCopy2D());
        /*
         * Further on in the code an SVD was calculated using the matrix of those selected columns.
         */
        SingularValue<Double> svd = SingularValue.R064.make(selectedColumns);
        svd.decompose(selectedColumns);

        /*
         * Let's do this in steps to see what actually happens
         */

        /*
         * This used to create a thin wrapper exposing additional "logical" API - without copying. Now this
         * method is deprecated and entirely redundant as it simply returns "this".
         */
        MatrixR064 logical = initialMatrix.logical();
        /*
         * A thin wrapper specifying that the number of columns should be limited to 5. (The rows limit is set
         * to -1 resulting in "no limitation". That's an alternative to setting the rows limit to the
         * full/actual number of rows.)
         */
        MatrixR064 limited = logical.limits(-1, 5);
        /*
         * Previously used to get the "matrix" back from the "logical builder" – now redundant and simply
         * returns "this".
         */
        MatrixR064 gotten = limited.get();
        /*
         * This extracts the elements as a 2D double array. That's another copy.
         */
        double[][] raw2D = gotten.toRawCopy2D();
        /*
         * This builds the final matrix. A third copy.
         */
        selectedColumns = MatrixR064.FACTORY.rows(raw2D);

        /*
         * I hope it's clear that extracting the double[][] was completely unnecessary resulting in 2
         * additional copies. Try to stop thinking in terms of double[][] ! It will prevent you from doing
         * that kind of thing. Instead simply...
         */
        selectedColumns = gotten;
        /*
         * And if you really do want another copy there is a more direct way to do it.
         */
        selectedColumns = MatrixR064.FACTORY.copy(gotten);

        /*
         * The SVD instance has its own internal storage. When you perform a decomposition another copy is
         * created.
         */
        svd.decompose(gotten);

        /*
         * The ojAlgo design/feature set allows to get rid of the explicit copy used to calculate the SVD. A
         * selected columns copy needs to be created somehow, but that copy can be made directly to the
         * internal storage of the SVD instance. Simply feed the thin wrapper representing the columns
         * selection of the original matrix to the SVD instance.
         */
        svd.decompose(limited);

        /*
         * This, original, code...
         */
        MatrixR064 unnecessary = MatrixR064.FACTORY.rows(initialMatrix.logical().limits(-1, 5).get().toRawCopy2D());
        svd.decompose(unnecessary);
        /*
         * ...is equivalent to this. That is less to type, less to execute, and less (essentially no) data
         * copying.
         */
        svd.decompose(initialMatrix.limits(-1, 5));
        /*
         * Prior to v50 you had to do it this way:
         */
        svd.decompose(initialMatrix.logical().limits(-1, 5));

        /*
         * Thinking in terms of double[][] will cause developers to create matrices this way.
         */
        double[][] newData = new double[7][9];
        for (int i = 0; i < newData.length; i++) {
            for (int j = 0; j < newData[i].length; j++) {
                newData[i][j] = RANDOM.doubleValue();
            }
        }
        MatrixR064 newMatrix = MatrixR064.FACTORY.rows(newData);

        /*
         * The correct, ojAlgo, way to do it would be something like this. (Apart from being less code, this
         * avoids that intermediate copy.)
         */
        newMatrix = MatrixR064.FACTORY.makeFilled(7, 9, RANDOM);
        /*
         * Since this example is using the immutable MatrixR064 class, manually/explicitly setting the
         * elements needs to be done using an intermediate mutable receiver. (Don't worry, there's no
         * unnecessary copying.)
         */
        DenseReceiver receiver = MatrixR064.FACTORY.newDenseBuilder(7, 9);
        receiver.fillAll(RANDOM);
        newMatrix = receiver.get();

        svd.decompose(newMatrix);

        /*
         * If you don't need the new matrix afterwards you can do it even more directly. Actually "getting"
         * the matrix can be handled internally by the SVD instance.
         */
        svd.decompose(receiver);

        /*
         * Integrating with other tools or existing code it is common to have data already in a double[][]. To
         * create an ojAlgo data structure this still does not need to be copied. There exists an ojAlgo class
         * that uses double[][] internally and allows wrapping existing data.
         */
        RawStore wrapped = RawStore.wrap(raw2D);
        svd.decompose(wrapped);
    }

}

Summary

  1. Stop thinking in terms of, converting to/from, double[][] and just use the ojAlgo data structures directly
  2. ojAlgo has special features that enable writing exceptionally memory efficient code. These feature are not well known/understood. To get the most out of ojAlgo you should become familiar with them.

Example – Memory Efficiency

ElementConsumersAndSuppliers.java

import static org.ojalgo.function.constant.PrimitiveMath.DIVIDE;
import static org.ojalgo.function.constant.PrimitiveMath.SUBTRACT;

import org.ojalgo.OjAlgoUtils;
import org.ojalgo.matrix.decomposition.QR;
import org.ojalgo.matrix.store.ElementsSupplier;
import org.ojalgo.matrix.store.R064Store;
import org.ojalgo.matrix.store.TransformableRegion;
import org.ojalgo.netio.BasicLogger;
import org.ojalgo.random.Normal;
import org.ojalgo.random.Uniform;

/**
 * A small example to show what can be done with {@link ElementsSupplier}s and {@link TransformableRegion}s.
 * To realy "see" the benefit of using these you should take this example, scale up the matrix dimensions and
 * extend the loop. Then run the code in a profiler and monitor memory allocation and garbage collection.
 *
 * @see https://www.ojalgo.org/2021/08/common-mistake/
 * @see https://github.com/optimatika/ojAlgo/wiki/Element-Consumers-and-Suppliers
 */
public class ElementConsumersAndSuppliers {

    static int ITERATIONS = 3;
    static int MATRIX_SIZE_SCALE = 1;

    public static void main(final String[] args) {

        BasicLogger.debug();
        BasicLogger.debug(ElementConsumersAndSuppliers.class);
        BasicLogger.debug(OjAlgoUtils.getTitle());
        BasicLogger.debug(OjAlgoUtils.getDate());
        BasicLogger.debug();

        Normal mormalDistribution = Normal.standard();
        Uniform uniformDistribution = Uniform.standard();

        // Assume you have the matrices [A],[B],[C] and [D]
        R064Store mtrxA = R064Store.FACTORY.make(5 * MATRIX_SIZE_SCALE, 7 * MATRIX_SIZE_SCALE);
        R064Store mtrxB = R064Store.FACTORY.make(7 * MATRIX_SIZE_SCALE, 9 * MATRIX_SIZE_SCALE);
        R064Store mtrxC = R064Store.FACTORY.make(5 * MATRIX_SIZE_SCALE, 9 * MATRIX_SIZE_SCALE);
        R064Store mtrxD = R064Store.FACTORY.make(9 * MATRIX_SIZE_SCALE, 5 * MATRIX_SIZE_SCALE);

        QR<Double> decompQR = QR.R064.make();

        // Imagine the matrices involved to be gigantic and the iterations (loop) many
        for (int l = 1; l <= ITERATIONS; l++) {
            // Inside this loop no array/vector/matrix memory is allocated,
            // only thin wrappers and placeholders

            BasicLogger.debug();
            BasicLogger.debug("Iteration {}", l);

            mtrxA.fillAll(uniformDistribution);
            mtrxB.fillAll(mormalDistribution);

            // [E] = [A][B]
            ElementsSupplier<Double> placeholderE = mtrxB.premultiply(mtrxA);
            // Note that the multiplication is not yet performed - it's a placeholder for a deferred operation.
            BasicLogger.debug("The matrix [E] = [A][B] will have {} rows and {} columns.", placeholderE.countRows(), placeholderE.countColumns());

            mtrxC.fillAll(uniformDistribution);

            // [F] = [A][B] - [C]
            ElementsSupplier<Double> placeholderF = placeholderE.onMatching(SUBTRACT, mtrxC);
            // Still, nothing is actually done - no calculations and no memory/array allocation.

            // [G] = [F]t
            ElementsSupplier<Double> placeholderG = placeholderF.transpose();
            // Still nothing...

            // [H] = [G] ./ 2.0 (divide each element by 2.0)
            ElementsSupplier<Double> placeholderH = placeholderG.onAll(DIVIDE.by(2.0));
            // Still nothing...
            BasicLogger.debug("The matrix [H] will have {} rows and {} columns.", placeholderH.countRows(), placeholderH.countColumns());

            // Now it happens all at once!
            placeholderH.supplyTo(mtrxD);
            // ...and you get to decide where to put the results.

            if (placeholderH.count() <= 100) {
                BasicLogger.debugMatrix("Resulting D", mtrxD);
            }

            // You can supply it directly to a decomposition's internal storage
            // (as you decompose it)
            decompQR.decompose(placeholderH);
            // The matrix that was decomposed never "existed" other than as an
            // intermediate state of the matrix decomposition
            // The elements of the intermediate placeholder H was supplied/written
            // directly to the decomposition's internal storage, and then further
            // worked on to perform the actual decomposition.

            BasicLogger.debug();
            BasicLogger.debug("Is it full rank? {}", decompQR.isFullRank());
        }

    }

}

Console Output

class ElementConsumersAndSuppliers
ojAlgo
2021-08-24


Iteration 1
The matrix [E] = [A][B] will have 5 rows and 9 columns.
The matrix [H] will have 9 rows and 5 columns.
Resulting D
  1.342196 -0.080722  0.942705  1.326401    0.0343
  0.952738  0.004302  0.680059  1.420024  0.749506
  -0.42315 -0.099356  0.602592 -0.128131  0.197995
 -0.874579 -1.989023 -1.479145 -1.963914 -1.480745
  0.403442  0.120605   0.38933  0.028234  0.202146
 -0.377875 -0.346078  0.662662  0.603232 -0.249229
   0.71807 -0.053794  0.956435  1.288518  0.184963
 -0.347932 -0.410419 -0.410582 -0.175535 -0.176471
 -1.163575 -0.363667 -0.779835  -0.91265 -0.809488

Is it full rank? true

Iteration 2
The matrix [E] = [A][B] will have 5 rows and 9 columns.
The matrix [H] will have 9 rows and 5 columns.
Resulting D
 -0.875651 -1.121959 -1.061803  -0.37018 -0.581922
 -0.021286  0.438606  0.439868 -0.348148 -0.243602
  0.115525  0.604465  0.305096 -0.068395 -0.584497
  0.755863  1.222269   0.56302  0.486334  0.292702
  0.203498  0.421592  0.160509  0.535417  0.587083
  0.779344  0.502896  0.338471  0.741886   0.28547
  1.077329  0.066448  0.018712  0.466584  0.771742
 -0.889879 -1.223156 -0.555143 -0.583784 -0.297454
  0.269486 -0.064799  0.055584  0.161476  0.126919

Is it full rank? true

Iteration 3
The matrix [E] = [A][B] will have 5 rows and 9 columns.
The matrix [H] will have 9 rows and 5 columns.
Resulting D
  0.201721 -0.647097  0.051457 -1.154906 -0.779558
 -0.623876 -1.614637 -0.886076 -1.052066 -1.868837
 -1.012378 -0.796751 -0.001973 -0.228281 -0.680911
 -0.211605  0.002338  0.524016  0.490409 -0.281095
  0.091177 -0.117529  0.522852  0.167912  0.201598
  0.263462  1.078214  0.049482  -0.31292  0.781295
 -0.484798  1.122343  0.285031  0.254435  1.012248
 -0.376376 -0.008645 -0.693459 -0.621088 -0.587739
  1.441024 -0.648415 -0.946848 -0.967546 -0.875984

Is it full rank? true

In addition it’s very much recommended to read the posts Linear Algebra Introduction.

Linear Algebra Introduction