Sparse and Special Structure Matrices

ojAlgo’s MatrixStore interface largely defines what you can do with matrices. There are many implementations of that interface provided. They store the elements differently and exploit structure to perform operations more efficiently. If you have special requirements, implementing your own MatrixStore is very easy.

Custom MatrixStore Example

CreateCustomMatrixStore.java
import static org.ojalgo.function.constant.PrimitiveMath.ONE;
import static org.ojalgo.function.constant.PrimitiveMath.ZERO;

import org.ojalgo.OjAlgoUtils;
import org.ojalgo.matrix.decomposition.QR;
import org.ojalgo.matrix.store.MatrixStore;
import org.ojalgo.matrix.store.PhysicalStore;
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.type.Stopwatch;

/**
 * Creating custom matrix stores exploiting special structures.
 *
 * @see https://www.ojalgo.org/2020/09/sparse-and-special-structure-matrices/
 * @see https://github.com/optimatika/ojAlgo/wiki/Create-Custom-MatrixStore
 */
public class CreateCustomMatrixStore {

    /**
     * In addition to the required 4 methods, you probably should also implement these.
     */
    static class BetterPerformingEye extends MostBasicEye {

        public BetterPerformingEye(final int rowsCount, final int columnsCount) {
            super(rowsCount, columnsCount);
        }

        /**
         * For a primitive valued implementation, at least, you should implement this method
         */
        @Override
        public double doubleValue(final int row, final int col) {
            return row == col ? ONE : ZERO;
        }

        /**
         * May allow the standard matrix multiplication code to exploit structure, certainly possible here.
         */
        @Override
        public int firstInColumn(final int col) {
            return Math.min(col, this.getRowDim());
        }

        /**
         * @see #firstInColumn(int)
         */
        @Override
        public int firstInRow(final int row) {
            return Math.min(row, this.getColDim());
        }

        /**
         * @see #firstInColumn(int)
         */
        @Override
        public int limitOfColumn(final int col) {
            return Math.min(col + 1, this.getRowDim());
        }

        /**
         * @see #firstInColumn(int)
         */
        @Override
        public int limitOfRow(final int row) {
            return Math.min(row + 1, this.getColDim());
        }

        /**
         * Custom/optimised copy functionality.
         */
        @Override
        public void supplyTo(final TransformableRegion<Double> receiver) {
            receiver.reset();
            receiver.fillDiagonal(ONE);
        }

    }

    /**
     * There are only 4 methods you have to implement in a custom MatrixStore, but there's a whole lot you can
     * to with that implementation. It's fully functional!
     */
    static class MostBasicEye implements MatrixStore<Double> {

        private final int myColumnsCount;
        private final int myRowsCount;

        public MostBasicEye(final int rowsCount, final int columnsCount) {
            super();
            myRowsCount = rowsCount;
            myColumnsCount = columnsCount;
        }

        @Override
        public Double get(final int row, final int col) {
            return row == col ? ONE : ZERO;
        }

        @Override
        public int getColDim() {
            return myColumnsCount;
        }

        @Override
        public int getRowDim() {
            return myRowsCount;
        }

        @Override
        public PhysicalStore.Factory<Double, R064Store> physical() {
            return R064Store.FACTORY;
        }

    }

    private static int DIM = 2000;

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

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

        MatrixStore<Double> mostBasicEye = new MostBasicEye(DIM + DIM, DIM);
        MatrixStore<Double> identityAndZero = R064Store.FACTORY.makeIdentity(DIM).below(DIM);
        MatrixStore<Double> betterPerformingEye = new BetterPerformingEye(DIM + DIM, DIM);

        Stopwatch stopwatch = new Stopwatch();

        BasicLogger.debug();
        BasicLogger.debug("Compare multiplication");

        PhysicalStore<Double> right = R064Store.FACTORY.makeFilled(DIM, DIM, new Normal());
        PhysicalStore<Double> product = R064Store.FACTORY.make(DIM + DIM, DIM);

        stopwatch.reset();
        mostBasicEye.multiply(right, product);
        BasicLogger.debug("MostBasicEye multiplied in {}", stopwatch.stop());

        stopwatch.reset();
        identityAndZero.multiply(right, product);
        BasicLogger.debug("Built-in ojAlgo multiplied in {}", stopwatch.stop());

        stopwatch.reset();
        betterPerformingEye.multiply(right, product);
        BasicLogger.debug("BetterPerformingEye multiplied in {}", stopwatch.stop());

