Implementation:Haifengl Smile DenseMatrix Decomposition

Overview

The DenseMatrix Decomposition implementation provides methods on the DenseMatrix class for computing matrix factorizations: LU, QR, SVD, Eigenvalue (EVD), and Cholesky decompositions. Each method delegates to the corresponding LAPACK routine via the Java Foreign Function & Memory API. All decompositions are destructive -- they overwrite the calling matrix in-place. Use copy() before decomposing if the original is needed.

Type

API Doc

Source Location

base/src/main/java/smile/tensor/DenseMatrix.java:L990-1332

Import Statements

import smile.tensor.DenseMatrix;
import smile.tensor.ScalarType;
import smile.tensor.LU;
import smile.tensor.QR;
import smile.tensor.SVD;
import smile.tensor.EVD;
import smile.tensor.Cholesky;
import smile.linalg.SVDJob;
import smile.linalg.EVDJob;

External Dependencies

LAPACK Routine	Decomposition	Float64	Float32
LU factorization	`lu()`	`dgetrf_`	`sgetrf_`
QR factorization	`qr()`	`dgeqrf_`	`sgeqrf_`
SVD (divide-and-conquer)	`svd()`	`dgesdd_`	`sgesdd_`
Symmetric EVD	`eigen()` on symmetric	`dsyevd_`	`ssyevd_`
General EVD	`eigen()` on general	`dgeev_`	`sgeev_`
Cholesky factorization	`cholesky()`	`dpotrf_`	`spotrf_`

All routines are accessed via FFM through smile.linalg.lapack.clapack_h.

API Signatures

public abstract class DenseMatrix implements Matrix, Serializable {

    /**
     * LU decomposition. Overwrites this matrix.
     * @return LU record containing factored matrix, pivot vector, and info code.
     */
    public LU lu()

    /**
     * QR decomposition. Overwrites this matrix.
     * @return QR record containing factored matrix and tau vector.
     */
    public QR qr()

    /**
     * Compact SVD. Overwrites this matrix. Computes singular vectors.
     * @return SVD record with s, U, Vt.
     */
    public SVD svd()

    /**
     * SVD with optional singular vectors. Overwrites this matrix.
     * @param vectors whether to compute U and Vt.
     * @return SVD record.
     */
    public SVD svd(boolean vectors)

    /**
     * Eigenvalue decomposition. Overwrites this matrix.
     * Computes right eigenvectors only (no left).
     * @return EVD record.
     */
    public EVD eigen()

    /**
     * Eigenvalue decomposition with control over eigenvector computation.
     * Overwrites this matrix.
     * @param vl compute left eigenvectors.
     * @param vr compute right eigenvectors.
     * @return EVD record.
     */
    public EVD eigen(boolean vl, boolean vr)

    /**
     * Cholesky decomposition for SPD matrices. Overwrites this matrix.
     * Matrix must have UPLO set (be marked as symmetric).
     * @return Cholesky record.
     */
    public Cholesky cholesky()
}

Return Types

Method	Return Type	Record Components	Notes
`lu()`	`LU`	`DenseMatrix lu`, `int[] ipiv`, `int info`	`info > 0` means singular; `ipiv` uses 1-based LAPACK indexing
`qr()`	`QR`	`DenseMatrix qr`, `Vector tau`	`tau` holds Householder reflector scalars; use `Q()` and `R()` to extract factors
`svd()`	`SVD`	`int m`, `int n`, `Vector s`, `DenseMatrix U`, `DenseMatrix Vt`	Compact SVD: $U$ is $m \times k$ , $V^{T}$ is $k \times n$ , $k = \min (m, n)$
`svd(false)`	`SVD`	`int m`, `int n`, `Vector s`, `null`, `null`	Singular values only (faster)
`eigen()`	`EVD`	`Vector wr`, `Vector wi`, `DenseMatrix Vl`, `DenseMatrix Vr`	For symmetric: `wi` is null, `Vl == Vr`
`cholesky()`	`Cholesky`	`DenseMatrix lu`	Lower/upper triangular factor depending on UPLO

Implementation Details

LAPACK Workspace Query Pattern

All LAPACK routines in Smile follow a two-pass workspace query pattern:

Query pass: Call the routine with lwork = -1. LAPACK writes the optimal workspace size into work[0].
Compute pass: Allocate a work vector of the optimal size, then call the routine again.

// Example: QR workspace query (from source)
Vector work = vector(1);
int[] lwork = { -1 };
// Query optimal workspace
dgeqrf_(m_, n_, qr.memory, lda_, tau.memory, work.memory, lwork_, info_);

// Allocate optimal workspace
work = qr.vector((int) work.get(0));
lwork[0] = work.size();
// Compute QR
dgeqrf_(m_, n_, qr.memory, lda_, tau.memory, work.memory, lwork_, info_);

SVD Job Control

The svd() method uses SVDJob.COMPACT for the divide-and-conquer SVD (dgesdd):

COMPACT ('S'): Computes the first $\min (m, n)$ columns of $U$ and rows of $V^{T}$
NO_VECTORS ('N'): Computes only singular values (no $U$ or $V^{T}$ )

Symmetric EVD Dispatch

