Solving Linear Systems

A fundamental problem in linear algebra is solving systems of linear equations. A linear system is any equation than can be expressed in this format:

A*x = b

where A is m by n, x is n by o, and b is m by o. Most of the time o=1. The best way to solve these equations depends on the structure of the matrix A. For example, if it's square and positive definite then Cholesky decomposition is the way to go. On the other hand if it is tall m > n, then QR is the way to go.

Each of the three interfaces (Procedural, SimpleMatrix, Equations) provides high level ways to solve linear systems which don't require you to specify the underlying algorithm. While convenient, these are not always the best approach in high performance situations. They create/destroy memory and don't provide you with access to their full functionality. If the best performance is needed then you should use a LinearSolver or one of its derived interfaces for a specific family of algorithms.

First a description is provided on how to solve linear systems using Procedural, SimpleMatrix, and then Equations. After that an overview of LinearSolver is presented.

High Level Interfaces

All high level interfaces essentially use the same code at the low level, which is the Procedural interface. This means that they have the same strengths and weaknesses. Their strength is simplicity. They will automatically select LU and QR decomposition, depending on the matrix's shape.

You should use the lower level LinearSolver if any of the following are true:

• Your matrix can some times be singular
• You wish to perform a pseudo inverse
• You need to avoid creating new memory
• You need to select a specific decomposition
• You need access to the low level decomposition

The case of singular or nearly singular matrices is worth discussing more. All of these high level approaches do attempt to detect singular matrices. The problem is that they aren't reliable and no tweaking of thresholds will make them reliable. If you are in a situation where you need to come up with a solution and it might be singular then you really need to know what you are doing. If a system is singular it means there are an infinite number of solutions.

Procedural

The way to solve linear systems in the Procedural interface is with CommonOps.solve(). Make sure you check it's return value to see if it failed! It might fail if the matrix is singular or nearly singular.

if( !CommonOps_DDRM.solve(A,b,x) ) {
  throw new IllegalArgument("Singular matrix");
}

SimpleMatrix

try {
  SimpleMatrix x = A.solve(b);
} catch ( SingularMatrixException e ) {
  throw new IllegalArgument("Singular matrix");
}

SingularMatrixException is a RuntimeException and you technically don't have to catch it. If you don't catch it, it will take down your whole application if the matrix is singular!

Equations

eq.process("x=b/A");

If it's singular it will throw a RuntimeException.

Low level Linear Solvers

Low level linear solvers in EJML all implement the LinearSolver interface. It provides a lot more power than the high level interfaces but is also more difficult to use and require more diligence. For example, you can no longer assume that it won't modify the input matrices!

LinearSolver

The LinearSolver interface is designed to be easy to use and to provide most of the power that comes from directly using a decomposition would provide.

public interface LinearSolver< T extends Matrix> {

    public boolean setA( T A );

    public T getA();

    public double quality();

    public void solve( T B , T X );

    public void invert( T A_inv );

    public boolean modifiesA();

    public boolean modifiesB();

    public <D extends DecompositionInterface>D getDecomposition();
}

Each linear solver implementation is built around a different decomposition. The best way to create a new LinearSolver instance is with LinearSolverFactory_DDRM. It provides an easy way to select the correct solver without plowing through the documentation.

Two steps are required to solve a system with a LinearSolver, as is shown below:

LinearSolver<DMatrixRMaj> solver = LinearSolverFactory_DDRM.qr(A.numRows,A.numCols);
if( !solver.setA(A) ) {
  throw new IllegalArgument("Singular matrix");
}

if( solver.quality() <= 1e-8 )
  throw new IllegalArgument("Nearly singular matrix");

solver.solve(b,x);

As with the high-level interfaces you can't trust algorithms such as QR, LU, or Cholesky to detect singular matrices! Sometimes they will work and sometimes they will not. Even adjusting the quality threshold won't fix the problem in all situations.

Additional capabilities included in LinearSolver are:

• invert()
  • Will invert a matrix more efficiently than solve() can.
• quality()
  • Returns a positive number which if it is small indicates a singular or nearly singular system system. Much faster to compute than the SVD.
• modifiesA() and modifiesB()
  • To reduce memory requirements, most LinearSolvers will modify the 'A' and store the decomposition inside of it. Some do the same for 'B' These function tell the user if the inputs are being modified or not.
• getDecomposition()
  • Provides access to the internal decomposition used.

LinearSolverSafe

If the input matrices 'A' and 'B' should not be modified then LinearSolverSafe is a convenient way to ensure that precondition:

LinearSolver<DMatrixRMaj> solver = LinearSolverFactory_DDRM.leastSquares();
solver = new LinearSolverSafe<DMatrixRMaj>(solver);

Pseudo Inverse

EJML provides two different pseudo inverses. One is SVD based and the other QRP. QRP stands for QR with column pivots. QRP can be thought of as a light weight SVD, much faster to compute but doesn't handle singular matrices quite as well.

LinearSolver<DMatrixRMaj> pinv = LinearSolverFactory_DDRM.pseudoInverse(true);

This will create an SVD based pseudo inverse. Otherwise if you specify false then it will create a QRP pseudo-inverse.

AdjustableLinearSolver

In situations where rows from the linear system are added or removed (see Example Polynomial Fitting) an AdjustableLinearSolver_DDRM can be used to efficiently resolve the modified system. AdjustableLinearSolver_DDRM is an extension of LinearSolver that adds addRowToA() and removeRowFromA(). These functions add and remove rows from A respectively. After being involved the solution can be recomputed by calling solve() again.

AdjustableLinearSolver_DDRM solver = LinearSolverFactory_DDRM.adjustable();

if( !solver.setA(A) ) {
  throw new IllegalArgument("Singular matrix");
}

solver.solve(b,x);

// add a row
double row[] = new double[N];

... code ...

solver.addRowToA(row,2);

.... adjust b and x ....

solver.solve(b,x); 

// remove a row
solver.removeRowFromA(7);

.... adjust b and x ....

solver.solve(b,x);