Linear algebra examples

Introduction

The purpose of this section is to give you tips on how to use the linear algebra package LAPACK provided by Finastra.

First, you must consider the context where you want to use this library. If you want to perform several linear algebra’s operations on a huge range of small matrices in parallel, this library is not suitable and you should use the APACHE library. The LAPACK package is designed to perform linear algebra tasks on very large matrices, with each task taking advantage of the parallelism. The goal is to run algorithms as fast as possible.

The main available features are the following ones:

In order to use the functions, you have to declare and initialize a Lapack object. The Lapack interface provides a set of linear algebra functions. The getLapack() method provides a default implementation of these functions.

Here is an example:

final Lapack lapack = getLapack();

Matrix Multiplication

Matrix multiplication is one of the most important functions in linear algebra. That is why it is very important to have the most efficient algorithm, since it is one of the main pillars of the linear algebra. The Matrix multiplication algorithm suits very well to the parallelism. The use of the matrix multiplication is the easiest of all. Here is an example:

final Lapack lapack = getLapack();
 
final double[][] A = new double[size][size];
final double[][] B = new double[size][size];
 
// Do your stuff with A and B
 
final double[][] C = lapack.multiplyMatrices(A,B); // Performs the matrix multiplication of A and B

Cholesky Decomposition

In order to be able to use the Cholesky Decomposition, you have to use an Hermitian, symmetric, positive-definite matrix. Since the matrix is symmetric, an upper or lower Cholesky Decomposition can be done, using the upper or the lower triangular part of the matrix. The choice can be made with the use of a flag. You can use it as follows:

final Lapack lapack = getLapack();
 
final double[][] m = new double[size][size];
// Do your stuff with m
 
final CholeskyDecompositionResult resultLower = lapack.choleskyDecomposition(m, true); // Performs a lower Cholesky Decomposition
final CholeskyDecompositionResult resultUpper = lapack.choleskyDecomposition(m, false); // Performs an upper Cholesky Decomposition

The result of a Cholesky decomposition is a CholeskyDecompositionResult object. The class has two private variables, as you can see:

public class CholeskyDecompositionResult {
    // Constructor and getters...   
    private final double[][] res; // Stores the result of the decomposition.
    private final boolean isLower; // Says if the decomposition is lower or upper.
}

The LAPACK provides a way to solve a linear system using a Cholesky decomposition. You can choose to solve the linear system with a CholeskyDecompositionResult object (upper or lower), or directly with a matrix if you do not care about the decomposition result. In this case, you have to specify if you want to use a lower or upper decomposition. Here is a simple example to solve linear system using a Cholesky decomposition with LAPACK:

final Lapack lapack = getLapack();
 
final double[][] A = new double[size][size];
final double[][] b = new double[size][linearSystemSize];
 
// Do your stuff with A and b
 
final CholeskyDecompositionResult resultUpper = lapack.choleskyDecomposition(A, false);
final double[][] x = lapack.solveCholesky(resultUpper , b); // Solves the linear system with an upper Cholesky Decomposition
final double[][] xBis = lapack.solveCholesky(A, true, b); // Solves the linear system with a lower Cholesky Decomposition

The upper Cholesky decomposition is slightly faster than the lower Cholesky decomposition.

LU Decomposition

LU decomposition factors a matrix as the product of two matrices, a lower triangular matrix L and an upper triangular matrix U. Since the algorithm implemented uses permutations, the decomposition of \(A\) with an LU decomposition is equal to \(P x L x U\), with \(P\) a permutation matrix. Here is an example:

final Lapack lapack = getLapack();
 
final double[][] m = new double[size][size];
// Do your stuff with m
 
final LUDecompositionResult result = lapack.luDecomposition(m); // Performs the LU Decomposition

The result of an LU Decomposition is an LUDecompositionResult. The class has three private variables, as you can see:

public class LUDecompositionResult {
    // Constructor and getters...
    private final double[][] l; // The lower triangular matrix.
    private final double[][] u; // The upper triangular matrix.
    private final int[] pivot; // The permutation matrix stored in the vector form (pivot[4] = 10 says that there is a permutation between the fourth and the tenth columns)
}

LAPACK provides a way to solve a linear system using an LU decomposition. You can choose to solve a system using the LU Decomposition with an LUDecompositionResult object or directly with a matrix. Here is a simple example to solve linear system using an LU Decomposition with LAPACK:

final Lapack lapack = getLapack();
 
final double[][] A = new double[size][size];
final double[][] b = new double[size][linearSystemSize];
 
// Do your stuff with A and b
 