        BasicLogger.debug();
        BasicLogger.debug("Compare QR decomposition");

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

        stopwatch.reset();
        decompQR.compute(mostBasicEye);
        BasicLogger.debug("MostBasicEye decomposed in {}", stopwatch.stop());

        stopwatch.reset();
        decompQR.compute(identityAndZero);
        BasicLogger.debug("Built-in ojAlgo decomposed in {}", stopwatch.stop());

        stopwatch.reset();
        decompQR.compute(betterPerformingEye);
        BasicLogger.debug("BetterPerformingEye decomposed in {}", stopwatch.stop());

        // Actually decomposing should take exactly the same time for the 3 alternatives.
        // ...and the input is already in the decomposed form.
        // The time difference is the time it takes to copy the input matrix
        // to the decomposer's internal storage.

    }

}

Here’s an example how to implement a simple eye matrix. From this it should be possible to deduce how to implement any other special structered matrix.

There are 2 implementations described in the example:

  1. MostBasicEye just implements the 4 methods that all MatrixStore implementations have to implement. As you can see these are trivial to implement, and this class is already fully functional.
  2. BetterPerformingEye extends MostBasicEye and adds a few methods that make it perform better. It’s just a handfull additional methods that in most cases would be trivial to implement but doing so can make a big difference.

To further improve things you may override any/all additional methods. This would be specific to your use case.

class CreateCustomMatrixStore
ojAlgo
2020-09-14


Compare multiplication
MostBasicEye multiplied in 6.852596385E9ns
Built-in ojAlgo multiplied in 8.480996E7ns
BetterPerformingEye multiplied in 5.2626592E7ns

Compare QR decomposition
MostBasicEye decomposed in 1.69935516E8ns
Built-in ojAlgo decomposed in 1.53854866E8ns
BetterPerformingEye decomposed in 2.1318605E7ns

Sparse Matrix Example

SparseMatrices.java
import static org.ojalgo.type.CalendarDateUnit.MILLIS;

import java.time.Instant;
import java.util.Random;

import org.ojalgo.OjAlgoUtils;
import org.ojalgo.array.ArrayR064;
import org.ojalgo.array.LongToNumberMap;
import org.ojalgo.array.SparseArray;
import org.ojalgo.matrix.store.MatrixStore;
import org.ojalgo.matrix.store.R064Store;
import org.ojalgo.matrix.store.SparseStore;
import org.ojalgo.matrix.task.iterative.ConjugateGradientSolver;
import org.ojalgo.netio.BasicLogger;
import org.ojalgo.series.BasicSeries;
import org.ojalgo.type.Stopwatch;

/**
 * Example use of SparseStore and other special MatrixStore implementations.
 *
 * @see https://www.ojalgo.org/2020/09/sparse-and-special-structure-matrices/
 * @see https://github.com/optimatika/ojAlgo/wiki/Sparse-Matrices
 */
public class SparseMatrices {

    private static String NON_ZEROS = "{} non-zeroes out of {} matrix elements calculated in {}";
    private static Random RANDOM = new Random();

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

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

        int dim = 100_000;

        SparseStore<Double> mtrxA = SparseStore.R064.make(dim, dim);
        SparseStore<Double> mtrxB = SparseStore.R064.make(dim, dim);
        SparseStore<Double> mtrxC = SparseStore.R064.make(dim, dim);
        MatrixStore<Double> mtrxZ = R064Store.FACTORY.makeZero(dim, dim);
        MatrixStore<Double> mtrxI = R064Store.FACTORY.makeIdentity(dim);

        // 5 matrices * 100k rows * 100k cols * 8 bytes per element => would be more than 372GB if dense
        // This program runs with default settings of any JVM

        for (int i = 0; i < dim; i++) {
            int j = RANDOM.nextInt(dim);
            double val = RANDOM.nextDouble();
            mtrxA.set(i, j, val);
        } // Each row of A contains 1 non-zero element at random column

        for (int j = 0; j < dim; j++) {
            int i = RANDOM.nextInt(dim);
            double val = RANDOM.nextDouble();
            mtrxB.set(i, j, val);
        } // Each column of B contains 1 non-zero element at random row

        Stopwatch stopwatch = new Stopwatch();

