Implementation:Haifengl Smile Decomposition Solvers

Overview

The Decomposition Solvers implementation provides solve() methods on each decomposition record type (LU, QR, Cholesky, SVD) as well as the iterative BiconjugateGradient solver. Each direct solver delegates to the appropriate LAPACK back-solve routine, while the BiCG solver uses only matrix-vector products (mv/tv) and is compatible with any Matrix implementation (dense, sparse, or implicit).

Type

API Doc

Source Locations

Class	File	Lines
`LU`	`base/src/main/java/smile/tensor/LU.java`	L97-166
`QR`	`base/src/main/java/smile/tensor/QR.java`	L115-215
`Cholesky`	`base/src/main/java/smile/tensor/Cholesky.java`	L105-162
`SVD`	`base/src/main/java/smile/tensor/SVD.java`	L254-299
`BiconjugateGradient`	`base/src/main/java/smile/tensor/BiconjugateGradient.java`	L47-179

Import Statements

import smile.tensor.DenseMatrix;
import smile.tensor.Vector;
import smile.tensor.Matrix;
import smile.tensor.LU;
import smile.tensor.QR;
import smile.tensor.SVD;
import smile.tensor.Cholesky;
import smile.tensor.BiconjugateGradient;
import smile.tensor.Preconditioner;

API Signatures

LU Solver

public record LU(DenseMatrix lu, int[] ipiv, int info) implements Serializable {
    /** Solve A * x = b. Returns solution vector. */
    public Vector solve(double[] b)
    public Vector solve(float[] b)

    /** Solve A * X = B. B is overwritten with solution X. */
    public void solve(DenseMatrix B)

    /** Returns the matrix inverse. */
    public DenseMatrix inverse()

    /** Returns the matrix determinant. */
    public double det()

    /** Returns true if the matrix is singular. */
    public boolean isSingular()
}

QR Solver

public record QR(DenseMatrix qr, Vector tau) implements Serializable {
    /** Solve least squares min ||B - A*X||. Returns solution vector. */
    public Vector solve(double[] b)
    public Vector solve(float[] b)

    /** Solve least squares min ||B - A*X||. B is overwritten. */
    public void solve(DenseMatrix B)

    /** Extract the orthogonal factor Q. */
    public DenseMatrix Q()

    /** Extract the upper triangular factor R. */
    public DenseMatrix R()

    /** Convert to Cholesky decomposition of A'A. */
    public Cholesky toCholesky()
}

Cholesky Solver

public record Cholesky(DenseMatrix lu) implements Serializable {
    /** Solve A * x = b. Returns solution vector. */
    public Vector solve(double[] b)
    public Vector solve(float[] b)

    /** Solve A * X = B. B is overwritten with solution X. */
    public void solve(DenseMatrix B)

    /** Returns the matrix inverse. */
    public DenseMatrix inverse()

    /** Returns the matrix determinant. */
    public double det()

    /** Returns the log of matrix determinant. */
    public double logdet()
}

SVD Solver

public record SVD(int m, int n, Vector s, DenseMatrix U, DenseMatrix Vt) implements Serializable {
    /** Solve least squares min ||B - A*X||. Requires singular vectors. */
    public Vector solve(double[] b)
    public Vector solve(float[] b)
    public Vector solve(Vector b)

    /** Returns the pseudo inverse V * Sigma^{-1} * U'. */
    public DenseMatrix pinv()
}

Biconjugate Gradient Solver

public interface BiconjugateGradient {
    /** Solve A * x = b with Jacobi preconditioner. */
    static double solve(Matrix A, Vector b, Vector x)

    /** Solve A * x = b with custom preconditioner and tolerance. */
    static double solve(Matrix A, Vector b, Vector x,
                        Preconditioner P, double tol, int itol, int maxIter)
}

Implementation Details

LU Solve: LAPACK GETRS

The LU solver calls dgetrs_/sgetrs_ which solves $A x = b$ using the precomputed LU factorization with pivoting:

// From LU.solve(DenseMatrix B) - simplified
byte[] trans = { NO_TRANSPOSE.lapack() };
switch(lu.scalarType()) {
    case Float64 -> dgetrs_(trans_, n_, nrhs_, lu.memory, lda_, ipiv_, B.memory, ldb_, info_);
    case Float32 -> sgetrs_(trans_, n_, nrhs_, lu.memory, lda_, ipiv_, B.memory, ldb_, info_);
}

