Linear Algebra Introduction

ojAlgo's linear algebra has two layers. The top one, MatrixR064 and its siblings for other number types, is designed to be easy to use. The matrices are immutable, and they choose, keep and reuse the matrix decompositions for you. The lower layer is the matrix stores and the decompositions themselves, where you are in full control of numerics and performance.

This post starts at the top, solving a system of equations in a few lines, and ends with how to decide which layer is right for you.

All the code is from one example program, listed in full at the end together with its output. For the Maven dependency, see Get Started.

Create a Matrix and Solve [A][x] = [b]

double[][] data = { { 4, -2, 1 }, { -2, 4, -2 }, { 1, -2, 4 } };

MatrixR064 A = MatrixR064.FACTORY.rows(data);
MatrixR064 b = MatrixR064.FACTORY.column(11.0, -16.0, 17.0);

MatrixR064 x = A.solve(b);

The solution is x = [1, -2, 3], and the residual |Ax - b| is 0.

  • Every matrix class has a static FACTORY. rows(...) creates a matrix from a double[][], one array per row.
  • There is no separate vector class. A vector is simply a matrix with a single column or a single row – no distinction is made. column(...) creates one with a single column, row(...) one with a single row.
  • solve chooses the decomposition for you. It tries a fast one first (typically LU) and, if that reports the matrix as singular, uses the singular value decomposition instead.
  • The matrix keeps the decomposition and reuses it where it can. Solving again with another right-hand side, or asking for the inverse, determinant or rank, doesn't start from scratch.

Ask the Matrix About Itself

BasicLogger.debug("Determinant: {}", A.getDeterminant());
BasicLogger.debug("Rank: {}", A.getRank());
BasicLogger.debug("Condition number: {}", A.getCondition());
BasicLogger.debug("Singular values: {}", A.getSingularValues());
BasicLogger.debug("Symmetric: {}", A.isSymmetric());

List<Eigenpair> eigenpairs = A.getEigenpairs();
for (Eigenpair pair : eigenpairs) {
    BasicLogger.debug("Eigenvalue: {}", pair.value);
}

BasicLogger.debugMatrix("Inverse", A.invert());
Determinant: 36.0
Rank: 3
Condition number: 4.529210992451759
Singular values: { 7.372281323269012, 3.000000000000001, 1.6277186767309857 }
Symmetric: true
Eigenvalue: (7.372281323269014 + 0.0i)
Eigenvalue: (3.0 + 0.0i)
Eigenvalue: (1.6277186767309861 + 0.0i)
Inverse
 0.333333 0.166667 0.000000
 0.166667 0.416667 0.166667
 0.000000 0.166667 0.333333

Eigenvalues are ComplexNumbers, because a real matrix in general has complex eigenvalues. This matrix is symmetric, so they are all real, and each Eigenpair also holds the corresponding eigenvector.

Singular values are real and non-negative for any matrix, also a non-square one. The condition number is the ratio of the largest to the smallest, 7.37 / 1.63 = 4.53, and the numerical rank is the number of singular values that are not negligible. For a symmetric positive definite matrix, like this one, the singular values and the eigenvalues are the same.

The inverse is there when you need it, but to solve equations call solve rather than multiplying by the inverse. It is faster and more accurate.

Matrices Are Immutable

MatrixR064 doubled = A.multiply(2.0);

MatrixR064.DenseReceiver receiver = A.copy();
receiver.set(0, 0, 5.0);
MatrixR064 modified = receiver.get();
A(0,0) = 4.0, doubled(0,0) = 8.0, modified(0,0) = 5.0

A MatrixR064 never changes. Every operation returns a new matrix, and instances can be shared between threads without any synchronisation. Immutability is also what makes reusing decompositions safe: once computed, a decomposition stays valid for the life of the matrix. To change elements, copy() the matrix to a receiver, set what you need, and call get() to get a new matrix. To build a matrix element by element from scratch, start with newDenseBuilder (or newSparseBuilder) on the factory.

A new matrix is not necessarily a new copy of the elements. doubled above has no elements of its own: it reads those of A and multiplies them by 2 when asked. transpose() likewise reads the original elements in transposed order, and chained element-wise functions, A.onAll(f).onAll(g) (see the next section), stack up as views of A without any intermediate matrices. Nothing is calculated until the elements are needed, and since A never changes that is safe. Operations that have to calculate a new result anyway, such as adding or multiplying two matrices, produce a matrix that holds that result.

Element-wise Functions

onAll applies a function to every element and, like every operation, returns a new matrix. The functions come from PrimitiveMath (statically imported here):

MatrixR064 halved = A.onAll(DIVIDE.by(2.0));
MatrixR064 squared = A.onAll(POWER, 2);
MatrixR064 roots = A.onAll(ABS.andThen(SQRT));

PrimitiveFunction.Unary logistic = v -> 1.0 / (1.0 + Math.exp(-v));
MatrixR064 squashed = A.onAll(logistic);
A(0,1) = -2.0, halved -1.0, squared 4.0, sqrt of abs 1.4142135623730951, logistic 0.11920292202211755
  • DIVIDE is a binary function. by(2.0) fixes its second argument and turns it into a unary function.
  • POWER takes an integer parameter, given directly to onAll.
  • Functions compose with andThen.
  • To write your own, assign a lambda to a PrimitiveFunction.Unary.

The other number types have the same: RationalMath and RationalFunction.Unary for MatrixQ128, ComplexMath and ComplexFunction.Unary for MatrixC128, and so on. The arithmetic is then done in that number type.

Least Squares

When there are more equations than unknowns, solve returns the least squares solution. There is no separate method to learn. Fitting a straight line y = c0 + c1t to five points:

double[] t = { 0.0, 1.0, 2.0, 3.0, 4.0 };
double[] y = { 1.1, 2.9, 5.2, 7.1, 8.8 };

MatrixR064.DenseReceiver builder = MatrixR064.FACTORY.newDenseBuilder(t.length, 2);
for (int i = 0; i < t.length; i++) {
    builder.set(i, 0, 1.0);
    builder.set(i, 1, t[i]);
}
MatrixR064 T = builder.get();

MatrixR064 c = T.solve(MatrixR064.FACTORY.column(y));
Least squares coefficients [c0, c1]
 1.100000 1.960000

Exact Arithmetic with Rational Numbers

MatrixR064 holds double elements. The same API exists for other number types:

Class Elements
MatrixR064 double
MatrixR032 float storage
MatrixR128 Quadruple, about twice the precision of double
MatrixQ128 RationalNumber, exact fractions
MatrixC128 ComplexNumber
MatrixH256 Quaternion

The Hilbert matrix, with elements 1/(i + j + 1), is a classic example of an ill-conditioned matrix. Its inverse has integer elements, but rounding errors grow quickly with its size. With MatrixQ128 the code is the same and the result is exact:

int n = 5;

MatrixQ128.DenseReceiver hilbertBuilder = MatrixQ128.FACTORY.newDenseBuilder(n, n);
for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
        hilbertBuilder.set(i, j, RationalNumber.of(1, i + j + 1));
    }
}
MatrixQ128 hilbert = hilbertBuilder.get();

MatrixQ128 exactInverse = hilbert.invert();
MatrixR064 roundedInverse = MatrixR064.FACTORY.copy(hilbert).invert();
Exact inverse of the 5x5 Hilbert matrix
      25    -300    1050   -1400     630
    -300    4800  -18900   26880  -12600
    1050  -18900   79380 -117600   56700
   -1400   26880 -117600  179200  -88200
     630  -12600   56700  -88200   44100
Relative error in the double precision inverse: 4.688658107551693E-10

Rational arithmetic is much slower than double, so it is for small problems where exactness matters.

Which Layer? Simplicity or Full Control

MatrixR064 makes choices for you. It picks the decomposition, falls back to another one when the first fails, caches and reuses what it computed, and creates new matrices rather than modifying existing ones. That is best effort, and for most uses it is exactly what you want.

If you need full control over numerics and performance, use the matrix stores and the decompositions directly. That is the decision: are you fine with MatrixR064 making these choices, or do you need to make them yourself?

MatrixR064 is built on that lower layer: the MatrixStore/PhysicalStore interfaces in org.ojalgo.matrix.store (R064Store is the dense double implementation) and the decompositions in org.ojalgo.matrix.decomposition. The same solve, with the decomposition chosen explicitly:

