Sparse 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 just a little example of what it can do with sparse matrices, and a couple of the other special matrix types.

Example code

import java.time.Instant;
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.Stopwatch;

public class SparseMatrices {

    private static final String NON_ZEROES_OUT_OF_MATRIX_ELEMENTS_CALCULATED_IN = "{} 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_ZEROES_OUT_OF_MATRIX_ELEMENTS_CALCULATED_IN, mtrxC.nonzeros().estimateSize(), mtrxC.count(), stopwatch.stop());

        // 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_ZEROES_OUT_OF_MATRIX_ELEMENTS_CALCULATED_IN, mtrxC.nonzeros().estimateSize(), mtrxC.count(), stopwatch.stop());

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

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

        // 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_ZEROES_OUT_OF_MATRIX_ELEMENTS_CALCULATED_IN, mtrxC.nonzeros().estimateSize(), mtrxC.count(), stopwatch.stop());

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

        // 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());

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

        // SparseStore is backed by ojAlgo's SparseArray implementation
        final SparseArray<Double> array = SparseArray.factory(Primitive64Array.FACTORY, dim).make();

        // Other types that are backed by SparseArray include
        final LongToNumberMap<Double> map = LongToNumberMap.factory(Primitive64Array.FACTORY).make();
        final BasicSeries<Instant, Double> series = BasicSeries.INSTANT.build(Primitive64Array.FACTORY);
        final ConjugateGradientSolver solver = new ConjugateGradientSolver();
        // ...and others.

    }

}

Console output

SparseMatrices
ojAlgo
2017-07-20


Sparse-Sparse multiplication
99956 non-zeroes out of 10000000000 matrix elements calculated in 10625.42318MILLIS

Sparse-Identity multiplication
100000 non-zeroes out of 10000000000 matrix elements calculated in 3282.761063MILLIS

Sparse-Zero multiplication
0 non-zeroes out of 10000000000 matrix elements calculated in 31.930066MILLIS

Identity-Sparse multiplication
100000 non-zeroes out of 10000000000 matrix elements calculated in 18.62369MILLIS

Zero-Sparse multiplication
0 non-zeroes out of 10000000000 matrix elements calculated in 0.236331MILLIS

Scale logical structure
200000 x 200000 matrix scaled in 47.121008MILLIS