Sparse Matrices

An updated version of this page can be found at https://www.ojalgo.org/2020/09/sparse-and-special-structure-matrices/

---

ojAlgo is primarily known for its efficient multi-threaded algorthms on dense matrices, but ojAlgo also has sparse and other structured special matrix implementations. Here's just a little example of what it can do with sparse matrices, and 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.

Example code

import java.util.Random;

import org.ojalgo.OjAlgoUtils;
import org.ojalgo.array.LongToNumberMap;
import org.ojalgo.array.Primitive64Array;
import org.ojalgo.array.SparseArray;
import org.ojalgo.matrix.store.MatrixStore;
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.CalendarDateUnit;
import org.ojalgo.type.Stopwatch;

public class SparseMatrices {

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

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

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

        final int dim = 100_000;

        final SparseStore<Double> mtrxA = SparseStore.PRIMITIVE.make(dim, dim);
        final SparseStore<Double> mtrxB = SparseStore.PRIMITIVE.make(dim, dim);
        final SparseStore<Double> mtrxC = SparseStore.PRIMITIVE.make(dim, dim);
        final MatrixStore<Double> mtrxZ = MatrixStore.PRIMITIVE.makeZero(dim, dim).get();
        final MatrixStore<Double> mtrxI = MatrixStore.PRIMITIVE.makeIdentity(dim).get();

        // 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++) {
            final int j = RANDOM.nextInt(dim);
            final 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++) {
            final int i = RANDOM.nextInt(dim);
            final double val = RANDOM.nextDouble();
            mtrxB.set(i, j, val);
        } // Each column of B contains 1 non-zero element at random row

        final 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(CalendarDateUnit.MILLIS));

        // It's the left matrix 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(CalendarDateUnit.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(CalendarDateUnit.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(CalendarDateUnit.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(CalendarDateUnit.MILLIS));

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

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

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

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

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

        SparseArray.factory(Primitive64Array.FACTORY, dim).make();

        LongToNumberMap.factory(Primitive64Array.FACTORY).make();
        BasicSeries.INSTANT.build(Primitive64Array.FACTORY);
        new ConjugateGradientSolver();

    }

}

Console output

SparseMatrices
ojAlgo
2018-11-28


Sparse-Sparse multiplication
100229 non-zeroes out of 10000000000 matrix elements calculated in 74.31452ms

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

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

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

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

Scale logical structure
200000 x 200000 matrix scaled in 36.359936ms

To more precisely exploit special structures you may want to create a custom MatrixStore implemenation.