diff --git a/src/SVD.cpp b/src/SVD.cpp index 47990ec..53d95f7 100644 --- a/src/SVD.cpp +++ b/src/SVD.cpp @@ -8,10 +8,17 @@ #ifdef SVD_H_ // since the .cpp file has to be included by the .hpp file this // will evaluate to true #include "SVD.hpp" +#include #include // ============================================================================ -// SVD Building Block Implementations +// SVD Building Block Implementations (fully templated, heap-free) +// +// All block operations work on Matrix working buffers with runtime +// bounds. N is the maximum matrix dimension (max(rows, columns) of the +// SVD input). No dynamic allocation is performed anywhere in this file — +// all temporaries are fixed-size arrays whose bounds derive from the +// template parameter N (a compile-time constant per instantiation). // ============================================================================ float SVD::ComputeHouseholder(const float *x, uint8_t len, float *v, @@ -59,7 +66,8 @@ float SVD::ComputeHouseholder(const float *x, uint8_t len, float *v, return norm; } -void SVD::ApplyHouseholderLeft(Matrix<5, 5> &W, const float *v, +template +void SVD::ApplyHouseholderLeft(Matrix &W, const float *v, uint8_t startRow, uint8_t endRow) { uint8_t len = endRow - startRow + 1; @@ -74,8 +82,9 @@ void SVD::ApplyHouseholderLeft(Matrix<5, 5> &W, const float *v, float twoOverVv = 2.0f / vv; - // W = (I - 2vvᵀ) · W - for (uint8_t col = 0; col < 5; col++) { + // W = (I - 2vvᵀ) · W — applied across all N columns; zero-padded + // columns map to zero under the reflection, so this is a no-op there. + for (uint8_t col = 0; col < N; col++) { float dot = 0.0f; for (uint8_t i = 0; i < len; i++) { dot += v[i] * W[startRow + i][col]; @@ -87,7 +96,8 @@ void SVD::ApplyHouseholderLeft(Matrix<5, 5> &W, const float *v, } } -void SVD::ApplyHouseholderRight(Matrix<5, 5> &W, const float *v, +template +void SVD::ApplyHouseholderRight(Matrix &W, const float *v, uint8_t startCol, uint8_t endCol) { uint8_t len = endCol - startCol + 1; @@ -100,8 +110,9 @@ void SVD::ApplyHouseholderRight(Matrix<5, 5> &W, const float *v, float twoOverVv = 2.0f / vv; - // W = W · (I - 2vvᵀ) - for (uint8_t row = 0; row < 5; row++) { + // W = W · (I - 2vvᵀ) — applied across all N rows; zero-padded rows + // map to zero under the reflection, so this is a no-op there. + for (uint8_t row = 0; row < N; row++) { float dot = 0.0f; for (uint8_t i = 0; i < len; i++) { dot += W[row][startCol + i] * v[i]; @@ -126,11 +137,13 @@ void SVD::ApplyHouseholderRight(Matrix<5, 5> &W, const float *v, s = y / r; } -[[gnu::unused]] void SVD::ApplyGivensLeft(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c, +template +[[gnu::unused]] +void SVD::ApplyGivensLeft(Matrix &W, uint8_t i, uint8_t j, float c, float s, uint8_t startCol, uint8_t endCol) { // [c s] [row_i] = [new_row_i] // [-s c] [row_j] [new_row_j] - for (uint8_t col = startCol; col <= endCol && col < 5; col++) { + for (uint8_t col = startCol; col <= endCol && col < N; col++) { float t1 = W[i][col]; float t2 = W[j][col]; W[i][col] = c * t1 + s * t2; @@ -138,11 +151,13 @@ void SVD::ApplyHouseholderRight(Matrix<5, 5> &W, const float *v, } } -[[gnu::unused]] void SVD::ApplyGivensRight(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c, +template +[[gnu::unused]] +void SVD::ApplyGivensRight(Matrix &W, uint8_t i, uint8_t j, float c, float s, uint8_t startRow, uint8_t endRow) { // [col_i col_j] · [c -s] = [new_col_i new_col_j] // [s c] - for (uint8_t row = startRow; row <= endRow && row < 5; row++) { + for (uint8_t row = startRow; row <= endRow && row < N; row++) { float t1 = W[row][i]; float t2 = W[row][j]; W[row][i] = c * t1 + s * t2; @@ -154,15 +169,14 @@ void SVD::ApplyHouseholderRight(Matrix<5, 5> &W, const float *v, // Phase 1: Householder Bidiagonalization // ============================================================================ -void SVD::Bidiagonalize(Matrix<5, 5> &W, - uint8_t m, uint8_t q, uint8_t p, - Matrix<5, 5> &QL, - Matrix<5, 5> &QR) { - // Working matrix W is m×q (padded to 5×5). +template +void SVD::Bidiagonalize(Matrix &W, uint8_t m, uint8_t q, uint8_t p, + Matrix &QL, Matrix &QR) { + // Working matrix W is m×q (padded to N×N). // QL and QR are initialized to identity by the caller. // We reduce W to upper bidiagonal form B using Householder reflections. - float hhVec[5]; // Householder vector storage + float hhVec[N]; // Householder vector storage (N ≥ any len) for (uint8_t k = 0; k < p; k++) { // --- Left Householder on column k, rows k..m-1 --- @@ -173,7 +187,7 @@ void SVD::Bidiagonalize(Matrix<5, 5> &W, continue; // Extract the column segment W[k..k+len-1][k] - float x[5]; + float x[N]; for (uint8_t i = 0; i < len; i++) { x[i] = W[k + i][k]; } @@ -185,7 +199,7 @@ void SVD::Bidiagonalize(Matrix<5, 5> &W, if (alpha == 0.0f) continue; - // Apply H from left to W: W = H·W (columns k..q-1) + // Apply H from left to W: W = H·W (all columns; padding is no-op) SVD::ApplyHouseholderLeft(W, hhVec, k, k + len - 1); // Apply H from right to QL: QL = QL · H @@ -195,14 +209,15 @@ void SVD::Bidiagonalize(Matrix<5, 5> &W, // --- Right Householder on row k, columns k+1..q-1 --- // Zero out elements above the first superdiagonal in row k. // The Householder maps [W[k][k+1], ..., W[k][q-1]] to [gamma, 0, ..., 0], - // preserving the first superdiagonal element (now gamma) and zeroing the rest. + // preserving the first superdiagonal element (now gamma) and zeroing + // the rest. { int len = static_cast(q) - 1 - k; if (len <= 1) continue; // Need at least 2 elements to zero something out // Extract the row segment starting from column k+1 - float x[5]; + float x[N]; for (uint8_t i = 0; i < len; i++) { x[i] = W[k][k + 1 + i]; } @@ -227,7 +242,8 @@ void SVD::Bidiagonalize(Matrix<5, 5> &W, // Phase 2 helpers: block solving of the bidiagonal matrix // ============================================================================ -void SVD::DeflateBidiagonal(Matrix<5, 5> &W, uint8_t p, float tol) { +template +void SVD::DeflateBidiagonal(Matrix &W, uint8_t p, float tol) { // Zero out superdiagonal elements that are negligible relative to the // local diagonal scale. This deflates the bidiagonal matrix into // independent unreduced blocks, each of which can be solved on its own. @@ -243,7 +259,8 @@ void SVD::DeflateBidiagonal(Matrix<5, 5> &W, uint8_t p, float tol) { } } -bool SVD::BidiagonalIsDiagonal(const Matrix<5, 5> &W, uint8_t p, float tol) { +template +bool SVD::BidiagonalIsDiagonal(const Matrix &W, uint8_t p, float tol) { // True when every superdiagonal element of the p×p bidiagonal matrix // has been reduced to (numerically) zero, i.e. the diagonal holds the // singular values and no unreduced blocks remain. @@ -339,8 +356,9 @@ void SVD::SolveBidiagonalBlock2x2(float a, float b, float d, float Ublock[2][2], Ublock[1][1] = u2y; } -void SVD::JacobiEigenSymmetric(float T[5][5], uint8_t n, float evals[5], - float V[5][5]) { +template +void SVD::JacobiEigenSymmetric(Matrix &T, uint8_t n, float *evals, + Matrix &V) { // Cyclic Jacobi eigenvalue algorithm on symmetric n×n matrix T (in place). // On return: // - T is (near-)diagonal; its diagonal entries are the eigenvalues @@ -349,16 +367,16 @@ void SVD::JacobiEigenSymmetric(float T[5][5], uint8_t n, float evals[5], // as V ← V·J so that T·V = V·Λ) float jacTol = 1e-10f; - // V starts as the identity: eigenvector accumulator + // V starts as the identity (first n×n): eigenvector accumulator for (uint8_t i = 0; i < n; i++) for (uint8_t j = 0; j < n; j++) V[i][j] = (i == j) ? 1.0f : 0.0f; for (uint32_t jacIter = 0; jacIter < 100; jacIter++) { - // Check convergence over ALL off-diagonal entries, not just the - // tridiagonal band: cyclic Jacobi on a 3x3+ block fills non-band - // entries (e.g. T[0][2]) during sweeps, so a band-only test can - // declare convergence too early. + // Check convergence over ALL off-diagonal entries of the n×n part, + // not just the tridiagonal band: cyclic Jacobi on a 3x3+ block fills + // non-band entries (e.g. T[0][2]) during sweeps, so a band-only test + // can declare convergence too early. bool converged = true; for (uint8_t i = 0; i < n - 1 && converged; i++) { for (uint8_t j = i + 1; j < n; j++) { @@ -428,12 +446,14 @@ void SVD::JacobiEigenSymmetric(float T[5][5], uint8_t n, float evals[5], } } -void SVD::ApplyBlockFactorsToAccumulators(uint8_t blockStart, uint8_t blockSize, - const float Ublock[5][5], - const float Vblock[5][5], +template +void SVD::ApplyBlockFactorsToAccumulators(uint8_t blockStart, + uint8_t blockSize, + const Matrix &Ublock, + const Matrix &Vblock, uint8_t rowsQL, uint8_t rowsQR, - Matrix<5, 5> &QL, - Matrix<5, 5> &QR) { + Matrix &QL, + Matrix &QR) { // Fold the block SVD factors into the accumulated Householder // transformation matrices: // QL[:, blockStart..blockStart+blockSize-1] ← QL[:, ...] · Ublock @@ -445,16 +465,20 @@ void SVD::ApplyBlockFactorsToAccumulators(uint8_t blockStart, uint8_t blockSize, // for a transposed (wide) problem W = Aᵀ has n rows, so QL carries n // meaningful rows while in the normal case it carries m. // - // BUG FIX: Use temporary buffers to avoid in-place corruption. - // The old code updated QL[j][blockStart+i] while still reading from - // QL[j][blockStart+k] for later i values, corrupting subsequent columns. + // BUG FIX (preserved): Use temporary buffers to avoid in-place + // corruption. The old code updated QL[j][blockStart+i] while still + // reading from QL[j][blockStart+k] for later i values, corrupting + // subsequent columns. + // + // Only the block columns of QL/QR change; temps hold the full N×N + // worst case (compile-time sized, heap-free). - float newQL[5][5] = {{0}}; + float newQL[N][N] = {{0}}; for (uint8_t j = 0; j < rowsQL; j++) { for (uint8_t i = 0; i < blockSize; i++) { float sum = 0.0f; for (uint8_t k = 0; k < blockSize; k++) { - sum += QL[j][blockStart + k] * Ublock[k][i]; + sum += QL[j][blockStart + k] * Ublock.Get(k, i); } newQL[j][blockStart + i] = sum; } @@ -463,12 +487,12 @@ void SVD::ApplyBlockFactorsToAccumulators(uint8_t blockStart, uint8_t blockSize, for (uint8_t i = 0; i < blockSize; i++) QL[j][blockStart + i] = newQL[j][blockStart + i]; - float newQR[5][5] = {{0}}; + float newQR[N][N] = {{0}}; for (uint8_t j = 0; j < rowsQR; j++) { for (uint8_t i = 0; i < blockSize; i++) { float sum = 0.0f; for (uint8_t k = 0; k < blockSize; k++) { - sum += QR[j][blockStart + k] * Vblock[k][i]; + sum += QR[j][blockStart + k] * Vblock.Get(k, i); } newQR[j][blockStart + i] = sum; } @@ -478,31 +502,32 @@ void SVD::ApplyBlockFactorsToAccumulators(uint8_t blockStart, uint8_t blockSize, QR[j][blockStart + i] = newQR[j][blockStart + i]; } -void SVD::SolveBidiagonalBlockJacobi(Matrix<5, 5> &W, uint8_t blockStart, +template +void SVD::SolveBidiagonalBlockJacobi(Matrix &W, uint8_t blockStart, uint8_t blockSize, uint8_t rowsQL, - uint8_t rowsQR, Matrix<5, 5> &QL, - Matrix<5, 5> &QR, float tol) { + uint8_t rowsQR, Matrix &QL, + Matrix &QR, float tol) { // Full SVD of an unreduced upper-bidiagonal block of size > 2, computed // as the eigen-decomposition of the symmetric tridiagonal T = BᵀB: // // 1. Snapshot the ORIGINAL block diagonal/superdiagonal from W // 2. Form T = BᵀB (tridiagonal symmetric) // 3. JacobiEigenSymmetric → eigenvalues + eigenvector matrix V - // 4. Sort eigenvalues descending, reordering V - // 5. Singular values are the RESIDUAL norms σᵢ = ‖B·vᵢ‖ rather than - // sqrt(evals[i]): forming BᵀB squares the condition number, so - // float noise in T swamps the smallest eigenvalues of - // rank-deficient / near-deficient blocks (σ error ~1e-3 instead of - // ~1e-6). ‖B·vᵢ‖ stays accurate to ~eps·‖B‖. It also makes - // uᵢ = B·vᵢ/σᵢ unit-norm by construction; when σᵢ ≈ 0 (true rank - // deficiency) uᵢ is replaced by a Gram–Schmidt orthogonal - // completion so the accumulated QL/QR stay orthogonal. - // 6. Fold Ublock/Vblock into QL/QR via ApplyBlockFactorsToAccumulators - // 7. Only now write σ onto W's diagonal and zero the superdiagonals - (void)tol; // Jacobi convergence tolerance is internal + // 4. Sort eigenvalues descending, reordering V columns + // 5. Compute RESIDUAL singular values: σᵢ = ‖B_orig · vᵢ‖ + // (NOT sqrt(eigenvalue) — forming BᵀB squares the condition + // number, causing float noise to swamp true tiny eigenvalues for + // rank-deficient blocks) + // 6. Re-sort σ descending, keeping V and B·v consistent + // 7. Build Ublock: uᵢ = B_orig · vᵢ / σᵢ (unit norm); for σᵢ ≈ 0, + // use Gram-Schmidt orthogonal completion against prior U columns + // 8. Fold Ublock/Vblock into QL/QR via ApplyBlockFactorsToAccumulators + // 9. Write residual norms into W's diagonal and zero the block's + // superdiagonals - // Step 1: snapshot original block values (diagonal d[i], superdiag e[i]) - float d[5], e[4]; + // Step 1: snapshot the original bidiagonal block from W + float d[N]; // block diagonal + float e[N - 1]; // block superdiagonal for (uint8_t i = 0; i < blockSize; i++) { d[i] = W[blockStart + i][blockStart + i]; } @@ -513,13 +538,14 @@ void SVD::SolveBidiagonalBlockJacobi(Matrix<5, 5> &W, uint8_t blockStart, // Step 2: form T = BᵀB (tridiagonal) // T[i][i] = d[i]² + e[i−1]² (e[−1] = 0) // T[i][i+1] = d[i] · e[i] - float T[5][5] = {{0}}; + // (e[i] lives in column i+1 of B, so it contributes to T[i+1][i+1], + // NOT T[i][i] — do not add e[i]² here.) + Matrix T{0}; for (uint8_t i = 0; i < blockSize; i++) { - float diag = d[i] * d[i]; - if (i > 0) { - diag += e[i - 1] * e[i - 1]; - } - T[i][i] = diag; + float val = d[i] * d[i]; + if (i > 0) + val += e[i - 1] * e[i - 1]; + T[i][i] = val; if (i < blockSize - 1) { float off = d[i] * e[i]; T[i][i + 1] = off; @@ -527,437 +553,450 @@ void SVD::SolveBidiagonalBlockJacobi(Matrix<5, 5> &W, uint8_t blockStart, } } - // Steps 3-4: Jacobi eigen-decomposition, then sort descending - float evals[5] = {0}; - float V[5][5] = {{0}}; + // Step 3: Jacobi eigen decomposition of T + float evals[N]; + Matrix V{0}; SVD::JacobiEigenSymmetric(T, blockSize, evals, V); - for (uint8_t i = 0; i < blockSize - 1; i++) { + // After: evals[i] = T[i][i] (unsorted), V columns are eigenvectors. + // NOTE: T is overwritten in place, so the snapshot d/e from Step 1 is + // used to form B·v later, not the (destroyed) T. + + // Step 4: sort eigenvalues descending, reordering V columns. + // Stable: swap V columns as a whole so vᵢ stays paired with evals[i]. + for (uint8_t i = 0; i < blockSize; i++) { for (uint8_t j = i + 1; j < blockSize; j++) { if (evals[j] > evals[i]) { float tmpE = evals[i]; evals[i] = evals[j]; evals[j] = tmpE; - for (uint8_t k = 0; k < blockSize; k++) { - float tmpV = V[k][i]; - V[k][i] = V[k][j]; - V[k][j] = tmpV; + for (uint8_t r = 0; r < blockSize; r++) { + float tmpV = V[r][i]; + V[r][i] = V[r][j]; + V[r][j] = tmpV; } } } } - // Step 5: residual singular values σᵢ = ‖B·vᵢ‖ and unit-norm uᵢ - float Bv[5][5]; + // Step 5: residual singular values σᵢ = ‖B_orig · vᵢ‖. + // B is upper bidiagonal with diagonal d and superdiagonal e, so: + // (B·v)[k] = d[k]·v[k] + (k < bs-1 ? e[k]·v[k+1] : 0) + // Computing from the ORIGINAL block avoids the double conditioning of + // sqrt(eigenvalue-of-BᵀB), which destroys rank-deficient blocks. + float Bv[N][N]; // Bv[k][i] = (B·vᵢ)[k] for (uint8_t i = 0; i < blockSize; i++) { - for (uint8_t r = 0; r < blockSize; r++) { - float result = d[r] * V[r][i]; - if (r + 1 < blockSize) { - result += e[r] * V[r + 1][i]; - } - Bv[r][i] = result; + for (uint8_t k = 0; k < blockSize; k++) { + float val = d[k] * V[k][i]; + if (k < blockSize - 1) + val += e[k] * V[k + 1][i]; + Bv[k][i] = val; } } - float sigma[5]; + + float sigma[N]; for (uint8_t i = 0; i < blockSize; i++) { - float n = 0.0f; - for (uint8_t r = 0; r < blockSize; r++) { - n += Bv[r][i] * Bv[r][i]; - } - sigma[i] = sqrtf(n); + float s2 = 0.0f; + for (uint8_t k = 0; k < blockSize; k++) + s2 += Bv[k][i] * Bv[k][i]; + sigma[i] = sqrtf(s2); } - // Re-sort sigma descending, swapping V AND Bv columns consistently - // (residual norms can differ slightly in ordering from sqrt(evals)) - for (uint8_t i = 0; i < blockSize - 1; i++) { + + // Step 6: re-sort by σ descending (σ and evals may differ in order due + // to residual computation), moving V and B·v columns together. + for (uint8_t i = 0; i < blockSize; i++) { for (uint8_t j = i + 1; j < blockSize; j++) { if (sigma[j] > sigma[i]) { - float tmpS = sigma[i]; + float ts = sigma[i]; sigma[i] = sigma[j]; - sigma[j] = tmpS; - for (uint8_t k = 0; k < blockSize; k++) { - float tmpV = V[k][i]; - V[k][i] = V[k][j]; - V[k][j] = tmpV; - float tmpB = Bv[k][i]; - Bv[k][i] = Bv[k][j]; - Bv[k][j] = tmpB; + sigma[j] = ts; + for (uint8_t r = 0; r < blockSize; r++) { + float tv = V[r][i]; + V[r][i] = V[r][j]; + V[r][j] = tv; + float tb = Bv[r][i]; + Bv[r][i] = Bv[r][j]; + Bv[r][j] = tb; } } } } - float Ublock[5][5] = {{0}}; - // Pass 1: raw u-columns. For non-degenerate σ, uᵢ = B·vᵢ/σᵢ. + // Step 7: build Ublock (left singular vectors). + // For non-zero σᵢ, uᵢ = B·vᵢ / σᵢ is already unit norm (up to float + // error). For zero σᵢ (rank-deficient block), B·vᵢ ≈ 0 and any unit + // vector orthogonal to the other U columns completes the SVD: + // ||B·Uᵢ|| = ||B·vᵢ|| = σᵢ = 0, and uᵢ·uᵢ = 1. + // We use Gram-Schmidt orthogonal completion (with a basis-vector seed) + // to guarantee Ublock stays orthogonal even when σᵢ is tiny-but- + // non-zero (noise-dominated) — the 2×2 solver uses the simple swap + // trick (only one complement to worry about), but the Jacobi path can + // have many small σᵢ so we re-orthogonalize against ALL prior columns. + Matrix Ublock{0}; for (uint8_t i = 0; i < blockSize; i++) { if (sigma[i] > 1e-30f) { - for (uint8_t r = 0; r < blockSize; r++) { - Ublock[r][i] = Bv[r][i] / sigma[i]; - } - } - // Degenerate (sigma[i] ~ 0): leave zero, completed in pass 2 - } - // Pass 2: enforce a full orthonormal set. Even when B·vᵢ ≠ 0, the - // smallest-σ columns come from eigenvectors of the noise-dominated - // tail of T = BᵀB, so their u-directions are near-random and NOT - // mutually orthogonal. Gram–Schmidt against the previous columns and - // re-normalize (no-op for the well-conditioned columns); if the - // residual is ~0 (true rank deficiency) pick a basis-vector seed and - // complete orthonally instead. - for (uint8_t i = 0; i < blockSize; i++) { - for (uint8_t p = 0; p < i; p++) { - float dot = 0.0f; - for (uint8_t r = 0; r < blockSize; r++) { - dot += Ublock[r][i] * Ublock[r][p]; - } - for (uint8_t r = 0; r < blockSize; r++) { - Ublock[r][i] -= dot * Ublock[r][p]; - } - } - float nn = 0.0f; - for (uint8_t r = 0; r < blockSize; r++) { - nn += Ublock[r][i] * Ublock[r][i]; - } - if (nn > 1e-12f) { - float inv = 1.0f / sqrtf(nn); - for (uint8_t r = 0; r < blockSize; r++) { - Ublock[r][i] *= inv; + for (uint8_t k = 0; k < blockSize; k++) { + Ublock[k][i] = Bv[k][i] / sigma[i]; } } else { - // True rank deficiency: direction arbitrary. Seed with a basis - // vector, Gram–Schmidt against previous columns, normalize. - bool found = false; - for (uint8_t seed = 0; seed < blockSize && !found; seed++) { - float g[5]; - for (uint8_t r = 0; r < blockSize; r++) { - g[r] = (r == seed) ? 1.0f : 0.0f; - } - for (uint8_t p = 0; p < i; p++) { + // Seed with the first basis vector that has meaningful alignment + // with the null-space direction: prefer eᵢ (natural for a + // rank-deficient bidiagonal block), fall back to any eₖ. + float g[N] = {0}; + for (uint8_t k = 0; k < blockSize; k++) { + g[k] = (k == i % blockSize) ? 1.0f : 0.0f; + } + // Re-orthogonalize against all prior U columns (twice for float + // robustness) + for (uint8_t pass = 0; pass < 2; pass++) { + for (uint8_t j = 0; j < i; j++) { float dot = 0.0f; - for (uint8_t r = 0; r < blockSize; r++) { - dot += g[r] * Ublock[r][p]; - } - for (uint8_t r = 0; r < blockSize; r++) { - g[r] -= dot * Ublock[r][p]; - } - } - float gg = 0.0f; - for (uint8_t r = 0; r < blockSize; r++) { - gg += g[r] * g[r]; - } - if (gg > 1e-12f) { - float inv = 1.0f / sqrtf(gg); - for (uint8_t r = 0; r < blockSize; r++) { - Ublock[r][i] = g[r] * inv; - } - found = true; + for (uint8_t k = 0; k < blockSize; k++) + dot += g[k] * Ublock[k][j]; + for (uint8_t k = 0; k < blockSize; k++) + g[k] -= dot * Ublock[k][j]; } } - if (!found) { - // Degenerate fallback (cannot happen for i < blockSize): e_0 - for (uint8_t r = 0; r < blockSize; r++) { - Ublock[r][i] = (r == 0) ? 1.0f : 0.0f; - } + float norm = 0.0f; + for (uint8_t k = 0; k < blockSize; k++) + norm += g[k] * g[k]; + norm = sqrtf(norm); + if (norm < 1e-30f) { + // Degenerate: re-orthogonalization collapsed; force a unit vector + g[i % blockSize] = 1.0f; + norm = 1.0f; } + for (uint8_t k = 0; k < blockSize; k++) + Ublock[k][i] = g[k] / norm; } } - // Step 6: fold the factors into the accumulators + // Second pass: full Gram-Schmidt re-orthogonalization of ALL Ublock + // columns against each other. This kills accumulated error from the + // Jacobi eigensolver AND from the residual-σ normalization (when two + // columns are both noise-dominated they can end up nearly parallel). + // Only touches the i-th column using columns 0..i-1 which are already + // finalized, so in-place is safe here (unlike the QL/QR fold below). + for (uint8_t i = 0; i < blockSize; i++) { + for (uint8_t pass = 0; pass < 2; pass++) { + for (uint8_t j = 0; j < i; j++) { + float dot = 0.0f; + for (uint8_t k = 0; k < blockSize; k++) + dot += Ublock[k][i] * Ublock[k][j]; + for (uint8_t k = 0; k < blockSize; k++) + Ublock[k][i] -= dot * Ublock[k][j]; + } + } + // Re-normalize after subtraction (can shrink slightly) + float norm = 0.0f; + for (uint8_t k = 0; k < blockSize; k++) + norm += Ublock[k][i] * Ublock[k][i]; + norm = sqrtf(norm); + if (norm > 1e-30f) { + for (uint8_t k = 0; k < blockSize; k++) + Ublock[k][i] /= norm; + } else { + // Collapse: pick any basis direction not already used + for (uint8_t k = 0; k < blockSize; k++) + Ublock[k][i] = (k == i) ? 1.0f : 0.0f; + } + } + + // Step 8: fold block factors into the accumulated QL/QR SVD::ApplyBlockFactorsToAccumulators(blockStart, blockSize, Ublock, V, rowsQL, rowsQR, QL, QR); - // Step 7: W last — write singular values onto the diagonal and zero - // the block's superdiagonals + // Step 9: write back the (sorted) residual σᵢ into W's diagonal and + // zero the block's superdiagonal to mark the block as fully reduced. + // ExtractAndSortSingularValues reads the diagonal of W to get σᵢ; + // writing the residual (not sqrt(evals)) keeps the diagonal consistent + // with the uᵢ/vᵢ we just installed, and is the numerically robust + // choice for rank-deficient blocks. for (uint8_t i = 0; i < blockSize; i++) { W[blockStart + i][blockStart + i] = sigma[i]; - if (i < blockSize - 1) { + if (i < blockSize - 1) W[blockStart + i][blockStart + i + 1] = 0.0f; - } } } // ============================================================================ -// Phase 3: Extract and Sort Singular Values +// Phase 3: Extract singular values and assemble U, Σ, Vt // ============================================================================ -void SVD::ExtractAndSortSingularValues(Matrix<5, 5> &W, - Matrix<5, 1> &sigma, - uint8_t p, - Matrix<5, 5> &QL, - Matrix<5, 5> &QR) { +template +void SVD::ExtractAndSortSingularValues(Matrix &W, Matrix &sigma, + uint8_t p, Matrix &QL, + Matrix &QR) { // Extract singular values as absolute values of diagonal elements. - // If a diagonal element is negative, flip the sign of the corresponding - // column in QL to maintain U * Sigma * Vt = A. + // If a diagonal element is negative, flip the sign of the + // corresponding column in QL to maintain U · Σ · Vᵀ = A. This fires + // for diagonal entries that never went through a block solver (the + // solvers always write non-negative σ, so their blocks are no-ops). for (uint8_t i = 0; i < p; i++) { if (W[i][i] < 0.0f) { - // Flip sign of column i in QL - for (uint8_t k = 0; k < 5; k++) { + for (uint8_t k = 0; k < N; k++) { QL[k][i] = -QL[k][i]; } } sigma[i][0] = fabsf(W[i][i]); } - // Sort singular values in descending order and reorder U, V accordingly + // Sort in descending order, reordering U and V columns to match. + // Selection sort: find the max, swap it into position i. for (uint8_t i = 0; i < p - 1; i++) { + uint8_t maxIdx = i; + float maxVal = sigma.Get(i, 0); for (uint8_t j = i + 1; j < p; j++) { - if (sigma[j][0] > sigma[i][0]) { - // Swap singular values - float tmpS = sigma[i][0]; - sigma[i][0] = sigma[j][0]; - sigma[j][0] = tmpS; + float val = sigma.Get(j, 0); + if (val > maxVal) { + maxVal = val; + maxIdx = j; + } + } + if (maxIdx != i) { + // Swap sigma entries + float tmp = sigma.Get(i, 0); + sigma[i][0] = sigma.Get(maxIdx, 0); + sigma[maxIdx][0] = tmp; - // Swap columns of QL - for (uint8_t k = 0; k < 5; k++) { - float tmpQ = QL[k][i]; - QL[k][i] = QL[k][j]; - QL[k][j] = tmpQ; - } + // Swap corresponding columns of QL (left singular vectors) + for (uint8_t k = 0; k < N; k++) { + float tmpQL = QL.Get(k, i); + QL[k][i] = QL.Get(k, maxIdx); + QL[k][maxIdx] = tmpQL; + } - // Swap columns of QR - for (uint8_t k = 0; k < 5; k++) { - float tmpQ = QR[k][i]; - QR[k][i] = QR[k][j]; - QR[k][j] = tmpQ; - } + // Swap corresponding columns of QR (right singular vectors) + for (uint8_t k = 0; k < N; k++) { + float tmpQR = QR.Get(k, i); + QR[k][i] = QR.Get(k, maxIdx); + QR[k][maxIdx] = tmpQR; } } } -} -// ============================================================================ -// Phase 4: Assemble Final U and Vt Matrices -// ============================================================================ + } +template void SVD::AssembleUAndVt(uint8_t m, uint8_t n, uint8_t p, - bool transposeNeeded, - const Matrix<5, 5> &QL, - const Matrix<5, 5> &QR, - Matrix<5, 5> &U, - Matrix<5, 5> &Vt) { - // Initialize output matrices to zero - for (uint8_t i = 0; i < 5; i++) - for (uint8_t j = 0; j < 5; j++) { - U[i][j] = 0; - Vt[i][j] = 0; - } - - // ---- Compute Final U and Vt ---- - // After bidiagonalization, A = QL · W · QRᵀ (QL = product of left - // Householders in application order, QR = product of right Householders), - // and the block solvers fold the block SVD factors in: QL <- QL·U_block, - // QR <- QR·V_block. Hence: + bool transposeNeeded, const Matrix &QL, + const Matrix &QR, Matrix &U, + Matrix &Vt) { + // Assemble the final U (m×k) and Vt (k×n) from the accumulated + // Householder transformations. // - // Non-transposed (m ≥ n): A = (QL)·Σ·(QR)ᵀ - // U = QL (first m rows, first p columns) - // Vt = QRᵀ (first p rows of the n×n matrix) + // For the non-transpose case (m >= n): + // U = first k columns of QL + // Vt = transpose of first k rows of QR (i.e. Vt[i][j] = QR[j][i]) // - // Transposed (wide, m < n): we computed SVD of Aᵀ = QL·W·QRᵀ, so - // A = (QR)·Σ·(QL)ᵀ - // U = QR (first m rows, first p columns) — NOT QRᵀ - // Vt = QLᵀ (the FULL n×n transpose: QL is the left singular-vector - // matrix of Aᵀ and has n = rows(W) meaningful rows, so Vt needs all - // n rows, not just p) + // For the transpose case (m < n, we computed SVD of Aᵀ = V·Σ·Uᵀ): + // U = first k columns of QR (right transforms of Aᵀ = left of A) + // Vt = transpose of QL (full QLᵀ, all n rows) + // + // NOTE: QL and QR were accumulated with the convention + // B = QLᵀ · A · QR (see Bidiagonalize) + // so the left singular vectors of A are columns of QL (not QLᵀ), and + // the right singular vectors of A are columns of QR. - for (uint8_t i = 0; i < m; i++) { - for (uint8_t j = 0; j < n; j++) { - if (j < p) { - if (transposeNeeded) { - U[i][j] = QR.Get(i, j); - } else { - U[i][j] = QL.Get(i, j); - } - } else { - U[i][j] = 0; + if (!transposeNeeded) { + // U = QL[:, 0:p] + for (uint8_t i = 0; i < m; i++) { + for (uint8_t j = 0; j < p; j++) { + U[i][j] = QL.Get(i, j); + } + for (uint8_t j = p; j < n; j++) { + U[i][j] = 0.0f; } } - } - - for (uint8_t i = 0; i < n; i++) { - for (uint8_t j = 0; j < n; j++) { - if (transposeNeeded) { + // Vt = QRᵀ [0:p, :] + for (uint8_t i = 0; i < n; i++) { + for (uint8_t j = 0; j < n; j++) { + if (i < p) { + Vt[i][j] = QR.Get(j, i); + } else { + Vt[i][j] = 0.0f; + } + } + } + } else { + // U = QR[:, 0:p] (U is rows×columns = n×m; only first p columns + // meaningful, remaining columns zero) + for (uint8_t i = 0; i < n; i++) { + for (uint8_t j = 0; j < p; j++) { + U[i][j] = QR.Get(i, j); + } + for (uint8_t j = p; j < m; j++) { + U[i][j] = 0.0f; + } + } + // Vt = QLᵀ — the FULL m×m transpose. QL is the left factor of Aᵀ and + // has m = rows(Aᵀ) = columns(A) = N meaningful rows, so Vt (N×N) + // needs ALL m rows, not just n or p. Using a smaller bound here + // leaves trailing Vt rows zero and breaks both orthogonality and the + // U·Σ·Vᵀ = A reconstruction of the last columns of A. + for (uint8_t i = 0; i < m; i++) { + for (uint8_t j = 0; j < m; j++) { Vt[i][j] = QL.Get(j, i); - } else if (i < p) { - Vt[i][j] = QR.Get(j, i); - } else { - Vt[i][j] = 0; } } } } // ============================================================================ -// SVD Implementation - Golub-Kahan-Reinsch Algorithm +// Main SVD Function // ============================================================================ -/** - * @brief SVD for any m×n matrix using Householder bidiagonalization + - * implicit QR iteration on the bidiagonal form. - * - * Given A (m×n), computes U (m×k), Σ (k×k diagonal), Vᵀ (k×n) where - * k = min(m,n) and A = U·Σ·Vᵀ. - * - * We store results as: - * - U: Matrix — first k columns are meaningful - * - sigma: Matrix — first k entries are non-zero singular values - * - Vt: Matrix — first k rows are meaningful - * - * For m < n (wide matrices), we work with Aᵀ and swap roles of U and V. - */ + template void SVD::SVD(Matrix &matrixToDecompose, Matrix &U, Matrix &sigma, Matrix &Vt) { - static_assert(rows <= 5 && columns <= 5, - "SVD currently supports matrices up to 5×5"); + // N = max(rows, columns): the working-buffer dimension. All internal + // temporaries are Matrix on the stack (heap-free). See the + // header for the stack-usage estimate (≈ 11·N² floats at peak). + constexpr uint8_t N = (rows > columns) ? rows : columns; - uint8_t m = rows; - uint8_t n = columns; - uint8_t p = (m < n) ? m : n; // rank = min(m,n) + // For wide matrices (rows < columns), compute SVD of Aᵀ (which is + // tall), then swap back. This keeps all Householder logic in the + // tall case. + const bool transposeNeeded = rows < columns; + const uint8_t m = transposeNeeded ? columns : rows; + const uint8_t n = transposeNeeded ? rows : columns; + const uint8_t p = m < n ? m : n; // rank (min of m and n) - // For wide matrices (m < n), work with Aᵀ instead. - // SVD(A) = U·Σ·Vᵀ ⟺ SVD(Aᵀ) = V·Σ·Uᵀ - // So if we compute SVD(Aᵀ) = Ũ·Σ·Ṽᵀ, then U = Ṽ and Vt = Ũᵀ. - bool transposeNeeded = (m < n); + // Working matrices (N×N, heap-free stack storage) + Matrix W{0}; + Matrix QL{0}; + Matrix QR{0}; + Matrix UInternal{0}; + Matrix VtInternal{0}; + Matrix sigmaInternal{0}; - // Working matrix: always p×p or larger square - Matrix<5, 5> W{0}; - for (uint8_t i = 0; i < m; i++) { - for (uint8_t j = 0; j < n; j++) { - float val = matrixToDecompose.Get(i, j); - if (transposeNeeded) { - W[j][i] = val; // store Aᵀ - } else { - W[i][j] = val; + // Fill W with the working matrix (A or Aᵀ) + if (transposeNeeded) { + for (uint8_t i = 0; i < m; i++) { + for (uint8_t j = 0; j < n; j++) { + W[i][j] = matrixToDecompose.Get(j, i); + } + } + } else { + for (uint8_t i = 0; i < m; i++) { + for (uint8_t j = 0; j < n; j++) { + W[i][j] = matrixToDecompose.Get(i, j); } } } - // After bidiagonalization, W holds the bidiagonal matrix B. - // Q_L and Q_R accumulate the Householder transformations. - Matrix<5, 5> QL{0}, QR{0}; - for (uint8_t i = 0; i < 5; i++) { - QL[i][i] = 1; - QR[i][i] = 1; + // Initialize QL and QR to identity (N×N). NOTE: Matrix::Identity() is + // a static factory returning by value — a bare call would be a no-op. + for (uint8_t i = 0; i < N; i++) { + QL[i][i] = 1.0f; + QR[i][i] = 1.0f; } - // ---- Phase 1: Householder Bidiagonalization ---- - // For non-transpose (m ≥ n): W is m×n, bidiagonalize to get B (m×n) - // For transpose (m < n): W is n×m (= Aᵀ), bidiagonalize to get B (n×m) - // Pass the ACTUAL dimensions of W: Bidiagonalize needs the full row count - // so that the last left Householder (k = p-1, len = rowsW - k) folds the - // extra rows into the last diagonal element and zeros them out. - { - uint8_t rowsW = transposeNeeded ? n : m; // rows of W - uint8_t colsW = transposeNeeded ? m : n; // columns of W - SVD::Bidiagonalize(W, rowsW, colsW, p, QL, QR); - } + // Phase 1: Householder bidiagonalization + SVD::Bidiagonalize(W, m, n, p, QL, QR); - // ---- Phase 2: QR Iteration on Bidiagonal Matrix --> - // W now contains the upper bidiagonal matrix B. - // We apply QR iterations to converge superdiagonal elements to zero, - // leaving singular values on the diagonal. + // Phase 2: reduce the bidiagonal matrix to diagonal form by solving + // each unreduced block independently. // - // rowsQL / rowsQR are the meaningful row extents of the accumulators. - // QL is the LEFT factor of the bidiagonalized working matrix W, so it - // carries rowsW = (transposeNeeded ? n : m) meaningful rows: in the wide - // (transposed) case Vt = QLᵀ needs ALL n rows, so block factors must be - // applied over 0..n−1. QR is only ever read back over its first m rows - // (as U), but applying factors over all n rows is harmless and matches - // the full-row Householder application in Bidiagonalize. - uint8_t rowsQL = transposeNeeded ? n : m; - uint8_t rowsQR = n; - - uint32_t maxIter = 1000; - float tol = 1e-8f; - - for (uint32_t iter = 0; iter < maxIter; iter++) { - // Deflate: zero out negligible SUPERDIAGONAL elements - SVD::DeflateBidiagonal(W, p, tol); - - // If all superdiagonal elements are zero, we're done - if (SVD::BidiagonalIsDiagonal(W, p, tol)) - break; - - // Process all unreduced blocks in the matrix - bool processedAny = false; + // Deflation zeros out negligible superdiagonals, splitting the + // bidiagonal matrix into independent blocks. Each block of size 2 is + // solved in closed form; blocks larger than 2 use a cyclic Jacobi + // eigen-solve of BᵀB. + SVD::DeflateBidiagonal(W, p, 1e-8f); + if (!SVD::BidiagonalIsDiagonal(W, p, 1e-10f)) { uint8_t blockStart = 0; - - while (blockStart < p - 1) { - // Find end of current unreduced block + bool processedAny = false; + while (blockStart < p) { + // Find the end of the current unreduced block uint8_t blockEnd = blockStart; while (blockEnd < p - 1 && - fabsf(W[blockEnd][blockEnd + 1]) > - tol * fmaxf(fabsf(W[blockEnd][blockEnd]) + - fabsf(W[blockEnd + 1][blockEnd + 1]), - 1e-10f)) { + fabsf(W.Get(blockEnd, blockEnd + 1)) > 0.0f) { blockEnd++; } - - // blockStart..blockEnd is an unreduced block of size - // (blockEnd - blockStart + 1) uint8_t blockSize = blockEnd - blockStart + 1; if (blockSize == 2) { - // Handle 2×2 block directly using closed-form solution - float Ub2[2][2] = {{0}}, Vb2[2][2] = {{0}}; - float sig[2] = {0, 0}; - SVD::SolveBidiagonalBlock2x2(W[blockStart][blockStart], - W[blockStart][blockEnd], - W[blockEnd][blockEnd], Ub2, Vb2, sig); + // 2×2 block: closed-form SVD + float a = W.Get(blockStart, blockStart); + float b = W.Get(blockStart, blockStart + 1); + float d = W.Get(blockStart + 1, blockStart + 1); - float Ublock[5][5] = {{0}}, Vblock[5][5] = {{0}}; - for (uint8_t i = 0; i < 2; i++) + float Ublock2[2][2], Vblock2[2][2], sigma2[2]; + SVD::SolveBidiagonalBlock2x2(a, b, d, Ublock2, Vblock2, sigma2); + + // Expand the 2×2 block factors into N×N working buffers, then + // fold them into QL and QR (the snapshot-before-compute pattern + // in ApplyBlockFactorsToAccumulators avoids the in-place + // corruption of the original code). + Matrix Ublock{0}; + Matrix Vblock{0}; + for (uint8_t i = 0; i < 2; i++) { for (uint8_t j = 0; j < 2; j++) { - Ublock[i][j] = Ub2[i][j]; - Vblock[i][j] = Vb2[i][j]; + Ublock[i][j] = Ublock2[i][j]; + Vblock[i][j] = Vblock2[i][j]; } + } - // Apply block factors to the accumulators SVD::ApplyBlockFactorsToAccumulators(blockStart, 2, Ublock, Vblock, - rowsQL, rowsQR, QL, QR); + m, n, QL, QR); - // Store singular values on diagonal, zero the superdiagonal - W[blockStart][blockStart] = sig[0]; - W[blockEnd][blockEnd] = sig[1]; - W[blockStart][blockEnd] = 0; + // Write the sorted singular values back into W's diagonal and + // zero the block's superdiagonal. + W[blockStart][blockStart] = sigma2[0]; + W[blockStart + 1][blockStart + 1] = sigma2[1]; + W[blockStart][blockStart + 1] = 0.0f; } else if (blockSize > 2) { - // Larger blocks: SVD via eigendecomposition of BᵀB (Jacobi) - SVD::SolveBidiagonalBlockJacobi(W, blockStart, blockSize, rowsQL, - rowsQR, QL, QR, tol); + // Larger block: cyclic Jacobi eigen-solve of BᵀB + SVD::SolveBidiagonalBlockJacobi(W, blockStart, blockSize, m, n, QL, + QR, 1e-10f); } + // Move to the next block + blockStart = blockEnd + 1; processedAny = true; - blockStart = blockEnd + 1; // Move to next block - } - - if (!processedAny) { - // No unreduced blocks found, but superdiagonal is not all zero - // This can happen with numerical issues, just break - break; } + (void)processedAny; } - // ---- Phase 3: Extract and Sort Singular Values ---- - // Use internal 5×1 buffer for sigma - Matrix<5, 1> sigmaInternal{0}; - ExtractAndSortSingularValues(W, sigmaInternal, p, QL, QR); + // Phase 3: extract and sort singular values, assemble U and Vt + SVD::ExtractAndSortSingularValues(W, sigmaInternal, p, QL, QR); + SVD::AssembleUAndVt(m, n, p, transposeNeeded, QL, QR, UInternal, + VtInternal); - // ---- Phase 4: Assemble Final U and Vt ---- - // Use internal 5×5 buffers for U and Vt - Matrix<5, 5> UInternal{0}, VtInternal{0}; - AssembleUAndVt(m, n, p, transposeNeeded, QL, QR, UInternal, VtInternal); - - // Copy results to output parameters - for (uint8_t i = 0; i < columns; i++) { - sigma[i][0] = sigmaInternal.Get(i, 0); - } - for (uint8_t i = 0; i < rows; i++) { - for (uint8_t j = 0; j < columns; j++) { - U[i][j] = UInternal.Get(i, j); + // Copy results into the output matrices + if (transposeNeeded) { + // U (m×n = rows×columns) from UInternal (N×N) + for (uint8_t i = 0; i < rows; i++) { + for (uint8_t j = 0; j < columns; j++) { + U[i][j] = UInternal.Get(i, j); + } } - } - // Vt is a columns×columns (n×n) matrix: BOTH bounds must run over columns. - for (uint8_t i = 0; i < columns; i++) { - for (uint8_t j = 0; j < columns; j++) { - Vt[i][j] = VtInternal.Get(i, j); + // sigma (n×1 = columns×1) from sigmaInternal + for (uint8_t i = 0; i < columns; i++) { + sigma[i][0] = sigmaInternal.Get(i, 0); + } + // Vt (n×n = columns×columns) from VtInternal + for (uint8_t i = 0; i < columns; i++) { + for (uint8_t j = 0; j < columns; j++) { + Vt[i][j] = VtInternal.Get(i, j); + } + } + } else { + // U (m×n = rows×columns) from UInternal + for (uint8_t i = 0; i < rows; i++) { + for (uint8_t j = 0; j < columns; j++) { + U[i][j] = UInternal.Get(i, j); + } + } + // sigma (n×1 = columns×1) from sigmaInternal + for (uint8_t i = 0; i < columns; i++) { + sigma[i][0] = sigmaInternal.Get(i, 0); + } + // Vt (n×n = columns×columns) from VtInternal + for (uint8_t i = 0; i < columns; i++) { + for (uint8_t j = 0; j < columns; j++) { + Vt[i][j] = VtInternal.Get(i, j); + } } } } diff --git a/src/SVD.hpp b/src/SVD.hpp index 7c88362..dfa579f 100644 --- a/src/SVD.hpp +++ b/src/SVD.hpp @@ -3,6 +3,21 @@ /** * @brief library that uses Matrix.hpp and performs SVD on a matrix + * + * @note Fully templated: SVD works for ANY Matrix 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 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 { /** @@ -18,14 +33,28 @@ namespace SVD { * The decomposition satisfies: A ≈ U × diag(Σ) × Vᵀ * Singular values are returned in descending order. * - * @param U Output: left singular vectors (m×k orthogonal matrix) - * @param sigma Output: singular values as k×1 column vector, sorted descending - * @param Vt Output: right singular vectors transposed (k×n matrix) + * Output storage conventions: + * - U: Matrix — first k columns are meaningful + * (rows k..columns−1 are zero in the wide case) + * - sigma: Matrix — first k entries are the singular + * values; entries beyond k (wide matrices only) are zero + * - Vt: Matrix — first k rows are meaningful + * (zero-padded in the tall case) * - * @note This implementation uses the Golub-Kahan-Reinsch algorithm: - * 1. Householder bidiagonalization of A - * 2. Implicit QR iteration on the bidiagonal matrix - * 3. Accumulation of U and V factors throughout + * For wide matrices (rows < columns) the SVD is computed on Aᵀ and the + * factors are swapped back. + * + * @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 void SVD(Matrix &matrixToDecompose, Matrix &U, @@ -33,7 +62,14 @@ void SVD(Matrix &matrixToDecompose, Matrix &U, // ======================================================================== // 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 &matrixToDecompose, Matrix &U, * Given input vector x, computes normalized v and scalar alpha such that: * (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 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 * @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. * - * 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 startRow First row index * @param endRow Last row index */ -static void ApplyHouseholderLeft(Matrix<5, 5> &W, const float *v, +template +static void ApplyHouseholderLeft(Matrix &W, const float *v, uint8_t startRow, uint8_t endRow); /** * @brief Apply a Householder reflection from the right. * * 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 startCol First column index * @param endCol Last column index */ -static void ApplyHouseholderRight(Matrix<5, 5> &W, const float *v, +template +static void ApplyHouseholderRight(Matrix &W, const float *v, uint8_t startCol, uint8_t endCol); /** * @brief Reduce a matrix to upper bidiagonal form using Householder reflections. * * 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 - * the left and right transformation matrices in QL and QR respectively. + * W (m×q, where q ≥ p, stored in N×N working space) to upper bidiagonal + * form B (p×q), accumulating the left and right transformation matrices + * in QL and QR respectively. * * Algorithm (Golub-Kahan bidiagonalization): * For k = 0 to p-1: @@ -93,34 +136,34 @@ static void ApplyHouseholderRight(Matrix<5, 5> &W, const float *v, * The accumulated transformations satisfy: * 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 q Number of columns in the working matrix (q ≥ p) * @param p Rank = min(m, original_columns) — number of bidiagonalization steps - * @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, - * output: QR such that W·QR = B after left apply) + * @param QL Input/output: left Householder accumulation (initialized to identity) + * @param QR Input/output: right Householder accumulation (initialized to identity) */ -static void Bidiagonalize(Matrix<5, 5> &W, - uint8_t m, uint8_t q, uint8_t p, - Matrix<5, 5> &QL, - Matrix<5, 5> &QR); +template +static void Bidiagonalize(Matrix &W, uint8_t m, uint8_t q, uint8_t p, + Matrix &QL, Matrix &QR); /** * @brief Deflate a bidiagonal matrix by zeroing negligible superdiagonals. * * 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 - * the local diagonal scale (|W[i][i]| + |W[i+1][i+1]|). Deflating splits - * the matrix into independent unreduced blocks that can each be solved - * separately. + * superdiagonal element W[i][i+1] whose magnitude is negligible relative + * to the local diagonal scale (|W[i][i]| + |W[i+1][i+1]|). Deflating + * splits the matrix into independent unreduced blocks that can each be + * 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 tol Relative deflation tolerance (e.g. 1e-8f) */ -static void DeflateBidiagonal(Matrix<5, 5> &W, uint8_t p, float tol); +template +static void DeflateBidiagonal(Matrix &W, uint8_t p, float tol); /** * @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 * (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 tol Numerical zero threshold multiplier * @return true when all superdiagonal elements are ~0 */ -static bool BidiagonalIsDiagonal(const Matrix<5, 5> &W, uint8_t p, float tol); +template +static bool BidiagonalIsDiagonal(const Matrix &W, uint8_t p, float tol); /** * @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: * - 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·Λ) * - * @param T Input/output: symmetric matrix (5×5 storage, first n×n used, - * destroyed in place) - * @param n Matrix size (≤ 5) - * @param evals Output: eigenvalues, unsorted (evals[i] = T[i][i]) - * @param V Output: eigenvector matrix, columns are eigenvectors + * Convergence: relative off-diagonal tolerance 1e-10, hard-capped at + * 100 sweeps. + * + * @tparam N Working buffer size (≥ n) + * @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], - float V[5][5]); +template +static void JacobiEigenSymmetric(Matrix &T, uint8_t n, float *evals, + Matrix &V); /** * @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 * 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 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 Vblock Right 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 (first blockSize×blockSize used) * @param rowsQL Number of meaningful rows of QL * @param rowsQR Number of meaningful rows of QR * @param QL Input/output: left transformation accumulator * @param QR Input/output: right transformation accumulator */ +template static void ApplyBlockFactorsToAccumulators(uint8_t blockStart, uint8_t blockSize, - const float Ublock[5][5], - const float Vblock[5][5], + const Matrix &Ublock, + const Matrix &Vblock, uint8_t rowsQL, uint8_t rowsQR, - Matrix<5, 5> &QL, - Matrix<5, 5> &QR); + Matrix &QL, + Matrix &QR); /** * @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, * causing float noise to swamp true tiny eigenvalues for * 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, * use Gram-Schmidt orthogonal completion against prior U columns * 8. Fold Ublock/Vblock into QL/QR via ApplyBlockFactorsToAccumulators * 9. Write residual norms into W's diagonal and zero the block's * superdiagonals * - * @param W Input/output: bidiagonal matrix (5×5 working array); the block's - * diagonal holds the singular values and its superdiagonals are - * zeroed on return + * @tparam N Working buffer size (≥ blockSize) + * @param W Input/output: bidiagonal matrix; the block's diagonal holds + * the singular values and its superdiagonals are zeroed on return * @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 rowsQR Number of meaningful rows of QR * @param QL Input/output: left transformation accumulator * @param QR Input/output: right transformation accumulator * @param tol (unused: Jacobi convergence tolerance is internal) */ -static void SolveBidiagonalBlockJacobi(Matrix<5, 5> &W, uint8_t blockStart, +template +static void SolveBidiagonalBlockJacobi(Matrix &W, uint8_t blockStart, uint8_t blockSize, uint8_t rowsQL, - uint8_t rowsQR, Matrix<5, 5> &QL, - Matrix<5, 5> &QR, float tol); + uint8_t rowsQR, Matrix &QL, + Matrix &QR, float tol); /** * @brief Extract singular values from bidiagonal matrix diagonal and sort. * * Extracts absolute values of diagonal elements of W as singular values, * 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) - * @param sigma Output: sorted singular values (5×1 column vector, only first p used) + * @tparam N Working buffer size + * @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 QL Input/output: left transformation matrix (modified during sort) * @param QR Input/output: right transformation matrix (modified during sort) */ -static void ExtractAndSortSingularValues(Matrix<5, 5> &W, - Matrix<5, 1> &sigma, - uint8_t p, - Matrix<5, 5> &QL, - Matrix<5, 5> &QR); +template +static void ExtractAndSortSingularValues(Matrix &W, Matrix &sigma, + uint8_t p, Matrix &QL, + Matrix &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. * * 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 n Number of columns in original matrix * @param p Rank = min(m, n) * @param transposeNeeded True if we computed SVD(Aᵀ) instead of SVD(A) - * @param QL Left Householder accumulation (5×5) - * @param QR Right Householder accumulation (5×5) - * @param U Output: left singular vectors (m×n matrix, only first p columns used) - * @param Vt Output: right singular vectors transposed (n×n matrix, only first p rows used) + * @param QL Left Householder accumulation (N×N) + * @param QR Right Householder accumulation (N×N) + * @param U Output: left singular vectors (N×N, first m×p used) + * @param Vt Output: right singular vectors transposed (N×N, first p×n used) */ +template static void AssembleUAndVt(uint8_t m, uint8_t n, uint8_t p, - bool transposeNeeded, - const Matrix<5, 5> &QL, - const Matrix<5, 5> &QR, - Matrix<5, 5> &U, - Matrix<5, 5> &Vt); + bool transposeNeeded, const Matrix &QL, + const Matrix &QR, Matrix &U, + Matrix &Vt); /** * @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). * + * @tparam N Working buffer size * @param W Input/output: matrix to transform * @param i First 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 endCol Last column to transform */ -static void ApplyGivensLeft(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c, +template +static void ApplyGivensLeft(Matrix &W, uint8_t i, uint8_t j, float c, 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). * + * @tparam N Working buffer size * @param W Input/output: matrix to transform * @param i First 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 endRow Last row to transform */ -static void ApplyGivensRight(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c, +template +static void ApplyGivensRight(Matrix &W, uint8_t i, uint8_t j, float c, float s, uint8_t startRow, uint8_t endRow); } // namespace SVD #ifndef SVD_H_ #include "SVD.cpp" -#endif // SVD_H_ \ No newline at end of file +#endif diff --git a/unit-tests/matrix-tests.cpp b/unit-tests/matrix-tests.cpp index f32035a..32e6c56 100644 --- a/unit-tests/matrix-tests.cpp +++ b/unit-tests/matrix-tests.cpp @@ -1160,4 +1160,141 @@ TEST_CASE("SVD: 1×2 Row Vector", "Matrix") { float reconErr = svdReconstructionError(A, U, sigma, Vt); REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f)); -} \ No newline at end of file +} +// ============================================================================ +// 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)); +} diff --git a/unit-tests/svd-build-blocks-tests.cpp b/unit-tests/svd-build-blocks-tests.cpp index 37932d5..fd03ef7 100644 --- a/unit-tests/svd-build-blocks-tests.cpp +++ b/unit-tests/svd-build-blocks-tests.cpp @@ -1523,8 +1523,18 @@ TEST_CASE("SVD Building Block: JacobiEigenSymmetric", "[Matrix][SVD]") { S_orig[i][j] = mats[c][i][j]; float evals[5] = {0}; + // JacobiEigenSymmetric operates on Matrix — 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}}; - 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 float sorted[5] = {0}; diff --git a/unit-tests/svd-reference-values.py b/unit-tests/svd-reference-values.py index c274210..d9e8504 100644 --- a/unit-tests/svd-reference-values.py +++ b/unit-tests/svd-reference-values.py @@ -401,6 +401,37 @@ def main(): ("Zero 3x3", np.zeros((3,3))), ("Col vector 2x1", 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: