Cynopolis/Vector3D
1
0
Fork 0
Code Issues 3 Pull Requests Actions Packages Projects Releases Wiki Activity

Templatize SVD on N: support any Matrix<R,C> with stack-only buffers

Browse Source 
Replace the fixed 5x5 SVD implementation with template <uint8_t N>
building blocks over Matrix<N,N> working buffers (N = max(rows, cols)),
removing the 5x5 size limit. Public API (SVD::SVD) is unchanged and the
whole path stays heap-free: peak stack is ~11*N^2 floats, budgeted in
the SVD.hpp header doc.

Fixes found while porting/validating the templated rewrite:
- QL.Identity() was a no-op (static factory returns by value); init
  the Householder accumulators with an explicit diagonal loop
- restore the QL column sign-flip in ExtractAndSortSingularValues for
  negative unsolved W diagonal entries (Householder sign flips)
- T = B^T B tridiagonal formula: T[i][i] = d[i]^2 + e[i-1]^2 only
  (e[i] contributes to T[i+1][i+1], not T[i][i])
- wide-matrix Vt assembly: Vt = QL^T must be filled over the full
  m x m (m = columns of A), not just the top n x n

Tests:
- matrix-tests: add large-size instantiation cases beyond the old limit
  (tall 7x5 N=7, square 6x6 N=6, wide 5x8 N=8 transpose path with full
  orthogonal 8x8 Vt, tall 6x4 near rank-deficient N=6 deflation path),
  all checked against numpy/scipy float32 references