The solve(double[] b) method creates a new vector, copies b into it, and calls the in-place solver. It does not overwrite the input array.

QR Solve: ORMQR + TRTRS

The QR solver performs two LAPACK calls:

dormqr_/sormqr_: Multiplies $B$ by $Q^{T}$ from the left using the Householder representation
dtrtrs_/strtrs_: Solves the upper triangular system $R x = Q^{T} b$

The result vector is sliced to the first $n$ elements (for overdetermined systems where $m > n$ ):

public Vector solve(double[] b) {
    Vector x = qr.vector(qr.m);
    for (int i = 0; i < qr.m; i++) x.set(i, b[i]);
    solve(x);
    return x.slice(0, qr.n);  // return only the first n elements
}

SVD Solve: Manual Pseudoinverse Application

The SVD solver does not call a LAPACK solve routine. Instead it manually computes $x = V_{r} Σ_{r}^{- 1} U_{r}^{T} b$ where the subscript $r$ denotes restriction to the first $r$ significant singular values:

public Vector solve(Vector b) {
    int r = rank();
    DenseMatrix Ur = r == U.ncol() ? U : U.submatrix(0, 0, m, r);
    DenseMatrix Vr = r == Vt.nrow() ? Vt : Vt.submatrix(0, 0, r, n);
    Vector Utb = Ur.vector(r);
    Ur.tv(b, Utb);              // U_r' * b
    for (int i = 0; i < r; i++) {
        Utb.set(i, Utb.get(i) / s.get(i));  // Sigma_r^{-1} * U_r' * b
    }
    return Vr.tv(Utb);          // V_r * Sigma_r^{-1} * U_r' * b
}

Cholesky Solve: LAPACK POTRS

The Cholesky solver calls dpotrs_/spotrs_:

switch(lu.scalarType()) {
    case Float64 -> dpotrs_(uplo_, n_, nrhs_, lu.memory, lda_, B.memory, ldb_, info_);
    case Float32 -> spotrs_(uplo_, n_, nrhs_, lu.memory, lda_, B.memory, ldb_, info_);
}

BiCG: Iterative Krylov Method

The BiCG solver uses only Matrix.mv() and Matrix.tv() operations, making it compatible with any matrix type:

static double solve(Matrix A, Vector b, Vector x, Preconditioner P,
                    double tol, int itol, int maxIter) {
    // Initialize residual r = b - A*x
    A.mv(x, r);
    for (j = 0; j < n; j++) {
        double br = b.get(j) - r.get(j);
        r.set(j, br);
        rr.set(j, br);
    }

    for (int iter = 1; iter <= maxIter; iter++) {
        P.solve(rr, zz);         // Preconditioner solve
        // ... compute search directions ...
        A.mv(p, z);              // Matrix-vector product
        A.tv(pp, zz);            // Transpose matrix-vector product
        // ... update solution and residual ...
        if (err <= tol) break;
    }
    return err;
}

The default preconditioner is Jacobi (diagonal scaling): Preconditioner.Jacobi(A).

Usage Examples

Solving with LU

import smile.tensor.DenseMatrix;
import smile.tensor.LU;
import smile.tensor.Vector;

double[][] data = {{2, 1, 1}, {4, 3, 3}, {8, 7, 9}};
DenseMatrix A = DenseMatrix.of(data);
double[] b = {1, 1, 1};

LU lu = A.copy().lu();
if (!lu.isSingular()) {
    Vector x = lu.solve(b);
    System.out.println("Solution: " + x);
    System.out.println("Determinant: " + lu.det());
}

// Multiple right-hand sides
DenseMatrix B = DenseMatrix.of(new double[][]{{1, 0}, {1, 0}, {1, 0}});
lu.solve(B);  // B is overwritten with solution X

Solving Least Squares with QR

import smile.tensor.DenseMatrix;
import smile.tensor.QR;
import smile.tensor.ScalarType;

// Overdetermined system: 100 equations, 3 unknowns
DenseMatrix A = DenseMatrix.randn(ScalarType.Float64, 100, 3);
double[] b = new double[100];
for (int i = 0; i < 100; i++) b[i] = 2.0 * i + 1.0;

QR qr = A.copy().qr();
Vector x = qr.solve(b);
System.out.println("Least squares solution: " + x);
System.out.println("x has " + x.size() + " elements"); // 3

