This document describes the comprehensive property-based test suite for the Cholesky factorization implementation in sparse-linear-algebra.
The Cholesky factorization decomposes a positive definite matrix A into a product L L^H, where L is a lower triangular matrix and L^H is the conjugate transpose of L. For real matrices, this becomes L L^T.
The original Cholesky implementation had two bugs for complex-valued matrices:
Changed from using **2 to proper complex conjugate multiplication:
-- Before (incorrect for complex numbers):
let l = sum (fmap (**2) lrow)
-- After (correct for both real and complex):
let l = sum (fmap (\x -> x * conj x) lrow)This ensures that for complex numbers, we compute x * conj(x) which gives the squared magnitude, rather than x^2 which would be incorrect.
Changed from using contractSub (which doesn't use conjugation) to using inner product:
-- Before (incorrect for complex numbers):
inn = contractSub ll ll i j (j - 1)
-- After (correct for both real and complex):
lrow_i = ifilterSV (\k _ -> k <= j - 1) (extractRow ll i)
lrow_j = ifilterSV (\k _ -> k <= j - 1) (extractRow ll j)
inn = lrow_i <.> lrow_jThe inner product operator <.> correctly uses conjugation for complex numbers: v <.> w = sum (v_i * conj(w_i)).
These fixes enable the Cholesky factorization to work correctly for both real symmetric positive definite matrices and complex Hermitian positive definite matrices.
The test suite consists of two main categories:
- 3x3 SPD matrix: A small dense SPD matrix for basic verification
- 5x5 tridiagonal SPD matrix: A sparse tridiagonal matrix to test sparsity handling
-
Reconstruction Property:
L L^T = A- Verifies that the factorization correctly reconstructs the original matrix
- Uses Frobenius norm to measure reconstruction error
-
Lower Triangular Property:
L is lower triangular- Ensures the factor L has no non-zero elements above the diagonal
-
Positive Diagonal Property: All diagonal elements of
Lare positive- Validates that diagonal entries
L[i,i] > 0for all i
- Validates that diagonal entries
- 2x2 HPD matrix: A small Hermitian matrix for basic verification
- 3x3 HPD matrix: A slightly larger Hermitian matrix
-
Reconstruction Property:
L L^H = A- Verifies that the factorization correctly reconstructs the original Hermitian matrix
- Uses Frobenius norm to measure reconstruction error
-
Lower Triangular Property:
L is lower triangular- Ensures the factor L has no non-zero elements above the diagonal
-
Positive Real Diagonal Property: All diagonal elements of
Lhave positive real parts- Validates that
realPart(L[i,i]) > 0for all i - For Hermitian matrices, the diagonal of L is typically real and positive
- Validates that
Generates random SPD matrices using the construction:
SPD = M^T M + εI
where M is a random matrix and ε = 0.1 is a small constant ensuring positive definiteness.
Generates random HPD matrices using the construction:
HPD = M^H M + εI
where M is a random complex matrix and ε = 0.1 + 0i.
The property-based tests use QuickCheck's sized combinator to generate matrices of varying dimensions (typically 2x2 to ~10x10), ensuring the factorization works correctly across different scales.
stack testOr with make:
make testThe test suite provides extensive coverage of:
- ✅ Real-valued symmetric positive definite matrices (sparse and dense)
- ✅ Complex-valued Hermitian positive definite matrices (dense)
- ✅ Various matrix sizes (2x2 through ~10x10)
- ✅ Mathematical properties (reconstruction, triangularity, positive diagonal)
- ✅ Edge cases via QuickCheck property-based testing
-
The Cholesky factorization requires the input matrix to be positive definite. The algorithm will throw a
NeedsPivotingexception if a zero diagonal element is encountered. -
For numerical stability, the test generators add a small diagonal perturbation (ε = 0.1) to ensure matrices are strictly positive definite.
-
All tests use the
nearZerofunction from theNumeric.Epsmodule to handle floating-point comparison tolerances appropriately.
- Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press.
- The Cholesky-Banachiewicz algorithm implemented proceeds row-by-row from the upper-left corner.