Implementation:Haifengl Smile SVD EVD Results

Overview

The SVD/EVD Results implementation provides the accessor methods and utility operations on the SVD, EVD, and ARPACK types. These include extracting singular values and vectors, computing derived quantities (rank, condition number, null space, range, pseudoinverse), sorting eigenvalues, constructing diagonal matrices, and using ARPACK for large sparse eigenproblems.

Type

API Doc

Source Locations

Class	File	Key Methods
`SVD`	`base/src/main/java/smile/tensor/SVD.java`	`diag()`, `norm()`, `rank()`, `nullity()`, `condition()`, `range()`, `nullspace()`, `pinv()`
`EVD`	`base/src/main/java/smile/tensor/EVD.java`	`diag()`, `sort()`
`ARPACK`	`base/src/main/java/smile/tensor/ARPACK.java`	`syev()`, `eigen()`, `svd()`

Import Statements

import smile.tensor.SVD;
import smile.tensor.EVD;
import smile.tensor.ARPACK;
import smile.tensor.ARPACK.SymmOption;
import smile.tensor.ARPACK.AsymmOption;
import smile.tensor.DenseMatrix;
import smile.tensor.Vector;
import smile.tensor.Matrix;

External Dependencies

Library	Routines	Integration
ARPACK	`dsaupd_c`, `dseupd_c`, `ssaupd_c`, `sseupd_c` (symmetric)	FFM via `smile.linalg.arpack.arpack_h`
ARPACK	`dnaupd_c`, `dneupd_c`, `snaupd_c`, `sneupd_c` (asymmetric)	FFM via `smile.linalg.arpack.arpack_h`

API Signatures

SVD Record

public record SVD(int m, int n, Vector s, DenseMatrix U, DenseMatrix Vt) implements Serializable {
    /** Constructor for singular values only (no vectors). */
    public SVD(int m, int n, Vector s)

    /** Constructor inferring m, n from U and Vt. */
    public SVD(Vector s, DenseMatrix U, DenseMatrix Vt)

    /** Diagonal matrix of singular values (m x n). */
    public DenseMatrix diag()

    /** L2 matrix norm: largest singular value sigma_1. */
    public double norm()

    /** Effective numerical rank (number of significant singular values). */
    public int rank()

    /** Dimension of null space: min(m,n) - rank. */
    public int nullity()

    /** L2 condition number: sigma_1 / sigma_k. */
    public double condition()

    /** Orthonormal basis for range(A): first r columns of U. */
    public DenseMatrix range()

    /** Orthonormal basis for null(A): last nullity rows of Vt, transposed. */
    public DenseMatrix nullspace()

    /** Pseudo inverse: V * Sigma^{-1} * U'. */
    public DenseMatrix pinv()
}

EVD Record

public record EVD(Vector wr, Vector wi, DenseMatrix Vl, DenseMatrix Vr) implements Serializable {
    /** Constructor for symmetric EVD (wi is null, Vl == Vr). */
    public EVD(Vector w, DenseMatrix V)

    /**
     * Block diagonal eigenvalue matrix.
     * Diagonal: real parts. Sub/super-diagonal: imaginary parts.
     */
    public DenseMatrix diag()

    /**
     * Sort eigenvalues by descending magnitude.
     * Reorders eigenvectors accordingly.
     * @return a new sorted EVD.
     */
    public EVD sort()
}

ARPACK Interface

public interface ARPACK {

    /** Compute nev eigenvalues of a symmetric matrix. */
    static EVD syev(Matrix A, SymmOption which, int nev)
    static EVD syev(Matrix A, SymmOption which, int nev, int ncv, double tol)

    /** Compute nev eigenvalues of a general (asymmetric) matrix. */
    static EVD eigen(Matrix A, AsymmOption which, int nev)
    static EVD eigen(Matrix A, AsymmOption which, int nev, int ncv, double tol)

    /** Compute k largest singular triples via eigenproblem on A'A. */
    static SVD svd(Matrix A, int k)
    static SVD svd(Matrix A, int k, int ncv, double tol)

    /** Which eigenvalues of symmetric matrix to compute. */
    enum SymmOption { LA, SA, LM, SM, BE }

    /** Which eigenvalues of asymmetric matrix to compute. */
    enum AsymmOption { LM, SM, LR, SR, LI, SI }
}

Implementation Details

SVD Rank Computation

The rank is determined by counting singular values above a numerically derived threshold:

private double rcond() {
    return 0.5 * Math.sqrt(m + n + 1) * s.get(0) * MathEx.EPSILON;
}

public int rank() {
    int r = 0;
    double tol = rcond();
    for (int i = 0; i < s.size(); i++) {
        if (s.get(i) > tol) r++;
    }
    return r;
}

This uses machine epsilon ( $ϵ \approx 2.2 \times 1 0^{- 16}$ for Float64) scaled by the matrix dimensions and largest singular value.

SVD Pseudoinverse

The pseudoinverse $A^{+} = V Σ^{- 1} U^{T}$ inverts only the significant singular values:

public DenseMatrix pinv() {
    int k = s.size();
    Vector sigma = s.zeros(k);
    int r = rank();
    for (int i = 0; i < r; i++) {
        sigma.set(i, 1.0 / s.get(i));
    }
    // Compute V' * diag(1/sigma) * U' using the private adb() helper
    return adb(TRANSPOSE, Vt, sigma, TRANSPOSE, U);
}

EVD Sort Algorithm

The sort() method sorts eigenvalues by descending magnitude using QuickSort on the squared magnitudes:

public EVD sort() {
    int n = wr.size();
    double[] w = new double[n];
    if (wi != null) {
        for (int i = 0; i < n; i++) {
            w[i] = -(wr.get(i) * wr.get(i) + wi.get(i) * wi.get(i));
        }
    } else {
        for (int i = 0; i < n; i++) {
            w[i] = -(wr.get(i) * wr.get(i));
        }
    }
    int[] index = QuickSort.sort(w);
    // Reorder wr, wi, Vl, Vr according to index...
    return new EVD(wr2, wi2, Vl2, Vr2);
}

The negative sign ensures descending order when QuickSort sorts in ascending order.

EVD Diagonal Matrix

For the general case with complex eigenvalues, the diagonal matrix has a block structure:

public DenseMatrix diag() {
    DenseMatrix D = wr.diagflat();
    if (wi != null) {
        int n = wr.size();
        for (int i = 0; i < n; i++) {
            double wii = wi.get(i);
            if (wii > 0) {
                D.set(i, i + 1, wii);   // positive imaginary part above diagonal
            } else if (wii < 0) {
                D.set(i, i - 1, wii);   // negative imaginary part below diagonal
            }
        }
    }
    return D;
}

ARPACK Reverse Communication

ARPACK uses a reverse communication interface: the library repeatedly returns control to the caller with ido flags requesting a matrix-vector product, rather than calling a user-supplied function pointer.

do {
    // ARPACK computes the next Arnoldi step
    dsaupd_c(ido_, bmat_, n, bwhich_, nev, tol, resid.memory, ncv,
             V.memory, ldv, iparam_, ipntr_, workd.memory, workl.memory,
             workl.size(), info_);

    if (ido[0] == -1 || ido[0] == 1) {
        // ARPACK requests: compute A * x and store in y
        // ipntr[0]-1 and ipntr[1]-1 are offsets into workd
        A.mv(workd, ipntr[0] - 1, ipntr[1] - 1);
    }
} while (ido[0] == -1 || ido[0] == 1);

This design enables ARPACK to work with any matrix type (dense, sparse, or implicit) since it only needs the matrix-vector product operation.

ARPACK SVD via Eigenproblem

ARPACK computes the top- $k$ SVD by solving the symmetric eigenproblem $A^{T} A v = σ^{2} v$ using the AtA implicit matrix wrapper:

static SVD svd(Matrix A, int k, int ncv, double tol) {
    Matrix ata = new AtA(A);  // implicit A'A without forming it
    EVD eigen = syev(ata, SymmOption.LM, k, ncv, tol);

    Vector s = eigen.wr();
    for (int i = 0; i < s.size(); i++) {
        s.set(i, Math.sqrt(s.get(i)));  // eigenvalues of A'A -> singular values
    }

    // Recover U from V: u_j = A * v_j / sigma_j
    // (or V from U if m < n)
    // ...
}

Usage Examples

Extracting SVD Components

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();

// Singular values (Vector)
Vector sigma = svd.s();
System.out.println("Top 5 singular values:");
for (int i = 0; i < 5; i++) {
    System.out.printf("  sigma[%d] = %.6f%n", i, sigma.get(i));
}

// Diagonal matrix of singular values (100 x 50)
DenseMatrix S = svd.diag();

// Matrix properties
System.out.println("L2 norm: " + svd.norm());
System.out.println("Rank: " + svd.rank());
System.out.println("Nullity: " + svd.nullity());
System.out.println("Condition number: " + svd.condition());

// Range and null space
DenseMatrix R = svd.range();
System.out.println("Range space dimension: " + R.ncol());

DenseMatrix N = svd.nullspace();
if (N != null) {
    System.out.println("Null space dimension: " + N.ncol());
}

// Pseudoinverse
DenseMatrix Apinv = svd.pinv();
System.out.println("Pseudo-inverse size: " + Apinv.nrow() + " x " + Apinv.ncol());

Working with EVD Results

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

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

// Sort eigenvalues in descending order
EVD sorted = evd.sort();
System.out.println("Eigenvalues (descending):");
for (int i = 0; i < sorted.wr().size(); i++) {
    System.out.printf("  lambda[%d] = %.6f%n", i, sorted.wr().get(i));
}

// Access eigenvectors
DenseMatrix V = sorted.Vr();
System.out.println("Eigenvector for largest eigenvalue: " + V.column(0));

// Reconstruct diagonal eigenvalue matrix
DenseMatrix Lambda = sorted.diag();

// General non-symmetric EVD
double[][] genData = {{0, 1}, {-2, -3}};
DenseMatrix G = DenseMatrix.of(genData);
EVD gevd = G.copy().eigen(false, true);  // right eigenvectors only

System.out.println("Real parts: " + gevd.wr());
if (gevd.wi() != null) {
    System.out.println("Imaginary parts: " + gevd.wi());
}

// Block diagonal matrix with complex eigenvalue structure
DenseMatrix D = gevd.diag();

ARPACK for Large Sparse Eigenproblems

import smile.tensor.ARPACK;
import smile.tensor.ARPACK.SymmOption;
import smile.tensor.DenseMatrix;
import smile.tensor.EVD;
import smile.tensor.SVD;
import smile.tensor.ScalarType;
import smile.linalg.UPLO;

// Large symmetric matrix (e.g., graph Laplacian)
DenseMatrix L = DenseMatrix.randn(ScalarType.Float64, 500, 500);
DenseMatrix A = L.ata();  // 500x500 SPD matrix

// Compute top 10 eigenvalues (largest algebraic)
EVD evd = ARPACK.syev(A, SymmOption.LA, 10);
System.out.println("Top 10 eigenvalues:");
for (int i = 0; i < 10; i++) {
    System.out.printf("  lambda[%d] = %.6f%n", i, evd.wr().get(i));
}

// Compute smallest 5 eigenvalues
EVD evdSmall = ARPACK.syev(A, SymmOption.SA, 5);

// With custom parameters: ncv=30 Arnoldi vectors, tol=1e-10
EVD evdCustom = ARPACK.syev(A, SymmOption.LM, 10, 30, 1e-10);

ARPACK SVD for Large Matrices

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