        BasicLogger.debug();
        BasicLogger.debug("Sparse-Sparse multiplication");
        stopwatch.reset();
        mtrxA.multiply(mtrxB, mtrxC);
        BasicLogger.debug(NON_ZEROS, mtrxC.nonzeros().estimateSize(), mtrxC.count(), stopwatch.stop(MILLIS));

        // It's the left matrix "this" that decides the multiplication algorithm,
        // and it knows nothing about the input/right matrix other than that it implements the required interface.
        // It could be another sparse matrix as in the example above. It could be a full/dense matrix. Or, it could something else...

        // Let's try an identity matrix...

        BasicLogger.debug();
        BasicLogger.debug("Sparse-Identity multiplication");
        stopwatch.reset();
        mtrxA.multiply(mtrxI, mtrxC);
        BasicLogger.debug(NON_ZEROS, mtrxC.nonzeros().estimateSize(), mtrxC.count(), stopwatch.stop(MILLIS));

        // ...or an all zeros matrix...

        BasicLogger.debug();
        BasicLogger.debug("Sparse-Zero multiplication");
        stopwatch.reset();
        mtrxA.multiply(mtrxZ, mtrxC);
        BasicLogger.debug(NON_ZEROS, mtrxC.nonzeros().estimateSize(), mtrxC.count(), stopwatch.stop(MILLIS));

        // What if we turn things around so that the identity or zero matrices become "this" (the left matrix)?

        BasicLogger.debug();
        BasicLogger.debug("Identity-Sparse multiplication");
        stopwatch.reset();
        mtrxI.multiply(mtrxB, mtrxC);
        BasicLogger.debug(NON_ZEROS, mtrxC.nonzeros().estimateSize(), mtrxC.count(), stopwatch.stop(MILLIS));

        BasicLogger.debug();
        BasicLogger.debug("Zero-Sparse multiplication");
        stopwatch.reset();
        mtrxZ.multiply(mtrxB, mtrxC);
        BasicLogger.debug(NON_ZEROS, mtrxC.nonzeros().estimateSize(), mtrxC.count(), stopwatch.stop(MILLIS));

        // Q: Why create identity and zero matrices?
        // A: The can be used as building blocks for larger logical structures.

        MatrixStore<Double> mtrxL = mtrxI.right(mtrxA).below(mtrxZ, mtrxB);

        // There's much more you can do with that logical builder...

        BasicLogger.debug();
        BasicLogger.debug("Scale logical structure");
        stopwatch.reset();
        MatrixStore<Double> mtrxScaled = mtrxL.multiply(3.14);
        BasicLogger.debug("{} x {} matrix scaled in {}", mtrxScaled.countRows(), mtrxScaled.countColumns(), stopwatch.stop(MILLIS));

        // By now we would have passed 1TB, if dense

        // Other data structures that are sparse, include:

        SparseArray<Double> sparseArray = SparseArray.factory(ArrayR064.FACTORY).make(dim);

        LongToNumberMap<Double> longToNumberMap = LongToNumberMap.factory(ArrayR064.FACTORY).make();
        BasicSeries<Instant, Double> series = BasicSeries.INSTANT.build(ArrayR064.FACTORY);
        ConjugateGradientSolver conjugateGradientSolver = new ConjugateGradientSolver();
    }

}

Among the provided MatrixStore implementations is a general purpose sparse matrix called SparseStore. Here’s a little example of what you can do with that, as well as a couple of the other special matrix types.

The example program below outputs various execution times, but the interesting part is not the speed. The main takeaway here is that the calculations are at all possible. If all the matrices created in the example had been dense implementations this program would require more than 1TB of memory to execute. Yet it runs with the default settings of a standard JVM (and finishes in a few seconds).

The execution times are significantly different for the various matrix multiplication alternatives in the example program. You’re encouraged to think about why/how that is, and then check the source code to verify your conclusions.

class SparseMatrices
ojAlgo
2020-09-14


Sparse-Sparse multiplication
100314 non-zeroes out of 10000000000 matrix elements calculated in 54.874324ms

Sparse-Identity multiplication
100000 non-zeroes out of 10000000000 matrix elements calculated in 27.532112ms

Sparse-Zero multiplication
0 non-zeroes out of 10000000000 matrix elements calculated in 14.865259ms

Identity-Sparse multiplication
100000 non-zeroes out of 10000000000 matrix elements calculated in 11.162651ms

Zero-Sparse multiplication
0 non-zeroes out of 10000000000 matrix elements calculated in 0.499874ms

Scale logical structure
200000 x 200000 matrix scaled in 30.004948ms