R064Store body = R064Store.FACTORY.copy(A);
R064Store rhs = R064Store.FACTORY.copy(b);

LU<Double> lu = LU.R064.make(body);

if (lu.compute(body)) {
    MatrixStore<Double> solution = lu.getSolution(rhs);
    BasicLogger.debugMatrix("Solution x, using R064Store and LU", solution);
}

compute returns false if the matrix is singular, and then it is up to you what to do. MatrixR064.solve would instead have moved on to the singular value decomposition.

This layer is for when you need to:

  • choose the decomposition, and decide what happens when it fails
  • get at the factors themselves: L and U, Q and R, singular vectors
  • work with mutable matrices and update them in place
  • reuse memory in loops that run many times, rather than create new matrices
  • exploit sparse or special structure

The two layers are interoperable. MatrixR064 and the stores implement the same interfaces – Access2D for reading the elements – so the factories and decompositions accept either. R064Store.FACTORY.copy(A) above works for that reason, lu.compute(A) would have worked just as well, and MatrixR064.FACTORY.copy(solution) takes you back up.

Still, work primarily in one of the layers. Don't move data back and forth between them, or to and from double[][], inside a loop. That is the Common Mistake.

Further reading:

Example Code

LinearAlgebraIntroduction.java
import static org.ojalgo.function.constant.PrimitiveMath.ABS;
import static org.ojalgo.function.constant.PrimitiveMath.DIVIDE;
import static org.ojalgo.function.constant.PrimitiveMath.POWER;
import static org.ojalgo.function.constant.PrimitiveMath.SQRT;

import java.util.List;

import org.ojalgo.OjAlgoUtils;
import org.ojalgo.function.PrimitiveFunction;
import org.ojalgo.matrix.MatrixQ128;
import org.ojalgo.matrix.MatrixR064;
import org.ojalgo.matrix.decomposition.Eigenvalue.Eigenpair;
import org.ojalgo.matrix.decomposition.LU;
import org.ojalgo.matrix.store.MatrixStore;
import org.ojalgo.matrix.store.R064Store;
import org.ojalgo.netio.BasicLogger;
import org.ojalgo.scalar.RationalNumber;

/**
 * Linear Algebra Introduction
 *
 * @see https://www.ojalgo.org/2026/10/linear-algebra-introduction/
 */
public class LinearAlgebraIntroduction {

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

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

        /* Create matrices and solve [A][x] = [b] */

        double[][] data = { { 4, -2, 1 }, { -2, 4, -2 }, { 1, -2, 4 } };

        MatrixR064 A = MatrixR064.FACTORY.rows(data);
        MatrixR064 b = MatrixR064.FACTORY.column(11.0, -16.0, 17.0);

        MatrixR064 x = A.solve(b);

        BasicLogger.debugMatrix("Solution x", x);
        BasicLogger.debug("Residual norm |Ax - b| = {}", A.multiply(x).subtract(b).norm());
        BasicLogger.debug();

        /* Ask the matrix about itself */

        BasicLogger.debug("Determinant: {}", A.getDeterminant());
        BasicLogger.debug("Rank: {}", A.getRank());
        BasicLogger.debug("Condition number: {}", A.getCondition());
        BasicLogger.debug("Singular values: {}", A.getSingularValues());
        BasicLogger.debug("Symmetric: {}", A.isSymmetric());

        List<Eigenpair> eigenpairs = A.getEigenpairs();
        for (Eigenpair pair : eigenpairs) {
            BasicLogger.debug("Eigenvalue: {}", pair.value);
        }

        BasicLogger.debugMatrix("Inverse", A.invert());
        BasicLogger.debug();

        /* Immutable: every operation returns a new matrix */

        MatrixR064 doubled = A.multiply(2.0);

        MatrixR064.DenseReceiver receiver = A.copy();
        receiver.set(0, 0, 5.0);
        MatrixR064 modified = receiver.get();

        BasicLogger.debug("A(0,0) = {}, doubled(0,0) = {}, modified(0,0) = {}", A.doubleValue(0, 0), doubled.doubleValue(0, 0), modified.doubleValue(0, 0));
        BasicLogger.debug();

        /* Element-wise functions: predefined, parameterised, composed or your own */

        MatrixR064 halved = A.onAll(DIVIDE.by(2.0));
        MatrixR064 squared = A.onAll(POWER, 2);
        MatrixR064 roots = A.onAll(ABS.andThen(SQRT));

        PrimitiveFunction.Unary logistic = v -> 1.0 / (1.0 + Math.exp(-v));
        MatrixR064 squashed = A.onAll(logistic);

        BasicLogger.debug("A(0,1) = {}, halved {}, squared {}, sqrt of abs {}, logistic {}", A.doubleValue(0, 1), halved.doubleValue(0, 1), squared.doubleValue(0, 1),
                roots.doubleValue(0, 1), squashed.doubleValue(0, 1));
        BasicLogger.debug();

        /* Least squares: fit a straight line y = c0 + c1 t */

        double[] t = { 0.0, 1.0, 2.0, 3.0, 4.0 };
        double[] y = { 1.1, 2.9, 5.2, 7.1, 8.8 };

        MatrixR064.DenseReceiver builder = MatrixR064.FACTORY.newDenseBuilder(t.length, 2);
        for (int i = 0; i < t.length; i++) {
            builder.set(i, 0, 1.0);
            builder.set(i, 1, t[i]);
        }
        MatrixR064 T = builder.get();

        MatrixR064 c = T.solve(MatrixR064.FACTORY.column(y));

        BasicLogger.debugMatrix("Least squares coefficients [c0, c1]", c.transpose());
        BasicLogger.debug();

        /* Exact arithmetic: the same API with rational numbers */

        int n = 5;

        MatrixQ128.DenseReceiver hilbertBuilder = MatrixQ128.FACTORY.newDenseBuilder(n, n);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                hilbertBuilder.set(i, j, RationalNumber.of(1, i + j + 1));
            }
        }
        MatrixQ128 hilbert = hilbertBuilder.get();

        MatrixQ128 exactInverse = hilbert.invert();
        MatrixR064 roundedInverse = MatrixR064.FACTORY.copy(hilbert).invert();

        MatrixR064 exactAsDouble = MatrixR064.FACTORY.copy(exactInverse);

        BasicLogger.debugMatrix("Exact inverse of the 5x5 Hilbert matrix", exactInverse);
        BasicLogger.debug("Relative error in the double precision inverse: {}", roundedInverse.subtract(exactAsDouble).norm() / exactAsDouble.norm());
        BasicLogger.debug();

        /* One level down: mutable stores and an explicitly chosen decomposition */

        R064Store body = R064Store.FACTORY.copy(A);
        R064Store rhs = R064Store.FACTORY.copy(b);

        LU<Double> lu = LU.R064.make(body);

        if (lu.compute(body)) {
            MatrixStore<Double> solution = lu.getSolution(rhs);
            BasicLogger.debugMatrix("Solution x, using R064Store and LU", solution);
        }
    }

}

Console Output

class LinearAlgebraIntroduction
ojAlgo
57.6

Solution x
  1.0
 -2.0
  3.0
Residual norm |Ax - b| = 0.0

Determinant: 36.0
Rank: 3
Condition number: 4.529210992451759
Singular values: { 7.372281323269012, 3.000000000000001, 1.6277186767309857 }
Symmetric: true
Eigenvalue: (7.372281323269014 + 0.0i)
Eigenvalue: (3.0 + 0.0i)
Eigenvalue: (1.6277186767309861 + 0.0i)
Inverse
 0.333333 0.166667 0.000000
 0.166667 0.416667 0.166667
 0.000000 0.166667 0.333333

A(0,0) = 4.0, doubled(0,0) = 8.0, modified(0,0) = 5.0

A(0,1) = -2.0, halved -1.0, squared 4.0, sqrt of abs 1.4142135623730951, logistic 0.11920292202211755

Least squares coefficients [c0, c1]
 1.100000 1.960000

Exact inverse of the 5x5 Hilbert matrix
      25    -300    1050   -1400     630
    -300    4800  -18900   26880  -12600
    1050  -18900   79380 -117600   56700
   -1400   26880 -117600  179200  -88200
     630  -12600   56700  -88200   44100
Relative error in the double precision inverse: 4.688658107551693E-10

Solution x, using R064Store and LU
  1.0
 -2.0
  3.0