Solving SPD System with Cholesky

import smile.tensor.DenseMatrix;
import smile.tensor.Cholesky;
import smile.tensor.ScalarType;
import smile.linalg.UPLO;

// Construct SPD matrix: A = B'B + I (ensure positive definite)
DenseMatrix B = DenseMatrix.randn(ScalarType.Float64, 50, 10);
DenseMatrix A = B.ata();  // 10x10, SPD, marked symmetric

Cholesky chol = A.copy().cholesky();
double[] b = new double[10];
for (int i = 0; i < 10; i++) b[i] = 1.0;

Vector x = chol.solve(b);
System.out.println("Cholesky solution: " + x);

Solving Rank-Deficient System with SVD

import smile.tensor.DenseMatrix;
import smile.tensor.SVD;
import smile.tensor.ScalarType;

DenseMatrix A = DenseMatrix.randn(ScalarType.Float64, 100, 50);
SVD svd = A.copy().svd();

double[] b = new double[100];
for (int i = 0; i < 100; i++) b[i] = Math.sin(i);

Vector x = svd.solve(b);
System.out.println("SVD solution (" + x.size() + " elements): " + x);
System.out.println("Effective rank: " + svd.rank());

// Pseudo inverse for multiple solves
DenseMatrix pinv = svd.pinv();

Iterative Solving with BiCG

import smile.tensor.Matrix;
import smile.tensor.Vector;
import smile.tensor.DenseMatrix;
import smile.tensor.BiconjugateGradient;
import smile.tensor.ScalarType;

// For large sparse systems
DenseMatrix A = DenseMatrix.randn(ScalarType.Float64, 100, 100);
// Make it diagonally dominant for convergence
for (int i = 0; i < 100; i++) A.set(i, i, A.get(i, i) + 100.0);

Vector b = Vector.column(new double[100]);
for (int i = 0; i < 100; i++) b.set(i, 1.0);
Vector x = Vector.column(new double[100]);  // initial guess: zeros

double err = BiconjugateGradient.solve(A, b, x);
System.out.println("BiCG estimated error: " + err);

Error Handling

Condition	Exception	Solver
LU singular matrix	`RuntimeException`	`LU.solve()`
Dimension mismatch (A vs B)	`IllegalArgumentException`	All solvers
ScalarType mismatch	`IllegalArgumentException`	`LU.solve()`, `QR.solve()`, `Cholesky.solve()`
Layout mismatch	`IllegalArgumentException`	`LU.solve()`
Singular vectors not available	`IllegalStateException`	`SVD.solve()`
LAPACK GETRS/POTRS/TRTRS error	`ArithmeticException`	Direct solvers
Invalid tolerance ( $\leq 0$ )	`IllegalArgumentException`	`BiconjugateGradient.solve()`
Invalid itol (not 1--4)	`IllegalArgumentException`	`BiconjugateGradient.solve()`

Performance Comparison

Solver	Factorization Cost	Solve Cost (single RHS)	Solve Cost (k RHS)
LU	$\frac{2}{3} n^{3}$	$2 n^{2}$	$2 k n^{2}$
Cholesky	$\frac{1}{3} n^{3}$	$n^{2}$	$k n^{2}$
QR	$\frac{4}{3} n^{2} m$	$2 m n + n^{2}$	$k (2 m n + n^{2})$
SVD	$O (m n^{2})$	$O (m n)$	$O (k m n)$
BiCG	None	$O (iter \cdot nnz)$	N/A (solve independently)

Environment

Heuristics

Page Connections

Double-click a node to navigate. Hold to expand connections.

Principle

Implementation

Heuristic

Environment

Implementation:Haifengl Smile Decomposition Solvers

Overview

Type

Source Locations

Import Statements

API Signatures

LU Solver

QR Solver

Cholesky Solver

SVD Solver

Biconjugate Gradient Solver

Implementation Details

LU Solve: LAPACK GETRS

QR Solve: ORMQR + TRTRS

SVD Solve: Manual Pseudoinverse Application

Cholesky Solve: LAPACK POTRS

BiCG: Iterative Krylov Method

Usage Examples

Solving with LU

Solving Least Squares with QR

Solving SPD System with Cholesky

Solving Rank-Deficient System with SVD

Iterative Solving with BiCG

Error Handling

Performance Comparison

Related

Environment

Heuristics

Categories

Page Connections