// Large matrix
DenseMatrix A = DenseMatrix.randn(ScalarType.Float64, 1000, 500);

// Compute top 20 singular triples
SVD svd = ARPACK.svd(A, 20);
System.out.println("Top 20 singular values:");
for (int i = 0; i < 20; i++) {
    System.out.printf("  sigma[%d] = %.6f%n", i, svd.s().get(i));
}

// Access the truncated U (1000 x 20) and Vt (20 x 500)
DenseMatrix U = svd.U();
DenseMatrix Vt = svd.Vt();
System.out.println("U: " + U.nrow() + " x " + U.ncol());
System.out.println("Vt: " + Vt.nrow() + " x " + Vt.ncol());

ARPACK for Asymmetric Eigenproblems

import smile.tensor.ARPACK;
import smile.tensor.ARPACK.AsymmOption;
import smile.tensor.DenseMatrix;
import smile.tensor.EVD;
import smile.tensor.ScalarType;

// Non-symmetric matrix
DenseMatrix A = DenseMatrix.randn(ScalarType.Float64, 200, 200);

// Compute 5 eigenvalues with largest magnitude
EVD evd = ARPACK.eigen(A, AsymmOption.LM, 5);
System.out.println("Eigenvalues (real part): " + evd.wr());
System.out.println("Eigenvalues (imaginary part): " + evd.wi());

ARPACK Parameters Reference

Parameter	Description	Default
`nev`	Number of eigenvalues to compute	Required, $0 < nev < n$
`ncv`	Number of Arnoldi vectors (subspace dimension)	$\min (3 \cdot nev, n)$
`tol`	Convergence tolerance	$1 0^{- 6}$
`maxIter`	Maximum iterations (set internally)	$10 \cdot n$
`bmat`	Problem type	`'I'` (standard eigenproblem)

Error Handling

Condition	Exception	Source
`rank()` or `condition()` on partial SVD	`UnsupportedOperationException`	`SVD`
`range()` or `nullspace()` without singular vectors	`IllegalStateException`	`SVD`
`pinv()` or `solve()` without singular vectors	`IllegalStateException`	`SVD`
Non-square matrix for ARPACK	`IllegalArgumentException`	`ARPACK`
Invalid `nev` ( $\leq 0$ or $\geq n$ )	`IllegalArgumentException`	`ARPACK`
ARPACK convergence failure	`ArithmeticException`	`ARPACK.syev()`, `ARPACK.eigen()`
Fewer eigenvalues computed than requested	`ArithmeticException`	`ARPACK`

Performance Comparison: Full vs. ARPACK

Method	Complexity	Memory	Use Case
`DenseMatrix.eigen()`	$O (n^{3})$	$O (n^{2})$	All eigenvalues, small to medium $n$
`ARPACK.syev()`	$O (n \cdot k \cdot ncv)$	$O (n \cdot ncv)$	$k$ eigenvalues, large $n$
`DenseMatrix.svd()`	$O (m n \cdot \min (m, n))$	$O (m n)$	All singular values, small to medium
`ARPACK.svd()`	$O (nnz \cdot k \cdot ncv)$	$O (n \cdot ncv)$	Top- $k$ singular values, large matrices

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 SVD EVD Results

Overview

Type

Source Locations

Import Statements

External Dependencies

API Signatures

SVD Record

EVD Record

ARPACK Interface

Implementation Details

SVD Rank Computation

SVD Pseudoinverse

EVD Sort Algorithm

EVD Diagonal Matrix

ARPACK Reverse Communication

ARPACK SVD via Eigenproblem

Usage Examples

Extracting SVD Components

Working with EVD Results

ARPACK for Large Sparse Eigenproblems

ARPACK SVD for Large Matrices

ARPACK for Asymmetric Eigenproblems

ARPACK Parameters Reference

Error Handling

Performance Comparison: Full vs. ARPACK

Related

Environment

Heuristics

Categories

Page Connections