- svd-build-blocks-tests: adapt Jacobi test to the Matrix<N,N> interface
- svd-reference-values.py: add the four new reference matrices
This commit is contained in:
Cynopolis
2026-08-20 10:22:49 -04:00
parent f8221dd9db
commit ab0cea104c
5 changed files with 737 additions and 459 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/SVD.cpp
+418 -379
View File
File diff suppressed because it is too large Load Diff
src/SVD.hpp
+139 -78
View File
@@ -3,6 +3,21 @@
/** /**
* @brief library that uses Matrix.hpp and performs SVD on a matrix * @brief library that uses Matrix.hpp and performs SVD on a matrix
*
* @note Fully templated: SVD works for ANY Matrix<R, C> with R, C in
* 1..255 (the uint8_t range of Matrix). There is no 5×5 limit.
*
* @note EMBEDDED CONSTRAINT — no heap. All working storage is stack
* allocated as templated Matrix<N,N> buffers where
* N = max(R, C). Peak stack usage per SVD call is
* ≈ 11·N² floats (≈ 44·N² bytes):
* N = 5 → ~1.1 KB
* N = 10 → ~4.4 KB
* N = 20 → ~18 KB
* N = 50 → ~110 KB
* N = 100 → ~440 KB
* N = 255 → ~2.9 MB
* Instantiate only the sizes that fit your call-stack budget.
*/ */
namespace SVD { namespace SVD {
/** /**
@@ -18,14 +33,28 @@ namespace SVD {
* The decomposition satisfies: A ≈ U × diag(Σ) × Vᵀ * The decomposition satisfies: A ≈ U × diag(Σ) × Vᵀ
* Singular values are returned in descending order. * Singular values are returned in descending order.
* *
* @param U Output: left singular vectors (m×k orthogonal matrix) * Output storage conventions:
* @param sigma Output: singular values as k×1 column vector, sorted descending * - U: Matrix<rows, columns> — first k columns are meaningful
* @param Vt Output: right singular vectors transposed (k×n matrix) * (rows k..columns−1 are zero in the wide case)
* - sigma: Matrix<columns, 1> — first k entries are the singular
* values; entries beyond k (wide matrices only) are zero
* - Vt: Matrix<columns, columns> — first k rows are meaningful
* (zero-padded in the tall case)
* *
* @note This implementation uses the Golub-Kahan-Reinsch algorithm: * For wide matrices (rows < columns) the SVD is computed on Aᵀ and the
* 1. Householder bidiagonalization of A * factors are swapped back.
* 2. Implicit QR iteration on the bidiagonal matrix *
* 3. Accumulation of U and V factors throughout * @tparam rows Number of rows in A (1..255)
* @tparam columns Number of columns in A (1..255)
* @param matrixToDecompose Input: the matrix A
* @param U Output: left singular vectors (rows×columns matrix)
* @param sigma Output: singular values (columns×1 vector, sorted descending)
* @param Vt Output: right singular vectors transposed (columns×columns)
*
* @note This implementation uses Householder bidiagonalization followed
* by block reduction: 2×2 blocks via closed form, larger blocks
* via cyclic Jacobi eigen-decomposition of BᵀB with residual
* singular values σᵢ = ‖B·vᵢ‖ (see docs/svd-refactor.md).
*/ */
template <uint8_t rows, uint8_t columns> template <uint8_t rows, uint8_t columns>
void SVD(Matrix<rows, columns> &matrixToDecompose, Matrix<rows, columns> &U, void SVD(Matrix<rows, columns> &matrixToDecompose, Matrix<rows, columns> &U,
@@ -33,7 +62,14 @@ void SVD(Matrix<rows, columns> &matrixToDecompose, Matrix<rows, columns> &U,
// ======================================================================== // ========================================================================
// SVD Building Block Functions (for unit testing) // SVD Building Block Functions (for unit testing)
// These operate on internal 5×5 working arrays for maximum flexibility. //
// Templated on the working-buffer size N. All block operations work on
// N×N matrices with runtime bounds (m, n, p, blockSize, ...) — the
// regions beyond the bounds are zero-padded working space.
//
// N is deduced from the Matrix arguments at the call site, e.g.
// Matrix<8, 8> W, QL, QR;
// SVD::Bidiagonalize(W, 6, 8, 6, QL, QR); // N = 8 deduced
// ======================================================================== // ========================================================================
/** /**
@@ -42,9 +78,9 @@ void SVD(Matrix<rows, columns> &matrixToDecompose, Matrix<rows, columns> &U,
* Given input vector x, computes normalized v and scalar alpha such that: * Given input vector x, computes normalized v and scalar alpha such that:
* (I - 2·v·vᵀ) · x = [alpha, 0, 0, ...]ᵀ * (I - 2·v·vᵀ) · x = [alpha, 0, 0, ...]ᵀ
* *
* @param x Input vector (up to 5 elements) * @param x Input vector (up to len elements)
* @param len Number of valid elements in x * @param len Number of valid elements in x
* @param v Output: normalized Householder vector (v[0] is the first element) * @param v Output: normalized Householder vector (length ≥ len)
* @param alpha Output: the resulting first element after reflection * @param alpha Output: the resulting first element after reflection
* @return The norm of the input vector x * @return The norm of the input vector x
*/ */
@@ -54,36 +90,43 @@ static float ComputeHouseholder(const float *x, uint8_t len, float *v,
/** /**
* @brief Apply a Householder reflection from the left. * @brief Apply a Householder reflection from the left.
* *
* Transforms W = (I - 2·v·vᵀ) · W where v operates on rows [startRow..endRow]. * Transforms W = (I - 2·v·vᵀ) · W where v operates on rows [startRow..endRow]
* and is applied across all N columns (zero-padded columns are a no-op).
* *
* @param W Input/output: matrix to transform (5×5 working array) * @tparam N Working buffer size
* @param W Input/output: matrix to transform
* @param v Householder vector (length = endRow - startRow + 1) * @param v Householder vector (length = endRow - startRow + 1)
* @param startRow First row index * @param startRow First row index
* @param endRow Last row index * @param endRow Last row index
*/ */
static void ApplyHouseholderLeft(Matrix<5, 5> &W, const float *v, template <uint8_t N>
static void ApplyHouseholderLeft(Matrix<N, N> &W, const float *v,
uint8_t startRow, uint8_t endRow); uint8_t startRow, uint8_t endRow);
/** /**
* @brief Apply a Householder reflection from the right. * @brief Apply a Householder reflection from the right.
* *
* Transforms W = W · (I - 2·v·vᵀ) where v operates on columns * Transforms W = W · (I - 2·v·vᵀ) where v operates on columns
* [startCol..endCol]. * [startCol..endCol] and is applied across all N rows (zero-padded rows
* are a no-op).
* *
* @param W Input/output: matrix to transform (5×5 working array) * @tparam N Working buffer size
* @param W Input/output: matrix to transform
* @param v Householder vector (length = endCol - startCol + 1) * @param v Householder vector (length = endCol - startCol + 1)
* @param startCol First column index * @param startCol First column index
* @param endCol Last column index * @param endCol Last column index
*/ */
static void ApplyHouseholderRight(Matrix<5, 5> &W, const float *v, template <uint8_t N>
static void ApplyHouseholderRight(Matrix<N, N> &W, const float *v,
uint8_t startCol, uint8_t endCol); uint8_t startCol, uint8_t endCol);
/** /**
* @brief Reduce a matrix to upper bidiagonal form using Householder reflections. * @brief Reduce a matrix to upper bidiagonal form using Householder reflections.
* *
* Applies a sequence of Householder reflections to reduce the input matrix * Applies a sequence of Householder reflections to reduce the input matrix
* W (m×q, where q ≥ p) to upper bidiagonal form B (p×q), accumulating * W (m×q, where q ≥ p, stored in N×N working space) to upper bidiagonal
* the left and right transformation matrices in QL and QR respectively. * form B (p×q), accumulating the left and right transformation matrices
* in QL and QR respectively.
* *
* Algorithm (Golub-Kahan bidiagonalization): * Algorithm (Golub-Kahan bidiagonalization):
* For k = 0 to p-1: * For k = 0 to p-1:
@@ -93,34 +136,34 @@ static void ApplyHouseholderRight(Matrix<5, 5> &W, const float *v,
* The accumulated transformations satisfy: * The accumulated transformations satisfy:
* QLᵀ · W_original · QR = B (upper bidiagonal) * QLᵀ · W_original · QR = B (upper bidiagonal)
* *
* @param W Input/output: matrix to bidiagonalize (5×5, must be at least p×q) * @tparam N Working buffer size (≥ m and ≥ q)
* @param W Input/output: matrix to bidiagonalize (first m×q used)
* @param m Number of rows in the working matrix * @param m Number of rows in the working matrix
* @param q Number of columns in the working matrix (q ≥ p) * @param q Number of columns in the working matrix (q ≥ p)
* @param p Rank = min(m, original_columns) — number of bidiagonalization steps * @param p Rank = min(m, original_columns) — number of bidiagonalization steps
* @param QL Input/output: left Householder accumulation (initialized to identity, * @param QL Input/output: left Householder accumulation (initialized to identity)
* output: QLᵀ such that QLᵀ·W = B) * @param QR Input/output: right Householder accumulation (initialized to identity)
* @param QR Input/output: right Householder accumulation (initialized to identity,
* output: QR such that W·QR = B after left apply)
*/ */
static void Bidiagonalize(Matrix<5, 5> &W, template <uint8_t N>
uint8_t m, uint8_t q, uint8_t p, static void Bidiagonalize(Matrix<N, N> &W, uint8_t m, uint8_t q, uint8_t p,
Matrix<5, 5> &QL, Matrix<N, N> &QL, Matrix<N, N> &QR);
Matrix<5, 5> &QR);
/** /**
* @brief Deflate a bidiagonal matrix by zeroing negligible superdiagonals. * @brief Deflate a bidiagonal matrix by zeroing negligible superdiagonals.
* *
* Scans the p×p upper-bidiagonal matrix stored in W and zeros out any * Scans the p×p upper-bidiagonal matrix stored in W and zeros out any
* superdiagonal element W[i][i+1] whose magnitude is negligible relative to * superdiagonal element W[i][i+1] whose magnitude is negligible relative
* the local diagonal scale (|W[i][i]| + |W[i+1][i+1]|). Deflating splits * to the local diagonal scale (|W[i][i]| + |W[i+1][i+1]|). Deflating
* the matrix into independent unreduced blocks that can each be solved * splits the matrix into independent unreduced blocks that can each be
* separately. * solved separately.
* *
* @param W Input/output: bidiagonal matrix (5×5 working array, first p×p used) * @tparam N Working buffer size
* @param W Input/output: bidiagonal matrix (first p×p used)
* @param p Size of the bidiagonal matrix (min(rows, columns)) * @param p Size of the bidiagonal matrix (min(rows, columns))
* @param tol Relative deflation tolerance (e.g. 1e-8f) * @param tol Relative deflation tolerance (e.g. 1e-8f)
*/ */
static void DeflateBidiagonal(Matrix<5, 5> &W, uint8_t p, float tol); template <uint8_t N>
static void DeflateBidiagonal(Matrix<N, N> &W, uint8_t p, float tol);
/** /**
* @brief Check whether a bidiagonal matrix has fully reduced to diagonal. * @brief Check whether a bidiagonal matrix has fully reduced to diagonal.
@@ -129,12 +172,14 @@ static void DeflateBidiagonal(Matrix<5, 5> &W, uint8_t p, float tol);
* matrix in W is (numerically) zero, i.e. the diagonal entries are the * matrix in W is (numerically) zero, i.e. the diagonal entries are the
* (unsorted) singular values and no unreduced blocks remain. * (unsorted) singular values and no unreduced blocks remain.
* *
* @param W Input: bidiagonal matrix (5×5 working array, first p×p used) * @tparam N Working buffer size
* @param W Input: bidiagonal matrix (first p×p used)
* @param p Size of the bidiagonal matrix (min(rows, columns)) * @param p Size of the bidiagonal matrix (min(rows, columns))
* @param tol Numerical zero threshold multiplier * @param tol Numerical zero threshold multiplier
* @return true when all superdiagonal elements are ~0 * @return true when all superdiagonal elements are ~0
*/ */
static bool BidiagonalIsDiagonal(const Matrix<5, 5> &W, uint8_t p, float tol); template <uint8_t N>
static bool BidiagonalIsDiagonal(const Matrix<N, N> &W, uint8_t p, float tol);
/** /**
* @brief Compute the full SVD of a 2×2 upper-bidiagonal block (pure). * @brief Compute the full SVD of a 2×2 upper-bidiagonal block (pure).
@@ -174,17 +219,22 @@ static void SolveBidiagonalBlock2x2(float a, float b, float d,
* *
* On return: * On return:
* - T's diagonal entries are the eigenvalues (off-diagonals ~0) * - T's diagonal entries are the eigenvalues (off-diagonals ~0)
* - evals[i] = T[i][i], UNSORTED * - evals[i] = T[i][i], UNSORTED, SIGNED (this is a general symmetric
* eigen solver, not just for PSD matrices like T = BᵀB)
* - columns of V are the corresponding eigenvectors (T·V = V·Λ) * - columns of V are the corresponding eigenvectors (T·V = V·Λ)
* *
* @param T Input/output: symmetric matrix (5×5 storage, first n×n used, * Convergence: relative off-diagonal tolerance 1e-10, hard-capped at
* destroyed in place) * 100 sweeps.
* @param n Matrix size (≤ 5) *
* @param evals Output: eigenvalues, unsorted (evals[i] = T[i][i]) * @tparam N Working buffer size (≥ n)
* @param V Output: eigenvector matrix, columns are eigenvectors * @param T Input/output: symmetric matrix (first n×n used, destroyed in place)
* @param n Matrix size
* @param evals Output: eigenvalues, unsorted, length ≥ n
* @param V Output: eigenvector matrix (first n×n used), columns are eigenvectors
*/ */
static void JacobiEigenSymmetric(float T[5][5], uint8_t n, float evals[5], template <uint8_t N>
float V[5][5]); static void JacobiEigenSymmetric(Matrix<N, N> &T, uint8_t n, float *evals,
Matrix<N, N> &V);
/** /**
* @brief Fold a block SVD's factors into the QL/QR accumulators. * @brief Fold a block SVD's factors into the QL/QR accumulators.
@@ -200,22 +250,27 @@ static void JacobiEigenSymmetric(float T[5][5], uint8_t n, float evals[5],
* (e.g. for a wide matrix W = Aᵀ, QL carries n = rows(W) meaningful * (e.g. for a wide matrix W = Aᵀ, QL carries n = rows(W) meaningful
* rows while QR is read back over its first m rows). * rows while QR is read back over its first m rows).
* *
* In-place update is done through temporary buffers (updating QL's block
* columns while still reading them corrupts the result).
*
* @tparam N Working buffer size
* @param blockStart First column/row index of the block in W * @param blockStart First column/row index of the block in W
* @param blockSize Size of the block (2, or > 2 for the Jacobi path) * @param blockSize Size of the block (2, or > 2 for the Jacobi path)
* @param Ublock Left singular-vector factor of the block (blockSize×blockSize in 5×5 storage) * @param Ublock Left singular-vector factor of the block (first blockSize×blockSize used)
* @param Vblock Right singular-vector factor of the block (blockSize×blockSize in 5×5 storage) * @param Vblock Right singular-vector factor of the block (first blockSize×blockSize used)
* @param rowsQL Number of meaningful rows of QL * @param rowsQL Number of meaningful rows of QL
* @param rowsQR Number of meaningful rows of QR * @param rowsQR Number of meaningful rows of QR
* @param QL Input/output: left transformation accumulator * @param QL Input/output: left transformation accumulator
* @param QR Input/output: right transformation accumulator * @param QR Input/output: right transformation accumulator
*/ */
template <uint8_t N>
static void ApplyBlockFactorsToAccumulators(uint8_t blockStart, static void ApplyBlockFactorsToAccumulators(uint8_t blockStart,
uint8_t blockSize, uint8_t blockSize,
const float Ublock[5][5], const Matrix<N, N> &Ublock,
const float Vblock[5][5], const Matrix<N, N> &Vblock,
uint8_t rowsQL, uint8_t rowsQR, uint8_t rowsQL, uint8_t rowsQR,
Matrix<5, 5> &QL, Matrix<N, N> &QL,
Matrix<5, 5> &QR); Matrix<N, N> &QR);
/** /**
* @brief Solve a bidiagonal block larger than 2×2 via Jacobi eigen of BᵀB. * @brief Solve a bidiagonal block larger than 2×2 via Jacobi eigen of BᵀB.
@@ -231,47 +286,49 @@ static void ApplyBlockFactorsToAccumulators(uint8_t blockStart,
* (NOT sqrt(eigenvalue) — forming BᵀB squares the condition number, * (NOT sqrt(eigenvalue) — forming BᵀB squares the condition number,
* causing float noise to swamp true tiny eigenvalues for * causing float noise to swamp true tiny eigenvalues for
* rank-deficient blocks) * rank-deficient blocks)
* 6. Re-sort σ descending, keeping V consistent * 6. Re-sort σ descending, keeping V and B·v consistent
* 7. Build Ublock: uᵢ = B_orig · vᵢ / σᵢ (unit norm); for σᵢ ≈ 0, * 7. Build Ublock: uᵢ = B_orig · vᵢ / σᵢ (unit norm); for σᵢ ≈ 0,
* use Gram-Schmidt orthogonal completion against prior U columns * use Gram-Schmidt orthogonal completion against prior U columns
* 8. Fold Ublock/Vblock into QL/QR via ApplyBlockFactorsToAccumulators * 8. Fold Ublock/Vblock into QL/QR via ApplyBlockFactorsToAccumulators
* 9. Write residual norms into W's diagonal and zero the block's * 9. Write residual norms into W's diagonal and zero the block's
* superdiagonals * superdiagonals
* *
* @param W Input/output: bidiagonal matrix (5×5 working array); the block's * @tparam N Working buffer size (≥ blockSize)
* diagonal holds the singular values and its superdiagonals are * @param W Input/output: bidiagonal matrix; the block's diagonal holds
* zeroed on return * the singular values and its superdiagonals are zeroed on return
* @param blockStart First column/row index of the block * @param blockStart First column/row index of the block
* @param blockSize Size of the block (> 2, ≤ 5) * @param blockSize Size of the block (> 2)
* @param rowsQL Number of meaningful rows of QL * @param rowsQL Number of meaningful rows of QL
* @param rowsQR Number of meaningful rows of QR * @param rowsQR Number of meaningful rows of QR
* @param QL Input/output: left transformation accumulator * @param QL Input/output: left transformation accumulator
* @param QR Input/output: right transformation accumulator * @param QR Input/output: right transformation accumulator
* @param tol (unused: Jacobi convergence tolerance is internal) * @param tol (unused: Jacobi convergence tolerance is internal)
*/ */
static void SolveBidiagonalBlockJacobi(Matrix<5, 5> &W, uint8_t blockStart, template <uint8_t N>
static void SolveBidiagonalBlockJacobi(Matrix<N, N> &W, uint8_t blockStart,
uint8_t blockSize, uint8_t rowsQL, uint8_t blockSize, uint8_t rowsQL,
uint8_t rowsQR, Matrix<5, 5> &QL, uint8_t rowsQR, Matrix<N, N> &QL,
Matrix<5, 5> &QR, float tol); Matrix<N, N> &QR, float tol);
/** /**
* @brief Extract singular values from bidiagonal matrix diagonal and sort. * @brief Extract singular values from bidiagonal matrix diagonal and sort.
* *
* Extracts absolute values of diagonal elements of W as singular values, * Extracts absolute values of diagonal elements of W as singular values,
* then sorts them in descending order while reordering columns of QL * then sorts them in descending order while reordering columns of QL
* and QR to maintain consistency. * and QR to maintain consistency. A negative diagonal element flips the
* sign of the corresponding QL column to keep A = U·Σ·Vᵀ.
* *
* @param W Input: bidiagonal matrix (5×5 working array) * @tparam N Working buffer size
* @param sigma Output: sorted singular values (5×1 column vector, only first p used) * @param W Input: bidiagonal matrix (first p×p used)
* @param sigma Output: sorted singular values (N×1 column vector, only first p used)
* @param p Number of singular values (min(rows, columns)) * @param p Number of singular values (min(rows, columns))
* @param QL Input/output: left transformation matrix (modified during sort) * @param QL Input/output: left transformation matrix (modified during sort)
* @param QR Input/output: right transformation matrix (modified during sort) * @param QR Input/output: right transformation matrix (modified during sort)
*/ */
static void ExtractAndSortSingularValues(Matrix<5, 5> &W, template <uint8_t N>
Matrix<5, 1> &sigma, static void ExtractAndSortSingularValues(Matrix<N, N> &W, Matrix<N, 1> &sigma,
uint8_t p, uint8_t p, Matrix<N, N> &QL,
Matrix<5, 5> &QL, Matrix<N, N> &QR);
Matrix<5, 5> &QR);
/** /**
* @brief Assemble final U and Vt matrices from QL/QR. * @brief Assemble final U and Vt matrices from QL/QR.
@@ -280,23 +337,23 @@ static void ExtractAndSortSingularValues(Matrix<5, 5> &W,
* transposed (Vt) from the accumulated Householder transformations. * transposed (Vt) from the accumulated Householder transformations.
* *
* For non-transpose case: U = QL[:,0:p], Vt = QR[:,0:p]ᵀ * For non-transpose case: U = QL[:,0:p], Vt = QR[:,0:p]ᵀ
* For transpose case: U = QR[:,0:p]ᵀ, Vt = QL[:,0:p]ᵀ * For transpose case: U = QR[:,0:p], Vt = full QLᵀ (all n rows)
* *
* @tparam N Working buffer size (≥ m and ≥ n)
* @param m Number of rows in original matrix * @param m Number of rows in original matrix
* @param n Number of columns in original matrix * @param n Number of columns in original matrix
* @param p Rank = min(m, n) * @param p Rank = min(m, n)
* @param transposeNeeded True if we computed SVD(Aᵀ) instead of SVD(A) * @param transposeNeeded True if we computed SVD(Aᵀ) instead of SVD(A)
* @param QL Left Householder accumulation (5×5) * @param QL Left Householder accumulation (N×N)
* @param QR Right Householder accumulation (5×5) * @param QR Right Householder accumulation (N×N)
* @param U Output: left singular vectors (m×n matrix, only first p columns used) * @param U Output: left singular vectors (N×N, first m×p used)
* @param Vt Output: right singular vectors transposed (n×n matrix, only first p rows used) * @param Vt Output: right singular vectors transposed (N×N, first p×n used)
*/ */
template <uint8_t N>
static void AssembleUAndVt(uint8_t m, uint8_t n, uint8_t p, static void AssembleUAndVt(uint8_t m, uint8_t n, uint8_t p,
bool transposeNeeded, bool transposeNeeded, const Matrix<N, N> &QL,
const Matrix<5, 5> &QL, const Matrix<N, N> &QR, Matrix<N, N> &U,
const Matrix<5, 5> &QR, Matrix<N, N> &Vt);
Matrix<5, 5> &U,
Matrix<5, 5> &Vt);
/** /**
* @brief Compute a Givens rotation that zeros out y. * @brief Compute a Givens rotation that zeros out y.
@@ -318,6 +375,7 @@ static void ComputeGivens(float x, float y, float &c, float &s);
* *
* Applies [c s; -s c] to rows i, j of W (columns startCol..endCol). * Applies [c s; -s c] to rows i, j of W (columns startCol..endCol).
* *
* @tparam N Working buffer size
* @param W Input/output: matrix to transform * @param W Input/output: matrix to transform
* @param i First row index * @param i First row index
* @param j Second row index * @param j Second row index
@@ -326,7 +384,8 @@ static void ComputeGivens(float x, float y, float &c, float &s);
* @param startCol First column to transform * @param startCol First column to transform
* @param endCol Last column to transform * @param endCol Last column to transform
*/ */
static void ApplyGivensLeft(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c, template <uint8_t N>
static void ApplyGivensLeft(Matrix<N, N> &W, uint8_t i, uint8_t j, float c,
float s, uint8_t startCol, uint8_t endCol); float s, uint8_t startCol, uint8_t endCol);
/** /**
@@ -334,6 +393,7 @@ static void ApplyGivensLeft(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c,
* *
* Applies [c -s; s c]ᵀ to columns i, j of W (rows startRow..endRow). * Applies [c -s; s c]ᵀ to columns i, j of W (rows startRow..endRow).
* *
* @tparam N Working buffer size
* @param W Input/output: matrix to transform * @param W Input/output: matrix to transform
* @param i First column index * @param i First column index
* @param j Second column index * @param j Second column index
@@ -342,10 +402,11 @@ static void ApplyGivensLeft(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c,
* @param startRow First row to transform * @param startRow First row to transform
* @param endRow Last row to transform * @param endRow Last row to transform
*/ */
static void ApplyGivensRight(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c, template <uint8_t N>
static void ApplyGivensRight(Matrix<N, N> &W, uint8_t i, uint8_t j, float c,
float s, uint8_t startRow, uint8_t endRow); float s, uint8_t startRow, uint8_t endRow);
} // namespace SVD } // namespace SVD
#ifndef SVD_H_ #ifndef SVD_H_
#include "SVD.cpp" #include "SVD.cpp"
#endif // SVD_H_ #endif
unit-tests/matrix-tests.cpp
+137
View File
@@ -1161,3 +1161,140 @@ TEST_CASE("SVD: 1×2 Row Vector", "Matrix") {
float reconErr = svdReconstructionError(A, U, sigma, Vt); float reconErr = svdReconstructionError(A, U, sigma, Vt);
REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f)); REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
} }
// ============================================================================
// SVD Tests — Large-Size Instantiations (N > 5)
//
// The SVD is templated on N = max(rows, cols) with stack-only buffers, so
// these cases exercise instantiations beyond the old 5×5 hard limit:
// 7×5 (N=7, tall), 6×6 (N=6, square), 5×8 (N=8, wide/transpose path),
// 6×4 (N=6, tall, near rank-deficiency → deflation path).
// Reference singular values: scipy.linalg.svd.
// ============================================================================
TEST_CASE("SVD: Tall 7×5 Matrix (N=7)", "Matrix") {
// Reference: scipy.linalg.svd
// σ = [7.9180769443, 4.6593687008, 4.2921645616, 2.6009010840, 1.9842770351]
Matrix<7, 5> A{-0.7528f, 2.7043f, 1.392f, 0.592f, -2.0639f,
-2.064f, -2.6515f, 2.1971f, 0.6067f, 1.2484f,
-2.8765f, 2.8195f, 1.9947f, -1.726f, -1.9091f,
-1.8996f, -1.1745f, 0.1485f, -0.4083f, -1.2526f,
0.6711f, -2.163f, -1.2471f, -0.8018f, -0.2636f,
1.7111f, -1.802f, 0.0854f, 0.5545f, -2.7213f,
0.6453f, -1.9769f, -2.6097f, 2.6933f, 2.7938f};
Matrix<7, 5> U{};
Matrix<5, 5> Vt{};
Matrix<5, 1> sigma{};
SVD::SVD(A, U, sigma, Vt);
REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(7.9180769443f, 1e-4f));
REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(4.6593687008f, 1e-4f));
REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(4.2921645616f, 1e-4f));
REQUIRE_THAT(sigma.Get(3, 0), Catch::Matchers::WithinRel(2.6009010840f, 1e-4f));
REQUIRE_THAT(sigma.Get(4, 0), Catch::Matchers::WithinRel(1.9842770351f, 1e-4f));
REQUIRE(isSortedDescending(sigma, 5));
REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
float reconErr = svdReconstructionError(A, U, sigma, Vt);
REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
}
TEST_CASE("SVD: Square 6×6 Matrix (N=6)", "Matrix") {
// Reference: scipy.linalg.svd (float32 inputs)
// σ = [5.018912792, 4.244967461, 2.505512476,
// 1.838801861, 0.9111995101, 0.4580149353]
Matrix<6, 6> A{1.2336f, -0.7815f, -1.6093f, 0.7369f, -0.2394f, -1.5118f,
-0.0193f, -1.8624f, 1.6373f, -0.9649f, 0.6501f, -0.7532f,
0.0803f, 0.1868f, -1.2606f, 1.8783f,
1.1005f, 1.758f, 1.5793f, 0.3916f, 1.6875f, -1.646f,
-1.2161f, -1.8191f, -0.6987f, -0.4453f, -0.9146f, 1.315f,
-0.573f, -0.8763f, 0.1708f, -1.4363f, 1.2088f, -1.7018f,
1.089f, 1.9475f};
Matrix<6, 6> U{}, Vt{};
Matrix<6, 1> sigma{};
SVD::SVD(A, U, sigma, Vt);
REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.018912792f, 1e-4f));
REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(4.244967461f, 1e-4f));
REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(2.505512476f, 1e-4f));
REQUIRE_THAT(sigma.Get(3, 0), Catch::Matchers::WithinRel(1.838801861f, 1e-4f));
REQUIRE_THAT(sigma.Get(4, 0), Catch::Matchers::WithinRel(0.9111995101f, 1e-4f));
REQUIRE_THAT(sigma.Get(5, 0), Catch::Matchers::WithinRel(0.4580149353f, 1e-4f));
REQUIRE(isSortedDescending(sigma, 6));
REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
float reconErr = svdReconstructionError(A, U, sigma, Vt);
REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
}
TEST_CASE("SVD: Wide 5×8 Matrix (N=8, transpose path)", "Matrix") {
// Reference: scipy.linalg.svd
// σ = [5.8027782929, 4.1105282764, 3.7755966048, 3.3208483982, 2.0321410547]
//
// Wide matrices take the Aᵀ transpose path; Vt must be the FULL 8×8
// orthogonal matrix (all 8 rows meaningful), not just the top 5.
Matrix<5, 8> A{-1.5064f, -2.4724f, 1.5773f, 1.0343f, 1.145f, 1.3564f, -2.1298f, -0.7077f,
-1.9207f, 1.8155f, 0.6165f, -0.8455f, -2.1822f, -0.9451f, -0.8741f, 1.148f,
0.6878f, 1.9361f, -0.1389f, -1.902f, 1.0662f, 1.3039f, 0.3064f, 1.3548f,
-0.031f, 0.1137f, -0.3623f, -2.3729f, -1.9605f, -2.3429f, 0.6821f, -0.9282f,
0.0429f, 2.0378f, -1.2535f, -0.4481f, 1.2778f, -1.356f, -2.1151f, -1.0512f};
Matrix<5, 8> U{};
Matrix<8, 8> Vt{};
Matrix<8, 1> sigma{};
SVD::SVD(A, U, sigma, Vt);
REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.8027782929f, 1e-4f));
REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(4.1105282764f, 1e-4f));
REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(3.7755966048f, 1e-4f));
REQUIRE_THAT(sigma.Get(3, 0), Catch::Matchers::WithinRel(3.3208483982f, 1e-4f));
REQUIRE_THAT(sigma.Get(4, 0), Catch::Matchers::WithinRel(2.0321410547f, 1e-4f));
// Remaining singular values must be at noise level
REQUIRE(sigma.Get(5, 0) < 1e-3f);
REQUIRE(sigma.Get(6, 0) < 1e-3f);
REQUIRE(sigma.Get(7, 0) < 1e-3f);
REQUIRE(isSortedDescending(sigma, 8));
REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
float reconErr = svdReconstructionError(A, U, sigma, Vt);
REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
}
TEST_CASE("SVD: Tall 6×4 Near Rank-Deficient (N=6, deflation path)", "Matrix") {
// Reference: scipy.linalg.svd
// σ = [5.9434060901, 3.2857910666, 0.3066158795, 6.48e-07]
//
// σ₄ ≈ 6.5e-7 forces the deflation logic to zero the last
// superdiagonal and isolate the trailing 1×1 block.
Matrix<6, 4> A{-0.086904f, 1.410225f, 1.308323f, 2.234762f,
0.022123f, 0.896751f, 0.324176f, 0.773607f,
-0.473015f, 1.555111f, 0.290059f, 1.157726f,
-0.78371f, 1.398884f, -1.930606f, -1.548717f,
0.201518f, -0.626835f, 0.976596f, 0.875294f,
-1.24206f, 1.60595f, -3.078089f, -2.73695f};
Matrix<6, 4> U{};
Matrix<4, 4> Vt{};
Matrix<4, 1> sigma{};
SVD::SVD(A, U, sigma, Vt);
REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.9434060901f, 1e-4f));
REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(3.2857910666f, 1e-4f));
REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(0.3066158795f, 1e-4f));
// Fourth singular value is at noise level (matrix is ~rank 3)
REQUIRE(sigma.Get(3, 0) < 1e-4f);
REQUIRE(isSortedDescending(sigma, 4));
REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
float reconErr = svdReconstructionError(A, U, sigma, Vt);
REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
}
unit-tests/svd-build-blocks-tests.cpp
+11 -1
View File
@@ -1523,8 +1523,18 @@ TEST_CASE("SVD Building Block: JacobiEigenSymmetric", "[Matrix][SVD]") {
S_orig[i][j] = mats[c][i][j]; S_orig[i][j] = mats[c][i][j];
float evals[5] = {0}; float evals[5] = {0};
// JacobiEigenSymmetric operates on Matrix<N,N> — copy the raw test
// data in, run the solver, copy the eigenvector matrix back out.
Matrix<5, 5> Tm{0};
for (int i = 0; i < 5; i++)
for (int j = 0; j < 5; j++)
Tm[i][j] = T[i][j];
Matrix<5, 5> Vm{0};
SVD::JacobiEigenSymmetric(Tm, ns[c], evals, Vm);
float V[5][5] = {{0}}; float V[5][5] = {{0}};
SVD::JacobiEigenSymmetric(T, ns[c], evals, V); for (int i = 0; i < 5; i++)
for (int j = 0; j < 5; j++)
V[i][j] = Vm[i][j];
// 1. Sorted eigenvalues match scipy // 1. Sorted eigenvalues match scipy
float sorted[5] = {0}; float sorted[5] = {0};
unit-tests/svd-reference-values.py
+31
View File
@@ -401,6 +401,37 @@ def main():
("Zero 3x3", np.zeros((3,3))), ("Zero 3x3", np.zeros((3,3))),
("Col vector 2x1", np.array([[3],[4]], dtype=np.float64)), ("Col vector 2x1", np.array([[3],[4]], dtype=np.float64)),
("Row vector 1x2", np.array([[3,4]], dtype=np.float64)), ("Row vector 1x2", np.array([[3,4]], dtype=np.float64)),
# Large-size instantiation cases (N > 5). Literals MUST match the
# C++ test matrices in unit-tests/matrix-tests.cpp exactly, and the
# C++ references use float32 inputs: cast to float32 before svd().
("Tall 7x5", np.array([
[-0.7528, 2.7043, 1.392, 0.592, -2.0639],
[-2.064, -2.6515, 2.1971, 0.6067, 1.2484],
[-2.8765, 2.8195, 1.9947, -1.726, -1.9091],
[-1.8996, -1.1745, 0.1485, -0.4083, -1.2526],
[0.6711, -2.163, -1.2471, -0.8018, -0.2636],
[1.7111, -1.802, 0.0854, 0.5545, -2.7213],
[0.6453, -1.9769, -2.6097, 2.6933, 2.7938]], dtype=np.float32)),
("Square 6x6", np.array([
[1.2336, -0.7815, -1.6093, 0.7369, -0.2394, -1.5118],
[-0.0193, -1.8624, 1.6373, -0.9649, 0.6501, -0.7532],
[0.0803, 0.1868, -1.2606, 1.8783, 1.1005, 1.758],
[1.5793, 0.3916, 1.6875, -1.646, -1.2161, -1.8191],
[-0.6987, -0.4453, -0.9146, 1.315, -0.573, -0.8763],
[0.1708, -1.4363, 1.2088, -1.7018, 1.089, 1.9475]], dtype=np.float32)),
("Wide 5x8", np.array([
[-1.5064, -2.4724, 1.5773, 1.0343, 1.145, 1.3564, -2.1298, -0.7077],
[-1.9207, 1.8155, 0.6165, -0.8455, -2.1822, -0.9451, -0.8741, 1.148],
[0.6878, 1.9361, -0.1389, -1.902, 1.0662, 1.3039, 0.3064, 1.3548],
[-0.031, 0.1137, -0.3623, -2.3729, -1.9605, -2.3429, 0.6821, -0.9282],
[0.0429, 2.0378, -1.2535, -0.4481, 1.2778, -1.356, -2.1151, -1.0512]], dtype=np.float32)),
("Tall 6x4 rank-def", np.array([
[-0.086904, 1.410225, 1.308323, 2.234762],
[0.022123, 0.896751, 0.324176, 0.773607],
[-0.473015, 1.555111, 0.290059, 1.157726],
[-0.78371, 1.398884, -1.930606, -1.548717],
[0.201518, -0.626835, 0.976596, 0.875294],
[-1.24206, 1.60595, -3.078089, -2.73695]], dtype=np.float32)),
] ]
for name, A in test_matrices: for name, A in test_matrices: