Matrix Decompositions

1. summary How to perform common matrix decompositions in EJML

Introduction

Matrix decomposition are used to reduce a matrix to a more simplic format which can be easily solved and used to extract characteristics from. Below is a list of matrix decompositions and data structures there are implementations for.

Decomposition	DMatrixRMaj	DMatrixRBlock	ZMatrixRMaj
LU	Yes		Yes
Cholesky L`*`LT and RT`*`R	Yes	Yes	Yes
Cholesky L`*`D`*`LT	Yes
QR	Yes	Yes	Yes
QR Column Pivot	Yes
Singular Value Decomposition (SVD)	Yes
Generalized Eigen Value	Yes
Symmetric Eigen Value	Yes	Yes
Bidiagonal	Yes
Tridiagonal	Yes	Yes	Yes
Hessenberg	Yes		Yes

Accessing Pivots and Other Internal Structures

Most of the time you don't need access to internal data structures for a decomposition, just the results of the decomposition. If you need access to information, like the row or column pivots, then you will want to use a decomposition interface. Those can be created in a DecompositionFactory. For example, LUDecomposition provides access to its row pivots.

Solving Using Decompositions

Decompositions, such as LU and QR, are used to solve a linear system. A common mistake in EJML is to directly decompose the matrix instead of using a LinearSolver. LinearSolvers simplify the process of solving a linear system, are very fast, and will automatically be updated as new algorithms are added. It is recommended that you use them whenever possible.

For more information on LinearSolvers see the wikipage at Solving Linear Systems.

SimpleMatrix

SimpleMatrix has easy to an use interface built in for SVD and EVD:

SimpleSVD svd = A.svd();
SimpleEVD evd = A.eig();

SimpleMatrix U = svd.getU();

where A is a SimpleMatrix.

As with most operators in SimpleMatrix, it is possible to chain a decompositions with other commands. For instance, to print the singular values in a matrix:

A.svd().getW().extractDiag().transpose().print();

Other decompositions can be performed by using accessing the internal DMatrixRMaj and using the decompositions shown in the following section below. The following is an example of applying a Cholesky decomposition.

CholeskyDecomposition_F64<DMatrixRMaj> chol = DecompositionFactory_DDRM.chol(A.numRows(),true);
    
if( !chol.decompose(A.getMatrix()))
   throw new RuntimeException("Cholesky failed!");
        
SimpleMatrix L = SimpleMatrix.wrap(chol.getT(null));

DecompositionFactory

The best way to create a matrix decomposition is by using DecompositionFactory_DDRM. Directly instantiating a decomposition is discouraged because of the added complexity. DecompositionFactory_DDRM is updated as new and faster algorithms are added.

public interface DecompositionInterface<T extends Matrix> {
    /**
     * Computes the decomposition of the input matrix.  Depending on the implementation
     * the input matrix might be stored internally or modified.  If it is modified then
     * the function {@link #inputModified()} will return true and the matrix should not be
     * modified until the decomposition is no longer needed.
     *
     * @param orig The matrix which is being decomposed.  Modification is implementation dependent.
     * @return Returns if it was able to decompose the matrix.
     */
    public boolean decompose( T orig );

    /**
     * Is the input matrix to {@link #decompose(org.ejml.data.DMatrixRMaj)} is modified during
     * the decomposition process.
     *
     * @return true if the input matrix to decompose() is modified.
     */
    public boolean inputModified();
}

Most decompositions in EJML implement DecompositionInterface. To decompose "A" matrix simply call decompose(A). It returns true if there are no error while decomposing and false otherwise. While in general you can trust the results if true is returned some algorithms can have faults that are not reported. This is true for all linear algebra libraries.

To minimize memory usage, most decompositions will modify the original matrix passed into decompose(). Call inputModified() to determine if the input matrix is modified or not. If it is modified, and you do not wish it to be modified, just pass in a copy of the original instead.

Below is an example of how to compute the SVD of a matrix:

    void decompositionExample( DenseMatrix64F A ) {
        SingularValueDecomposition_F64<DMatrixRMaj> svd = DecompositionFactory_DDRM.svd(A.numRows,A.numCols);

        if( !svd.decompose(A) )
            throw new RuntimeException("Decomposition failed");

        DMatrixRMaj U = svd.getU(null,false);
        DMatrixRMaj W = svd.getW(null);
        DMatrixRMaj V = svd.getV(null,false);

    }

Note how it checks the returned value from decompose.

In addition DecompositionFactory_DDRM provides functions for computing the quality of a decomposition. Being able measure the decomposition's quality is an important way to validate its correctness. It works by "reconstructing" the original matrix then computing the difference between the reconstruction and the original. Smaller the quality is the better the decomposition is. With an ideal value of being around 1e-15 in most cases.

if( DecompositionFactory_DDRM.quality(A,svd) > 1e-3 )
    throw new RuntimeException("Bad decomposition");

List of functions in DecompositionFactory_DDRM

Decomposition	Code
LU	DecompositionFactory_DDRM.lu()
QR	DecompositionFactory_DDRM.qr()
QRP	DecompositionFactory_DDRM.qrp()
Cholesky	DecompositionFactory_DDRM.chol()
Cholesky LDL	DecompositionFactory_DDRM.cholLDL()
SVD	DecompositionFactory_DDRM.svd()
Eigen	DecompositionFactory_DDRM.eig()

Helper Functions for SVD and Eigen

Two classes SingularOps_DDRM and EigenOps_DDRM have been provided for extracting useful information from these decompositions or to provide highly specialized ways of computing the decompositions. Below is a list of more common uses of these functions:

SingularOps_DDRM

• descendingOrder()
  • In EJML the ordering of the returned singular values is not in general guaranteed. This function will reorder the U,W,V matrices such that the singular values are in the standard largest to smallest ordering.
• nullSpace()
  • Computes the null space from the provided decomposition.
• rank()
  • Returns the matrix's rank.
• nullity()
  • Returns the matrix's nullity.

EigenOps_DDRM

• computeEigenValue()
  • Given an eigen vector compute its eigenvalue.
• computeEigenVector()
  • Given an eigenvalue compute its eigenvector.
• boundLargestEigenValue()
  • Returns a lower and upper bound for the largest eigenvalue.
• createMatrixD() and createMatrixV()
  • Reformats the results such that two matrices (D and V) contain the eigenvalues and eigenvectors respectively. This is similar to the format used by other libraries such as Jama.