final LUDecompositionResult result = lapack.luDecomposition(A);
final double[][] x = lapack.solveLU(result , b); // Solves the linear system with an LU Decomposition
final double[][] xBis = lapack.solveLU(A, b); // Solves the linear system with LU Decomposition without storing the result of the decomposition

LAPACK also provides a way to compute the determinant of a matrix given its LU decomposition:

final Lapack lapack = getLapack();
 
final double[][] A = new double[size][size];
// Do your stuff with A
final LUDecompositionResult result = lapack.luDecomposition(A);
final double determinant = lapack.getDeterminant(result); // Computes the determinant of the matrix A

QR Decomposition

The QR Decomposition factors a matrix as the product of two matrices, an orthogonal matrix \(Q\) and an upper triangular matrix \(R\). You can choose to compute a full QR decomposition, or only a reduced QR Decomposition. In the case of the reduced QR decomposition, the result is a matrix \(D\) and a vector \(\tau\). The elements on and above the diagonal of the matrix \(D\) contain the upper triangular matrix \(R\). The elements below the diagonal of \(D\), with the vector \(\tau\), represent the orthogonal matrix \(Q\) as a product of elementary reflectors. Here is an example:

final Lapack lapack = getLapack();
 
final double[][] m = new double[size][size];
// Do your stuff with m
 
final QRDecompositionResult result = lapack.qrDecomposition(m); // Performs a full QR Decomposition
final ReduceQRDecompositionResult reduceResult = lapack.reduceQrDecomposition(m); // Performs a reduce QR Decomposition

The result of a full QR Decomposition is a QRDecompositionResult. The class has two private variables, as you can see:

public class QRDecompositionResult{
    // Constructor and getters...
    private final double[][] q; // The orthogonal matrix Q
    private final double[][] r; // The upper triangular matrix R
}

The result of a reduced QR Decomposition is a ReduceQRDecompositionResult. The class has two private variables, as you can see:

public class ReduceQRDecompositionResult{
    // Constructor and getters...
    /*
     * Description from Netlib's LAPACK
     * The elements on and above the diagonal of the array contain the upper trapezoidal matrix R.
     * The elements below the diagonal, with the array TAU, represent the orthogonal matrix Q as a product of elementary reflectors
     * The matrix Q is represented as a product of elementary reflectors
     *
     * Q = H(1) H(2) . . . H(k), where k = min(m,n).
     *
     * Each H(i) has the form
     *
     * H(i) = I - tau * v * v'
     *
     * where tau is a real scalar, and v is a real vector with
     * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
     * and tau in TAU(i).
     */
    private final double[][] decomp; // The matrix D
    private final double[] tau; // The vector Tau
}

LAPACK provides a way to solve a linear system using a QR decomposition. You can choose to solve a system using the QR Decomposition with a QRDecompositionResult object, a ReduceQRDecompositionResult object or directly with a matrix. In this case, an implicit reduced QR decomposition is performed before solving the linear system:

final Lapack lapack = getLapack();
 
final double[][] A = new double[size][size];
final double[][] b = new double[size][linearSystemSize];
 
// Do your stuff with A and b
 
final QRDecompositionResult result = lapack.qrDecomposition(A); 
final ReduceQRDecompositionResult reduceResult = lapack.reduceQrDecomposition(A); 
final double[][] x1 = lapack.solveQR(result, b); // Solves the linear system with a full QR decomposition
final double[][] x2 = lapack.solveQR(reduceResult, b); // Solves the linear system with a reduce QR decomposition
final double[][] x3 = lapack.solveQR(A, b); // Solves the linear system with a reduce QR Decomposition without storing the result of the decomposition

LAPACK also provides a way to compute the determinant of a matrix given its reduced QR decomposition:

final Lapack lapack = getLapack();
 
final double[][] A = new double[size][size];
// Do your stuff with A
final ReduceQRDecompositionResult reduceResult = lapack.reduceQrDecomposition(A);
final double determinant = lapack.getDeterminant(reduceResult); // Computes the determinant of the matrix A
final double determinant2 = lapack.getDeterminant(A); // Computes the determinant of the matrix A with a reduce QR Decomposition without storing the result of the decomposition

Since the reduced QR decomposition does not compute \(Q\), it is twice faster than the full QR decomposition (there is the same number of operation to compute \(R\) and \(Q\), but \(Q\) have to be computed after \(R\) is finished). So if you do not not need the full QR decomposition, it is highly recommended to use the reduced QR Decomposition for the system solving. Even if the system solving using a reduced QR Decomposition is slightly slower than the system solving with full QR Decomposition, it is still a lot faster if you add the execution time of the decomposition.