The eigen() method checks isSymmetric() to choose the algorithm:

Case	LAPACK Routine	Algorithm
Symmetric (`uplo != null`)	`dsyevd_` / `ssyevd_`	Divide-and-conquer (faster for large $n$ )
General	`dgeev_` / `sgeev_`	QR iteration with Hessenberg reduction

For the symmetric case, the matrix itself becomes the eigenvector matrix on output, and its uplo flag is cleared.

Cholesky Precondition

The cholesky() method requires uplo != null:

public Cholesky cholesky() {
    if (uplo == null) {
        throw new IllegalArgumentException("The matrix is not symmetric");
    }
    // ... calls dpotrf_
}

Usage Examples

LU Decomposition

import smile.tensor.DenseMatrix;
import smile.tensor.LU;
import smile.tensor.ScalarType;

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

// Copy before decomposing (lu() is destructive)
LU lu = A.copy().lu();

System.out.println("Singular: " + lu.isSingular());
System.out.println("Determinant: " + lu.det());
DenseMatrix inv = lu.inverse();

QR Decomposition

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

DenseMatrix A = DenseMatrix.randn(ScalarType.Float64, 100, 50);
QR qr = A.copy().qr();

// Extract Q and R explicitly
DenseMatrix Q = qr.Q();   // 100 x 50 orthogonal matrix
DenseMatrix R = qr.R();   // 50 x 50 upper triangular

// Verify orthogonality: Q' * Q should be close to identity
DenseMatrix QtQ = Q.ata();
System.out.println(QtQ.get(0, 0)); // ~1.0
System.out.println(QtQ.get(0, 1)); // ~0.0

SVD

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

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

// Full SVD with singular vectors
SVD svd = A.copy().svd();
System.out.println("Largest singular value: " + svd.s().get(0));
System.out.println("Matrix rank: " + svd.rank());
System.out.println("Condition number: " + svd.condition());

// Singular values only (faster, no U or Vt)
SVD svdNoVec = A.copy().svd(false);
System.out.println("Largest singular value: " + svdNoVec.s().get(0));

Eigenvalue Decomposition

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

// Symmetric matrix
double[][] symData = {{4, 2, 1}, {2, 5, 3}, {1, 3, 6}};
DenseMatrix S = DenseMatrix.of(symData).withUplo(UPLO.LOWER);
EVD evd = S.copy().eigen();

// All eigenvalues are real for symmetric matrices
System.out.println("Eigenvalue 0: " + evd.wr().get(0));
System.out.println("Eigenvectors: " + evd.Vr());

// Sort eigenvalues in descending order
EVD sorted = evd.sort();

// General (non-symmetric) matrix
double[][] genData = {{1, 2, 3}, {0, 4, 5}, {0, 0, 6}};
DenseMatrix G = DenseMatrix.of(genData);
EVD gevd = G.copy().eigen(true, true);  // both left and right eigenvectors

Cholesky Decomposition

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

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

Cholesky chol = A.copy().cholesky();
System.out.println("Determinant: " + chol.det());
System.out.println("Log-determinant: " + chol.logdet());

Error Handling

Condition	Exception	Method
LAPACK GETRF returns `info < 0`	`ArithmeticException`	`lu()`
LAPACK GEQRF returns `info != 0`	`ArithmeticException`	`qr()`
LAPACK GESDD returns `info != 0`	`ArithmeticException`	`svd()`
LAPACK SYEVD/GEEV returns `info != 0`	`ArithmeticException`	`eigen()`
LAPACK POTRF returns `info != 0`	`ArithmeticException`	`cholesky()` (matrix not positive definite)
`cholesky()` on non-symmetric matrix	`IllegalArgumentException`	`cholesky()`
`eigen()` on non-square matrix	`IllegalArgumentException`	`eigen()`
Unsupported ScalarType	`UnsupportedOperationException`	All decompositions

Performance Notes

Decomposition	Complexity	LAPACK Driver	Notes
LU	$\frac{2}{3} n^{3}$	`dgetrf`	Partial pivoting
QR	$\frac{4}{3} n^{2} m$	`dgeqrf`	Householder reflections, workspace query
SVD	$O (m n \cdot \min (m, n))$	`dgesdd`	Divide-and-conquer (faster than `dgesvd`)
EVD (symmetric)	$O (n^{3})$	`dsyevd`	Divide-and-conquer
EVD (general)	$O (n^{3})$	`dgeev`	QR iteration
Cholesky	$\frac{1}{3} n^{3}$	`dpotrf`	SPD only, 2x faster than LU

Environment

Heuristics

Heuristic:Haifengl_Smile_Platform_Native_Library_Selection

Page Connections

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

Principle

Implementation

Heuristic

Environment

Implementation:Haifengl Smile DenseMatrix Decomposition

Overview

Type

Source Location

Import Statements

External Dependencies

API Signatures

Return Types

Implementation Details

LAPACK Workspace Query Pattern

SVD Job Control

Symmetric EVD Dispatch

Cholesky Precondition

Usage Examples

LU Decomposition

QR Decomposition

SVD

Eigenvalue Decomposition

Cholesky Decomposition

Error Handling

Performance Notes

Related

Environment

Heuristics

Categories

Page Connections