diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a1dcb05..531b8aa 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -56,4 +56,21 @@ target_link_libraries(matrix set_target_properties(matrix PROPERTIES LINKER_LANGUAGE CXX +) + +# SVD +add_library(svd + STATIC + SVD.cpp +) + +target_link_libraries(svd + PUBLIC + vector-3d-intf + PRIVATE +) + +set_target_properties(svd + PROPERTIES + LINKER_LANGUAGE CXX ) \ No newline at end of file diff --git a/src/SVD.cpp b/src/SVD.cpp new file mode 100644 index 0000000..f63eb92 --- /dev/null +++ b/src/SVD.cpp @@ -0,0 +1,540 @@ +// This #ifndef section makes clangd happy so that it can properly do type hints +// in this file +#ifndef SVD_H_ +#define SVD_H_ +#include "SVD.hpp" +#endif + +#ifdef SVD_H_ // since the .cpp file has to be included by the .hpp file this + // will evaluate to true +#include "SVD.hpp" +#include + +// ============================================================================ +// SVD Building Block Implementations +// ============================================================================ + +float SVD::ComputeHouseholder(const float *x, uint8_t len, float *v, + float &alpha) { + // Compute ||x|| + float norm = 0.0f; + for (uint8_t i = 0; i < len; i++) { + norm += x[i] * x[i]; + } + norm = sqrtf(norm); + + if (norm < 1e-30f) { + alpha = 0.0f; + for (uint8_t i = 0; i < len; i++) { + v[i] = 0.0f; + } + return 0.0f; + } + + // Choose sign to avoid cancellation: alpha has opposite sign of x[0] + alpha = (x[0] >= 0.0f) ? -norm : norm; + + // v = x - alpha * e1, then normalize + float v0 = x[0] - alpha; + + // Compute ||v||² directly: v0² + x₁² + ... + xₙ₋₁² + float vv = v0 * v0; + for (uint8_t i = 1; i < len; i++) { + vv += x[i] * x[i]; + } + + if (vv < 1e-30f) { + // Already aligned with e1 + for (uint8_t i = 0; i < len; i++) { + v[i] = (i == 0) ? 1.0f : 0.0f; + } + return norm; + } + + float scale = 1.0f / sqrtf(vv); + for (uint8_t i = 0; i < len; i++) { + v[i] = (i == 0) ? v0 * scale : x[i] * scale; + } + + return norm; +} + +void SVD::ApplyHouseholderLeft(Matrix<5, 5> &W, const float *v, + uint8_t startRow, uint8_t endRow) { + uint8_t len = endRow - startRow + 1; + + // Compute vᵀv (should be 2.0 for our normalized vectors, but compute + // explicitly) + float vv = 0.0f; + for (uint8_t i = 0; i < len; i++) { + vv += v[i] * v[i]; + } + if (vv < 1e-30f) + return; + + float twoOverVv = 2.0f / vv; + + // W = (I - 2vvᵀ) · W + for (uint8_t col = 0; col < 5; col++) { + float dot = 0.0f; + for (uint8_t i = 0; i < len; i++) { + dot += v[i] * W[startRow + i][col]; + } + dot *= twoOverVv; + for (uint8_t i = 0; i < len; i++) { + W[startRow + i][col] -= dot * v[i]; + } + } +} + +void SVD::ApplyHouseholderRight(Matrix<5, 5> &W, const float *v, + uint8_t startCol, uint8_t endCol) { + uint8_t len = endCol - startCol + 1; + + float vv = 0.0f; + for (uint8_t i = 0; i < len; i++) { + vv += v[i] * v[i]; + } + if (vv < 1e-30f) + return; + + float twoOverVv = 2.0f / vv; + + // W = W · (I - 2vvᵀ) + for (uint8_t row = 0; row < 5; row++) { + float dot = 0.0f; + for (uint8_t i = 0; i < len; i++) { + dot += W[row][startCol + i] * v[i]; + } + dot *= twoOverVv; + for (uint8_t i = 0; i < len; i++) { + W[row][startCol + i] -= dot * v[i]; + } + } +} + +void SVD::ComputeGivens(float x, float y, float &c, float &s) { + float r = sqrtf(x * x + y * y); + + if (r < 1e-30f) { + c = 1.0f; + s = 0.0f; + return; + } + + c = x / r; + s = y / r; +} + +void SVD::ApplyGivensLeft(Matrix<5, 5> &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++) { + float t1 = W[i][col]; + float t2 = W[j][col]; + W[i][col] = c * t1 + s * t2; + W[j][col] = -s * t1 + c * t2; + } +} + +void SVD::ApplyGivensRight(Matrix<5, 5> &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++) { + float t1 = W[row][i]; + float t2 = W[row][j]; + W[row][i] = c * t1 + s * t2; + W[row][j] = -s * t1 + c * t2; + } +} + +// ============================================================================ +// SVD Implementation - Golub-Kahan-Reinsch Algorithm +// ============================================================================ +/** + * @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"); + + uint8_t m = rows; + uint8_t n = columns; + uint8_t p = (m < n) ? m : n; // rank = min(m,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 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; + } + } + } + + // 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; + } + + // ---- Phase 1: Householder Bidiagonalization ---- + // Reduce W to upper bidiagonal form using Householder reflections. + // For a p×q matrix (p ≤ q after transpose), we do p steps: + // Step k: zero out subdiagonal in column k, then superdiagonal in row k + + float hhVec[5]; // Householder vector storage + + for (uint8_t k = 0; k < p; k++) { + // --- Left Householder on column k, rows k..min(m,n)-1 --- + { + uint8_t len = (m > n) ? m - k : n - k; + if (len <= 1) + continue; + + // Extract the column segment + float x[5]; + for (uint8_t i = 0; i < len; i++) { + x[i] = W[k + i][k]; + } + + // Compute Householder reflection: H·x = [α, 0, ..., 0]ᵀ + float norm = 0; + for (uint8_t i = 0; i < len; i++) + norm += x[i] * x[i]; + norm = sqrtf(norm); + + if (norm < 1e-30f) + continue; + + float alpha = (x[0] >= 0) ? -norm : norm; + float v0 = x[0] - alpha; + float vv = v0 * v0 + norm * norm - alpha * x[0]; + if (vv < 1e-30f) + continue; + + float scale = 1.0f / sqrtf(vv); + // Store Householder vector (first element is implicit 1, rest in hhVec) + hhVec[0] = + v0 * scale; // this is the first element of the reflected vector + for (uint8_t i = 1; i < len; i++) { + hhVec[i] = x[i] * scale; + } + + // Apply H from left to W: W = H·W + // For each column j, w[k+i][j] -= 2·v_i·(vᵀ·w_col) / (vᵀv) + float vvNorm = 1.0f + hhVec[0] * hhVec[0]; + for (uint8_t i = 1; i < len; i++) { + vvNorm += hhVec[i] * hhVec[i]; + } + + for (uint8_t j = k; j < n; j++) { + float dot = 0; + for (uint8_t i = 0; i < len; i++) { + dot += hhVec[i] * W[k + i][j]; + } + dot *= 2.0f / vvNorm; + for (uint8_t i = 0; i < len; i++) { + W[k + i][j] -= dot * hhVec[i]; + } + } + + // Apply H from right to QL: QL = QL · H + for (uint8_t j = k; j < m; j++) { + float dot = 0; + for (uint8_t i = 0; i < len; i++) { + dot += hhVec[i] * QL[j][k + i]; + } + dot *= 2.0f / vvNorm; + for (uint8_t i = 0; i < len; i++) { + QL[j][k + i] -= dot * hhVec[i]; + } + } + } + + // --- Right Householder on row k, columns k+1..min(m,n)-1 --- + { + uint8_t len = (p > 1) ? p - 1 - k : 0; + if (len <= 0) + continue; + + // Extract the row segment + float x[5]; + for (uint8_t i = 0; i < len; i++) { + x[i] = W[k][k + 1 + i]; + } + + // Compute Householder reflection + float norm = 0; + for (uint8_t i = 0; i < len; i++) + norm += x[i] * x[i]; + norm = sqrtf(norm); + + if (norm < 1e-30f) + continue; + + float alpha = (x[0] >= 0) ? -norm : norm; + float v0 = x[0] - alpha; + float vv = v0 * v0 + norm * norm - alpha * x[0]; + if (vv < 1e-30f) + continue; + + float scale = 1.0f / sqrtf(vv); + hhVec[0] = v0 * scale; + for (uint8_t i = 1; i < len; i++) { + hhVec[i] = x[i] * scale; + } + + // Compute vᵀv + float vvNorm = 1.0f + hhVec[0] * hhVec[0]; + for (uint8_t i = 1; i < len; i++) { + vvNorm += hhVec[i] * hhVec[i]; + } + + // Apply H from right to W: W = W·H + for (uint8_t i = 0; i < m; i++) { + float dot = 0; + for (uint8_t j = 0; j < len; j++) { + dot += hhVec[j] * W[i][k + 1 + j]; + } + dot *= 2.0f / vvNorm; + for (uint8_t j = 0; j < len; j++) { + W[i][k + 1 + j] -= dot * hhVec[j]; + } + } + + // Apply H from right to QR: QR = QR · H + for (uint8_t i = 0; i < n; i++) { + float dot = 0; + for (uint8_t j = 0; j < len; j++) { + dot += hhVec[j] * QR[i][k + 1 + j]; + } + dot *= 2.0f / vvNorm; + for (uint8_t j = 0; j < len; j++) { + QR[i][k + 1 + j] -= dot * hhVec[j]; + } + } + } + } + + // ---- Phase 2: Implicit QR Iteration on Bidiagonal Matrix ---- + // W now contains the upper bidiagonal matrix B. + // We apply implicit QR steps to diagonalize it. + + uint32_t maxIter = 1000; + float tol = 1e-8f; + + for (uint32_t iter = 0; iter < maxIter; iter++) { + // Deflate: zero out negligible subdiagonal elements + for (uint8_t i = p - 1; i > 0; i--) { + float test = fabsf(W[i][i - 1]); + float scale = fabsf(W[i - 1][i - 1]) + fabsf(W[i][i]); + if (test < tol * (scale + 1e-30f)) { + W[i][i - 1] = 0; + } + } + + // Find the smallest unreduced block [start..end] + uint8_t start = 0, end = p - 1; + for (uint8_t i = 0; i < p - 1; i++) { + if (fabsf(W[i + 1][i]) > + tol * (fabsf(W[i][i]) + fabsf(W[i + 1][i + 1]) + 1e-30f)) { + start = i + 1; + } + } + for (int8_t i = (int8_t)p - 2; i >= 0; i--) { + if (fabsf(W[i + 1][i]) > + tol * (fabsf(W[i][i]) + fabsf(W[i + 1][i + 1]) + 1e-30f)) { + end = (uint8_t)i; + break; + } + } + + // Check convergence of the block + if (start >= end) { + continue; + } + + bool blockConverged = true; + for (uint8_t i = start; i <= end; i++) { + if (i > start && fabsf(W[i][i - 1]) > tol * (fabsf(W[i - 1][i - 1]) + + fabsf(W[i][i]) + 1e-30f)) { + blockConverged = false; + break; + } + if (i < end && + fabsf(W[i][i + 1]) > + tol * (fabsf(W[i][i]) + fabsf(W[i + 1][i + 1]) + 1e-30f)) { + blockConverged = false; + break; + } + } + if (blockConverged) + continue; + + // Wilkinson shift from bottom 2×2 corner + float a = W[end - 1][end - 1]; + float b = W[end - 1][end]; + float c = W[end][end - 1]; + float d = W[end][end]; + + float trace = a + d; + float det = a * d - b * c; + float disc = trace * trace - 4.0f * det; + float shift; + if (disc >= 0) { + float sqrtDisc = sqrtf(disc); + float e1 = (trace + sqrtDisc) / 2.0f; + float e2 = (trace - sqrtDisc) / 2.0f; + shift = (fabsf(e1 - d) < fabsf(e2 - d)) ? e1 : e2; + } else { + shift = d; + } + + // --- Implicit QR step using Givens rotations --- + // First, apply Givens rotation from the left to zero out (W[start][start-1] + // - shift) For the bidiagonal structure, we process from top to bottom. + + float x = W[start][start] - shift; + float y = (start > 0) ? W[start][start - 1] : 0.0f; + + for (uint8_t i = start; i <= end; i++) { + float r = sqrtf(x * x + y * y); + if (r < 1e-30f) { + x = W[i][i]; + y = (i < end) ? W[i + 1][i] : 0.0f; + continue; + } + + float cs = x / r; + float sn = y / r; + + // Apply Givens from left to rows i, i+1 of W (columns i..p-1) + for (uint8_t j = i; j < p; j++) { + float t1 = W[i][j]; + float t2 = W[i + 1][j]; + W[i][j] = cs * t1 + sn * t2; + W[i + 1][j] = -sn * t1 + cs * t2; + } + + // Apply Givens from right to columns i, i+1 of W (rows 0..i) + if (i > start) { + for (uint8_t j = 0; j <= i; j++) { + float t1 = W[j][i]; + float t2 = W[j][i + 1]; + W[j][i] = cs * t1 + sn * t2; + W[j][i + 1] = -sn * t1 + cs * t2; + } + } + + // Accumulate right transformations into QR + for (uint8_t j = 0; j < n; j++) { + float t1 = QR[j][i]; + float t2 = QR[j][i + 1]; + QR[j][i] = cs * t1 + sn * t2; + QR[j][i + 1] = -sn * t1 + cs * t2; + } + + // Prepare next Givens rotation + x = W[i + 1][i]; + y = (i + 1 < end) ? W[i + 1][i + 1] : 0.0f; + } + } + + // ---- Phase 3: Extract Results ---- + // Singular values are the absolute values of diagonal elements of W + for (uint8_t i = 0; i < p; i++) { + sigma[i][0] = fabsf(W[i][i]); + } + + // Sort singular values in descending order and reorder U, V accordingly + for (uint8_t i = 0; i < p - 1; i++) { + for (uint8_t j = i + 1; j < p; j++) { + if (sigma[j][0] > sigma[i][0]) { + float tmpS = sigma[i][0]; + sigma[i][0] = sigma[j][0]; + sigma[j][0] = tmpS; + + // 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 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; + } + } + } + } + + // ---- Phase 4: Compute Final U and Vt ---- + // If transposeNeeded (wide matrix), we computed SVD(Aᵀ) = Ũ·Σ·Ṽᵀ + // Then U = Ṽ (= QR[:,0:p]) and Vt = Ũᵀ (= QL[:,0:p]ᵀ) + // Otherwise, SVD(A) = QL[:,0:p] · Σ · (QR[:,0:p])ᵀ + // So U = QL[:,0:p] and Vt = QR[:,0:p]ᵀ + + for (uint8_t i = 0; i < m; i++) { + for (uint8_t j = 0; j < n; j++) { + if (j < p) { + if (transposeNeeded) { + // U = QR[:, 0:p]ᵀ → U[i][j] = QR[j][i] + U[i][j] = QR[j][i]; + } else { + // U = QL[:, 0:p] + U[i][j] = QL[i][j]; + } + } else { + U[i][j] = 0; + } + } + } + + for (uint8_t i = 0; i < n; i++) { + for (uint8_t j = 0; j < m; j++) { + if (i < p && j < m) { + if (transposeNeeded) { + // Vt = QL[:, 0:p]ᵀ → Vt[i][j] = QL[j][i] + Vt[i][j] = QL[j][i]; + } else { + // Vt = QR[:, 0:p]ᵀ → Vt[i][j] = QR[j][i] + Vt[i][j] = QR[j][i]; + } + } else { + Vt[i][j] = 0; + } + } + } +} + +#endif \ No newline at end of file diff --git a/src/SVD.hpp b/src/SVD.hpp new file mode 100644 index 0000000..922a230 --- /dev/null +++ b/src/SVD.hpp @@ -0,0 +1,131 @@ +#pragma once +#include "Matrix.hpp" + +/** + * @brief library that uses Matrix.hpp and performs SVD on a matrix + */ +namespace SVD { +/** + * @brief Compute the Singular Value Decomposition (SVD) of this matrix. + * + * Decomposes A into U × Σ × Vᵀ where: + * - U is an m×k orthogonal matrix (left singular vectors) + * - Σ is a k×k diagonal matrix with non-negative singular values + * (stored as a k×1 column vector) + * - Vᵀ is a k×n orthogonal matrix (right singular vectors, transposed) + * - k = min(m, n) + * + * 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) + * + * @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 + */ +template +void SVD(Matrix &matrixToDecompose, Matrix &U, + Matrix &sigma, Matrix &Vt); + +// ======================================================================== +// SVD Building Block Functions (for unit testing) +// These operate on internal 5×5 working arrays for maximum flexibility. +// ======================================================================== + +/** + * @brief Compute a Householder reflector vector. + * + * 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 len Number of valid elements in x + * @param v Output: normalized Householder vector (v[0] is the first element) + * @param alpha Output: the resulting first element after reflection + * @return The norm of the input vector x + */ +static float ComputeHouseholder(const float *x, uint8_t len, float *v, + float &alpha); + +/** + * @brief Apply a Householder reflection from the left. + * + * Transforms W = (I - 2·v·vᵀ) · W where v operates on rows [startRow..endRow]. + * + * @param W Input/output: matrix to transform (5×5 working array) + * @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, + 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]. + * + * @param W Input/output: matrix to transform (5×5 working array) + * @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, + uint8_t startCol, uint8_t endCol); + +/** + * @brief Compute a Givens rotation that zeros out y. + * + * Computes c, s such that: + * [c s] [x] = [r] + * [-s c] [y] [0] + * where r = sqrt(x² + y²). + * + * @param x First element + * @param y Second element (to be zeroed) + * @param c Output: cosine of rotation angle + * @param s Output: sine of rotation angle + */ +static void ComputeGivens(float x, float y, float &c, float &s); + +/** + * @brief Apply a Givens rotation from the left to rows i and j. + * + * Applies [c s; -s c] to rows i, j of W (columns startCol..endCol). + * + * @param W Input/output: matrix to transform + * @param i First row index + * @param j Second row index + * @param c Cosine of rotation angle + * @param s Sine of rotation angle + * @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, + float s, uint8_t startCol, uint8_t endCol); + +/** + * @brief Apply a Givens rotation from the right to columns i and j. + * + * Applies [c -s; s c]ᵀ to columns i, j of W (rows startRow..endRow). + * + * @param W Input/output: matrix to transform + * @param i First column index + * @param j Second column index + * @param c Cosine of rotation angle + * @param s Sine of rotation angle + * @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, + 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 diff --git a/unit-tests/CMakeLists.txt b/unit-tests/CMakeLists.txt index a57babd..3ecc1b1 100644 --- a/unit-tests/CMakeLists.txt +++ b/unit-tests/CMakeLists.txt @@ -32,4 +32,24 @@ target_link_libraries(vector-3d-tests PRIVATE vector-3d Catch2::Catch2WithMain +) + +# SVD building block tests +add_executable(svd-build-blocks-tests svd-build-blocks-tests.cpp) + +target_link_libraries(svd-build-blocks-tests + PRIVATE + matrix + svd + Catch2::Catch2WithMain +) + +# SVD integration tests +add_executable(svd-integration-test svd-integration-test.cpp) + +target_link_libraries(svd-integration-test + PRIVATE + matrix + svd + Catch2::Catch2WithMain ) \ No newline at end of file diff --git a/unit-tests/matrix-tests.cpp b/unit-tests/matrix-tests.cpp index d9d67d4..1e125b5 100644 --- a/unit-tests/matrix-tests.cpp +++ b/unit-tests/matrix-tests.cpp @@ -4,6 +4,7 @@ // include the module you're going to test next #include "Matrix.hpp" +#include "SVD.hpp" // any other libraries #include @@ -637,4 +638,399 @@ TEST_CASE("Eigenvalues and Vectors", "Matrix") { REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(0.0f, 1e-4f)); REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(16.1168f, 1e-4f)); } +} + +// ============================================================================ +// SVD Tests — Reference values computed via scipy.linalg.svd (Python) +// ============================================================================ + +/** + * @brief Helper: compute Frobenius norm of a matrix. + */ +template +static float frobeniusNorm(const Matrix &M) { + float sum = 0; + for (uint8_t i = 0; i < rows; i++) { + for (uint8_t j = 0; j < columns; j++) { + float v = M.Get(i, j); + sum += v * v; + } + } + return sqrtf(sum); +} + +/** + * @brief Helper: compute reconstruction error ||A - UΣVᵀ||_F. + * + * Verifies the fundamental SVD identity A = U × diag(σ) × Vᵀ. + * For non-square matrices, only the first min(rows,cols) singular values + * contribute to the reconstruction. + */ +template +static float svdReconstructionError(const Matrix &A, + const Matrix &U, + const Matrix &sigma, + const Matrix &Vt) { + // Compute U × diag(σ): only first min(rows,cols) columns of U are used + constexpr uint8_t k = (rows < columns) ? rows : columns; + Matrix USigma{0}; + for (uint8_t i = 0; i < rows; i++) { + for (uint8_t j = 0; j < k; j++) { + USigma[i][j] = U.Get(i, j) * sigma.Get(j, 0); + } + } + + // Compute (UΣ) × Vᵀ: only first k rows of Vt are used + Matrix UVt{0}; + for (uint8_t i = 0; i < rows; i++) { + for (uint8_t j = 0; j < columns; j++) { + float sum = 0; + for (uint8_t p = 0; p < k; p++) { + sum += USigma[i][p] * Vt.Get(p, j); + } + UVt[i][j] = sum; + } + } + + // Compute ||A - UVᵀ||_F + Matrix diff{0}; + A.Sub(UVt, diff); + return frobeniusNorm(diff); +} + +/** + * @brief Helper: check orthogonality of the first k columns of M. + * Verifies M[:,0:k]ᵀ × M[:,0:k] ≈ I_k. + */ +template +static float orthogonalityError(const Matrix &M) { + constexpr uint8_t k = (rows < columns) ? rows : columns; + + // Compute Mᵀ × M (should be I_k in top-left) + Matrix Mt = M.Transpose(); + Matrix MtM{0}; + Mt.Mult(M, MtM); + + float err = 0; + for (uint8_t i = 0; i < k; i++) { + for (uint8_t j = 0; j < k; j++) { + float expected = (i == j) ? 1.0f : 0.0f; + err += (MtM.Get(i, j) - expected) * (MtM.Get(i, j) - expected); + } + } + return sqrtf(err); +} + +/** + * @brief Helper: check that singular values are sorted in descending order. + */ +template +static bool isSortedDescending(const Matrix &sigma, uint8_t count) { + for (uint8_t i = 0; i < count - 1; i++) { + if (sigma.Get(i + 1, 0) > sigma.Get(i, 0) + 1e-6f) { + return false; + } + } + return true; +} + +TEST_CASE("SVD: Simple 2x2 Matrix", "Matrix") { + // Reference: scipy.linalg.svd([[1,2],[3,4]]) + // σ = [5.4649857042, 0.3659661906] + Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f}; + Matrix<2, 2> U{}, Vt{}; + Matrix<2, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + // Verify singular values (verified with Python scipy.linalg.svd) + REQUIRE_THAT(sigma.Get(0, 0), + Catch::Matchers::WithinRel(5.4649857042f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), + Catch::Matchers::WithinRel(0.3659661906f, 1e-4f)); + + // Verify descending order + REQUIRE(isSortedDescending(sigma, 2)); + + // Verify U is orthogonal: UᵀU ≈ I + REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinRel(0.0f, 1e-4f)); + + // Verify Vt is orthogonal: VtVᵀ ≈ I + REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinRel(0.0f, 1e-4f)); + + // Verify reconstruction: A ≈ U Σ Vᵀ + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f)); +} + +TEST_CASE("SVD: Symmetric Positive Definite 2x2", "Matrix") { + // Reference: scipy.linalg.svd([[5,3],[3,5]]) + // σ = [8.0, 2.0] (eigenvalues since symmetric PD) + Matrix<2, 2> A{5.0f, 3.0f, 3.0f, 5.0f}; + Matrix<2, 2> U{}, Vt{}; + Matrix<2, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(8.0f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(2.0f, 1e-4f)); + + // For symmetric PD matrices, U ≈ V (up to sign) + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f)); +} + +TEST_CASE("SVD: Full-Rank 3x3 Matrix", "Matrix") { + // Reference: scipy.linalg.svd([[1,2,3],[4,5,6],[7,8,10]]) + // σ = [17.4125051668, 0.8751613501, 0.1968665211] + Matrix<3, 3> A{1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 10.0f}; + Matrix<3, 3> U{}, Vt{}; + Matrix<3, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE_THAT(sigma.Get(0, 0), + Catch::Matchers::WithinRel(17.4125051668f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), + Catch::Matchers::WithinRel(0.8751613501f, 1e-4f)); + REQUIRE_THAT(sigma.Get(2, 0), + Catch::Matchers::WithinRel(0.1968665211f, 1e-4f)); + + REQUIRE(isSortedDescending(sigma, 3)); + REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinRel(0.0f, 1e-4f)); + REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinRel(0.0f, 1e-4f)); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-3f)); +} + +TEST_CASE("SVD: Rank-Deficient 3x3 Matrix", "Matrix") { + // Reference: scipy.linalg.svd([[1,2,3],[4,5,6],[7,8,9]]) + // σ = [16.8481033526, 1.0683695146, ~0] (rank 2) + Matrix<3, 3> A{1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f}; + Matrix<3, 3> U{}, Vt{}; + Matrix<3, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE_THAT(sigma.Get(0, 0), + Catch::Matchers::WithinRel(16.8481033526f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), + Catch::Matchers::WithinRel(1.0683695146f, 1e-4f)); + // Third singular value should be ~0 (rank deficiency) + REQUIRE(sigma.Get(2, 0) < 1e-3f); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-3f)); +} + +TEST_CASE("SVD: Diagonal 3x3 Matrix", "Matrix") { + // For a diagonal matrix, σ = diagonal entries, U = V = I + Matrix<3, 3> A{10.0f, 0.0f, 0.0f, 5.0f, 0.0f, 0.0f, 0.0f, 0.0f, 2.0f}; + Matrix<3, 3> U{}, Vt{}; + Matrix<3, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(10.0f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(5.0f, 1e-4f)); + REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(2.0f, 1e-4f)); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f)); +} + +TEST_CASE("SVD: Tall Matrix (4×3)", "Matrix") { + // Reference: scipy.linalg.svd with full_matrices=False + // σ = [25.4624074360, 1.2906616758, ~0] (rank 2) + Matrix<4, 3> A{1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, + 7.0f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f}; + Matrix<4, 3> U{}; + Matrix<3, 3> Vt{}; // Vt is always n×n + Matrix<3, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE_THAT(sigma.Get(0, 0), + Catch::Matchers::WithinRel(25.4624074360f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), + Catch::Matchers::WithinRel(1.2906616758f, 1e-4f)); + REQUIRE(sigma.Get(2, 0) < 1e-3f); + + // U should be 4×3 with orthonormal columns + REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinRel(0.0f, 1e-3f)); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-3f)); +} + +TEST_CASE("SVD: Wide Matrix (3×5)", "Matrix") { + // Reference: scipy.linalg.svd with full_matrices=False + // σ = [35.1272233336, 2.4653966969, ~0] (rank 2) + Matrix<3, 5> A{1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, + 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f}; + Matrix<3, 5> U{}; + Matrix<5, 5> Vt{}; // Vt is always n×n + Matrix<5, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE_THAT(sigma.Get(0, 0), + Catch::Matchers::WithinRel(35.1272233336f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), + Catch::Matchers::WithinRel(2.4653966969f, 1e-4f)); + REQUIRE(sigma.Get(2, 0) < 1e-3f); + + // Vt should be 5×5 with orthonormal rows (first k) + REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinRel(0.0f, 1e-3f)); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-3f)); +} + +TEST_CASE("SVD: 5×5 Symmetric Tridiagonal", "Matrix") { + // Reference: scipy.linalg.svd for discrete Laplacian-like matrix + // σ = [3.7320508076, 3.0, 2.0, 1.0, 0.2679491924] + Matrix<5, 5> A{2.0f, -1.0f, 0.0f, 0.0f, 0.0f, -1.0f, 2.0f, -1.0f, 0.0f, + 0.0f, 0.0f, -1.0f, 2.0f, -1.0f, 0.0f, 0.0f, 0.0f, -1.0f, + 2.0f, -1.0f, 0.0f, 0.0f, 0.0f, -1.0f, 2.0f}; + Matrix<5, 5> U{}, Vt{}; + Matrix<5, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE_THAT(sigma.Get(0, 0), + Catch::Matchers::WithinRel(3.7320508076f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(3.0f, 1e-4f)); + REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(2.0f, 1e-4f)); + REQUIRE_THAT(sigma.Get(3, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f)); + REQUIRE_THAT(sigma.Get(4, 0), + Catch::Matchers::WithinRel(0.2679491924f, 1e-4f)); + + REQUIRE(isSortedDescending(sigma, 5)); + REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinRel(0.0f, 1e-3f)); + REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinRel(0.0f, 1e-3f)); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-3f)); +} + +TEST_CASE("SVD: Non-Square with Negative Values (2×3)", "Matrix") { + // Reference: scipy.linalg.svd([[0.5,-0.3,0.8],[-0.2,0.7,0.1]]) + // σ = [1.0384009867, 0.6646227432] + Matrix<2, 3> A{0.5f, -0.3f, 0.8f, -0.2f, 0.7f, 0.1f}; + Matrix<2, 3> U{}; + Matrix<3, 3> Vt{}; // Vt is always n×n + Matrix<3, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE_THAT(sigma.Get(0, 0), + Catch::Matchers::WithinRel(1.0384009867f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), + Catch::Matchers::WithinRel(0.6646227432f, 1e-4f)); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f)); +} + +TEST_CASE("SVD: Near-Singular 2×2 Matrix", "Matrix") { + // Condition number ≈ 1e6 — tests numerical stability + // Reference: scipy.linalg.svd([[1,0],[0,1e-6]]) + // σ = [1.0, 1e-6] + Matrix<2, 2> A{1.0f, 0.0f, 0.0f, 1e-6f}; + Matrix<2, 2> U{}, Vt{}; + Matrix<2, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1e-6f, 1e-2f)); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-6f)); +} + +TEST_CASE("SVD: Orthogonal Matrix (3×3)", "Matrix") { + // For an orthogonal matrix, all singular values should be 1. + // Rotation matrix about z-axis by 45° + float c = sqrtf(0.5f); // cos(45°) + float s = sqrtf(0.5f); // sin(45°) + Matrix<3, 3> A{c, -s, 0.0f, s, c, 0.0f, 0.0f, 0.0f, 1.0f}; + Matrix<3, 3> U{}, Vt{}; + Matrix<3, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + // All singular values should be 1 for an orthogonal matrix + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f)); + REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f)); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f)); +} + +TEST_CASE("SVD: Identity Matrix", "Matrix") { + // For I, σ = [1, 1, 1], U = V = I + Matrix<3, 3> A{1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f}; + Matrix<3, 3> U{}, Vt{}; + Matrix<3, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f)); + REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f)); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-6f)); +} + +TEST_CASE("SVD: Zero Matrix", "Matrix") { + // All singular values should be zero + Matrix<3, 3> A{0.0f}; + Matrix<3, 3> U{}, Vt{}; + Matrix<3, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE(sigma.Get(0, 0) < 1e-6f); + REQUIRE(sigma.Get(1, 0) < 1e-6f); + REQUIRE(sigma.Get(2, 0) < 1e-6f); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-6f)); +} + +TEST_CASE("SVD: 2×1 Column Vector", "Matrix") { + // For a column vector v, σ = ||v||, U = v/||v|| (with padding) + Matrix<2, 1> A{3.0f, 4.0f}; + Matrix<2, 1> U{}; + Matrix<1, 1> Vt{}; // Vt is always n×n + Matrix<1, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + // σ should be the Euclidean norm: ||[3,4]|| = 5 + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.0f, 1e-4f)); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f)); +} + +TEST_CASE("SVD: 1×2 Row Vector", "Matrix") { + // For a row vector vᵀ, σ = ||v||, Vt = v/||v|| (with padding) + Matrix<1, 2> A{3.0f, 4.0f}; + Matrix<1, 2> U{}; + Matrix<2, 2> Vt{}; // Vt is always n×n + Matrix<2, 1> sigma{}; + + SVD::SVD(A, U, sigma, Vt); + + // σ should be the Euclidean norm: ||[3,4]|| = 5 + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.0f, 1e-4f)); + + float reconErr = svdReconstructionError(A, U, sigma, Vt); + REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f)); } \ No newline at end of file diff --git a/unit-tests/svd-build-blocks-tests.cpp b/unit-tests/svd-build-blocks-tests.cpp new file mode 100644 index 0000000..857ee54 --- /dev/null +++ b/unit-tests/svd-build-blocks-tests.cpp @@ -0,0 +1,785 @@ +// include the unit test framework first +#include +#include + +// include the module you're going to test next +#include "Matrix.hpp" +#include "SVD.hpp" + +// any other libraries +#include +#include +#include + +// ============================================================================ +// Helper: Frobenius norm of a 5×5 matrix +// ============================================================================ +static float frobeniusNorm5(const Matrix<5, 5> &M) { + float sum = 0.0f; + for (uint8_t i = 0; i < 5; i++) { + for (uint8_t j = 0; j < 5; j++) { + float v = M.Get(i, j); + sum += v * v; + } + } + return sqrtf(sum); +} + +// ============================================================================ +// Helper: Check if a matrix is orthogonal (Mᵀ·M ≈ I) +// ============================================================================ +static bool isOrthogonal5(const Matrix<5, 5> &M, float tol = 1e-6f) { + Matrix<5, 5> Mt = M.Transpose(); + Matrix<5, 5> MtM{0}; + Mt.Mult(M, MtM); + + for (uint8_t i = 0; i < 5; i++) { + for (uint8_t j = 0; j < 5; j++) { + float expected = (i == j) ? 1.0f : 0.0f; + if (fabsf(MtM.Get(i, j) - expected) > tol) { + return false; + } + } + } + return true; +} + +// ============================================================================ +// TEST 1: ComputeHouseholder +// ============================================================================ +TEST_CASE("SVD Building Block: ComputeHouseholder", "[Matrix][SVD]") { + // Test case: [3, 4] should give alpha = -5 (norm), v normalized + // Reference: scipy.linalg.householder([3, 4]) → v ≈ [0.894427191, + // 0.447213596], α = -5 + { + float x[] = {3.0f, 4.0f}; + float v[5] = {0}; + float alpha = 0; + + float norm = SVD::ComputeHouseholder(x, 2, v, alpha); + + // Norm should be 5.0 + REQUIRE_THAT(norm, Catch::Matchers::WithinRel(5.0f, 1e-6f)); + + // Alpha should be -5 (negative norm) + REQUIRE_THAT(alpha, Catch::Matchers::WithinRel(-5.0f, 1e-6f)); + + // v should be normalized: ||v|| ≈ 1 + float vNorm = sqrtf(v[0] * v[0] + v[1] * v[1]); + REQUIRE_THAT(vNorm, Catch::Matchers::WithinRel(1.0f, 1e-6f)); + + // Verify H·x = [alpha, 0]: (I - 2vvᵀ)·x should give [-5, 0] + float hx0 = x[0] - 2.0f * v[0] * (v[0] * x[0] + v[1] * x[1]); + float hx1 = x[1] - 2.0f * v[1] * (v[0] * x[0] + v[1] * x[1]); + REQUIRE_THAT(hx0, Catch::Matchers::WithinRel(alpha, 1e-6f)); + REQUIRE_THAT(hx1, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + + // Test case: [1, 3] + // Reference: norm = √10 ≈ 3.16228, alpha = -√10 + { + float x[] = {1.0f, 3.0f}; + float v[5] = {0}; + float alpha = 0; + + float norm = SVD::ComputeHouseholder(x, 2, v, alpha); + + REQUIRE_THAT(norm, Catch::Matchers::WithinRel(sqrtf(10.0f), 1e-6f)); + REQUIRE_THAT(alpha, Catch::Matchers::WithinRel(-sqrtf(10.0f), 1e-6f)); + + // Verify H·x = [alpha, 0] + float dot = v[0] * x[0] + v[1] * x[1]; + float hx0 = x[0] - 2.0f * v[0] * dot; + float hx1 = x[1] - 2.0f * v[1] * dot; + REQUIRE_THAT(hx0, Catch::Matchers::WithinRel(alpha, 1e-6f)); + REQUIRE_THAT(hx1, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + + // Test case: [1, 2, 3] (3D) + // Reference: norm = √14 ≈ 3.74166 + { + float x[] = {1.0f, 2.0f, 3.0f}; + float v[5] = {0}; + float alpha = 0; + + float norm = SVD::ComputeHouseholder(x, 3, v, alpha); + + REQUIRE_THAT(norm, Catch::Matchers::WithinRel(sqrtf(14.0f), 1e-6f)); + REQUIRE_THAT(alpha, Catch::Matchers::WithinRel(-sqrtf(14.0f), 1e-6f)); + + // Verify v is normalized + float vNorm = sqrtf(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + REQUIRE_THAT(vNorm, Catch::Matchers::WithinRel(1.0f, 1e-6f)); + + // Verify H·x = [alpha, 0, 0] + float dot = v[0] * x[0] + v[1] * x[1] + v[2] * x[2]; + for (uint8_t i = 0; i < 3; i++) { + float hx_i = x[i] - 2.0f * v[i] * dot; + if (i == 0) { + REQUIRE_THAT(hx_i, Catch::Matchers::WithinRel(alpha, 1e-6f)); + } else { + REQUIRE_THAT(hx_i, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + } + } + + // Test case: [0, 0, 1] (already has leading zeros) + { + float x[] = {0.0f, 0.0f, 1.0f}; + float v[5] = {0}; + float alpha = 0; + + float norm = SVD::ComputeHouseholder(x, 3, v, alpha); + + REQUIRE_THAT(norm, Catch::Matchers::WithinRel(1.0f, 1e-6f)); + REQUIRE_THAT(alpha, Catch::Matchers::WithinRel(-1.0f, 1e-6f)); + + // Verify H·x = [-1, 0, 0] + float dot = v[0] * x[0] + v[1] * x[1] + v[2] * x[2]; + float hx0 = x[0] - 2.0f * v[0] * dot; + float hx1 = x[1] - 2.0f * v[1] * dot; + float hx2 = x[2] - 2.0f * v[2] * dot; + REQUIRE_THAT(hx0, Catch::Matchers::WithinRel(alpha, 1e-6f)); + REQUIRE_THAT(hx1, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + REQUIRE_THAT(hx2, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + + // Test case: [5, -3, 2, 1] (4D) + { + float x[] = {5.0f, -3.0f, 2.0f, 1.0f}; + float v[5] = {0}; + float alpha = 0; + + float norm = SVD::ComputeHouseholder(x, 4, v, alpha); + + REQUIRE_THAT(norm, Catch::Matchers::WithinRel(sqrtf(39.0f), 1e-6f)); + REQUIRE_THAT(alpha, Catch::Matchers::WithinRel(-sqrtf(39.0f), 1e-6f)); + + // Verify H·x = [alpha, 0, 0, 0] + float dot = v[0] * x[0] + v[1] * x[1] + v[2] * x[2] + v[3] * x[3]; + for (uint8_t i = 0; i < 4; i++) { + float hx_i = x[i] - 2.0f * v[i] * dot; + if (i == 0) { + REQUIRE_THAT(hx_i, Catch::Matchers::WithinRel(alpha, 1e-6f)); + } else { + REQUIRE_THAT(hx_i, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + } + } + + // Test case: zero vector + { + float x[] = {0.0f, 0.0f}; + float v[5] = {0}; + float alpha = 0; + + float norm = SVD::ComputeHouseholder(x, 2, v, alpha); + + REQUIRE_THAT(norm, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + REQUIRE(alpha == 0.0f); + } +} + +// ============================================================================ +// TEST 2: ApplyHouseholderLeft +// ============================================================================ +TEST_CASE("SVD Building Block: ApplyHouseholderLeft", "[Matrix][SVD]") { + // Test: Apply Householder to zero out column 0, rows 1:2 of a 3×3 matrix + // Input: [[1, 2, 3], [4, 5, 6], [7, 8, 9]] + // After applying HH on col 0 (rows 1:2): A[2,0] should be ~0 + { + Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0, 0, 4.0f, 5.0f, 6.0f, 0, + 0, 7.0f, 8.0f, 9.0f, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0}; + + // Compute Householder for column 0, rows 1:2 → vector [4, 7] + float x[] = {4.0f, 7.0f}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 2, v, alpha); + + // Apply from left + SVD::ApplyHouseholderLeft(W, v, 1, 2); + + // A[2,0] should be ~0 + REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + + // Verify orthogonality of the transformation: W = H·W_original + Matrix<5, 5> W_orig{1.0f, 2.0f, 3.0f, 0, 0, 4.0f, 5.0f, 6.0f, 0, + 0, 7.0f, 8.0f, 9.0f, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0}; + + // Compute H_left explicitly: I - 2*v*vᵀ (on rows 1:2) + Matrix<5, 5> H_left{0}; + for (uint8_t i = 0; i < 5; i++) { + H_left[i][i] = 1.0f; + } + // Apply -2*v*vᵀ to the sub-block + float vv = v[0] * v[0] + v[1] * v[1]; + for (uint8_t i = 1; i <= 2; i++) { + for (uint8_t j = 1; j <= 2; j++) { + H_left[i][j] -= 2.0f * v[i - 1] * v[j - 1] / vv; + } + } + + // Verify: W ≈ H_left · W_orig + Matrix<5, 5> HLeftW{0}; + H_left.Mult(W_orig, HLeftW); + float err = frobeniusNorm5(W - HLeftW); + REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-4f)); + + // Verify H_left is orthogonal + REQUIRE(isOrthogonal5(H_left)); + } + + // Test: Apply to a larger block (4 rows) + { + Matrix<5, 5> W{1.0f, 2.0f, 0, 0, 0, 3.0f, 4.0f, 0, 0, + 0, 5.0f, 6.0f, 0, 0, 0, 7.0f, 8.0f, 0, + 0, 0, 0, 0, 0, 0, 0}; + + // Householder on [3, 5, 7] (rows 1:3) + float x[] = {3.0f, 5.0f, 7.0f}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 3, v, alpha); + + SVD::ApplyHouseholderLeft(W, v, 1, 3); + + // A[2,0] and A[3,0] should be ~0 + REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + REQUIRE_THAT(W.Get(3, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } +} + +// ============================================================================ +// TEST 3: ApplyHouseholderRight +// ============================================================================ +TEST_CASE("SVD Building Block: ApplyHouseholderRight", "[Matrix][SVD]") { + // Test: Apply Householder to zero out row 0, cols 1:2 of a 3×3 matrix + // Input: [[1, 2, 3], [4, 5, 6], [7, 8, 9]] + // After applying HH on row 0 (cols 1:2): A[0,2] should be ~0 + { + Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0, 0, 4.0f, 5.0f, 6.0f, 0, + 0, 7.0f, 8.0f, 9.0f, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0}; + + // Householder for row 0, cols 1:2 → vector [2, 3] + float x[] = {2.0f, 3.0f}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 2, v, alpha); + + // Apply from right + SVD::ApplyHouseholderRight(W, v, 1, 2); + + // A[0,2] should be ~0 + REQUIRE_THAT(W.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + + // Verify W ≈ W_orig · H_right + Matrix<5, 5> W_orig{1.0f, 2.0f, 3.0f, 0, 0, 4.0f, 5.0f, 6.0f, 0, + 0, 7.0f, 8.0f, 9.0f, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0}; + + // Compute H_right = I - 2*v*vᵀ (on cols 1:2) + Matrix<5, 5> H_right{0}; + for (uint8_t i = 0; i < 5; i++) { + H_right[i][i] = 1.0f; + } + float vv = v[0] * v[0] + v[1] * v[1]; + for (uint8_t i = 1; i <= 2; i++) { + for (uint8_t j = 1; j <= 2; j++) { + H_right[i][j] -= 2.0f * v[i - 1] * v[j - 1] / vv; + } + } + + Matrix<5, 5> WOrigH{0}; + W_orig.Mult(H_right, WOrigH); + float err = frobeniusNorm5(W - WOrigH); + REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-4f)); + + // Verify H_right is orthogonal + REQUIRE(isOrthogonal5(H_right)); + } + + // Test: Apply to wider block (4 cols) + { + Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 4.0f, 0, 5.0f, 6.0f, 7.0f, 8.0f, + 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0}; + + // Householder on [2, 3, 4] (cols 1:3) + float x[] = {2.0f, 3.0f, 4.0f}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 3, v, alpha); + + SVD::ApplyHouseholderRight(W, v, 1, 3); + + // A[0,2] and A[0,3] should be ~0 + REQUIRE_THAT(W.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + REQUIRE_THAT(W.Get(0, 3), Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } +} + +// ============================================================================ +// TEST 4: ComputeGivens +// ============================================================================ +TEST_CASE("SVD Building Block: ComputeGivens", "[Matrix][SVD]") { + // Test case: [3, 4] → c = 0.6, s = 0.8 (3-4-5 triangle) + { + float c, s; + SVD::ComputeGivens(3.0f, 4.0f, c, s); + + REQUIRE_THAT(c, Catch::Matchers::WithinRel(0.6f, 1e-6f)); + REQUIRE_THAT(s, Catch::Matchers::WithinRel(0.8f, 1e-6f)); + + // Verify: [c s; -s c] · [3; 4] = [5; 0] + float r = c * 3.0f + s * 4.0f; + float z = -s * 3.0f + c * 4.0f; + REQUIRE_THAT(r, Catch::Matchers::WithinRel(5.0f, 1e-6f)); + REQUIRE_THAT(z, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + + // Verify c² + s² = 1 + REQUIRE_THAT(c * c + s * s, Catch::Matchers::WithinRel(1.0f, 1e-6f)); + } + + // Test case: [1, 0] → c = 1, s = 0 + { + float c, s; + SVD::ComputeGivens(1.0f, 0.0f, c, s); + REQUIRE_THAT(c, Catch::Matchers::WithinRel(1.0f, 1e-6f)); + REQUIRE_THAT(s, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + + // Test case: [0, 5] → c = 0, s = 1 + { + float c, s; + SVD::ComputeGivens(0.0f, 5.0f, c, s); + REQUIRE_THAT(c, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + REQUIRE_THAT(s, Catch::Matchers::WithinRel(1.0f, 1e-6f)); + + // Verify: [c s; -s c] · [0; 5] = [5; 0] + float r = c * 0.0f + s * 5.0f; + float z = -s * 0.0f + c * 5.0f; + REQUIRE_THAT(r, Catch::Matchers::WithinRel(5.0f, 1e-6f)); + REQUIRE_THAT(z, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + + // Test case: [-3, -4] → c = -0.6, s = -0.8 + { + float c, s; + SVD::ComputeGivens(-3.0f, -4.0f, c, s); + REQUIRE_THAT(c, Catch::Matchers::WithinRel(-0.6f, 1e-6f)); + REQUIRE_THAT(s, Catch::Matchers::WithinRel(-0.8f, 1e-6f)); + + // Verify: [c s; -s c] · [-3; -4] = [5; 0] + float r = c * (-3.0f) + s * (-4.0f); + float z = -s * (-3.0f) + c * (-4.0f); + REQUIRE_THAT(r, Catch::Matchers::WithinRel(5.0f, 1e-6f)); + REQUIRE_THAT(z, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + + // Test case: [1, -1] → c = 1/√2, s = -1/√2 (45°) + { + float c, s; + SVD::ComputeGivens(1.0f, -1.0f, c, s); + float invSqrt2 = 1.0f / sqrtf(2.0f); + REQUIRE_THAT(c, Catch::Matchers::WithinRel(invSqrt2, 1e-6f)); + REQUIRE_THAT(s, Catch::Matchers::WithinRel(-invSqrt2, 1e-6f)); + + // Verify: [c s; -s c] · [1; -1] = [√2; 0] + float r = c * 1.0f + s * (-1.0f); + float z = -s * 1.0f + c * (-1.0f); + REQUIRE_THAT(r, Catch::Matchers::WithinRel(sqrtf(2.0f), 1e-6f)); + REQUIRE_THAT(z, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + + // Test case: [0, 0] → c = 1, s = 0 (identity) + { + float c, s; + SVD::ComputeGivens(0.0f, 0.0f, c, s); + REQUIRE_THAT(c, Catch::Matchers::WithinRel(1.0f, 1e-6f)); + REQUIRE_THAT(s, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + + // Test case: [7, 24] → c = 7/25, s = 24/25 (7-24-25 triangle) + { + float c, s; + SVD::ComputeGivens(7.0f, 24.0f, c, s); + REQUIRE_THAT(c, Catch::Matchers::WithinRel(7.0f / 25.0f, 1e-6f)); + REQUIRE_THAT(s, Catch::Matchers::WithinRel(24.0f / 25.0f, 1e-6f)); + + float r = c * 7.0f + s * 24.0f; + float z = -s * 7.0f + c * 24.0f; + REQUIRE_THAT(r, Catch::Matchers::WithinRel(25.0f, 1e-6f)); + REQUIRE_THAT(z, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } +} + +// ============================================================================ +// TEST 5: ApplyGivensLeft +// ============================================================================ +TEST_CASE("SVD Building Block: ApplyGivensLeft", "[Matrix][SVD]") { + // Test: Apply Givens to zero out W[1,0] of a 2×2 matrix + // Input: [[3, 4], [1, 2]] + // Givens on rows 0,1 with x=W[0,0]=3, y=W[1,0]=1 + { + Matrix<5, 5> W{3.0f, 4.0f, 0, 0, 0, 1.0f, 2.0f, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + + float c, s; + SVD::ComputeGivens(3.0f, 1.0f, c, s); + + SVD::ApplyGivensLeft(W, 0, 1, c, s, 0, 4); + + // W[1,0] should be ~0 + REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + + // Verify W ≈ G · W_orig + Matrix<5, 5> W_orig{3.0f, 4.0f, 0, 0, 0, 1.0f, 2.0f, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + + // Givens rotation matrix (5×5) + Matrix<5, 5> G{0}; + for (uint8_t i = 0; i < 5; i++) { + G[i][i] = 1.0f; + } + G[0][0] = c; + G[0][1] = s; + G[1][0] = -s; + G[1][1] = c; + + Matrix<5, 5> GW{0}; + G.Mult(W_orig, GW); + float err = frobeniusNorm5(W - GW); + REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + + // Verify G is orthogonal + REQUIRE(isOrthogonal5(G)); + } + + // Test: Apply to larger range of columns + { + Matrix<5, 5> W{3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 1.0f, 2.0f, 3.0f, 4.0f, + 5.0f, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0}; + + float c, s; + SVD::ComputeGivens(3.0f, 1.0f, c, s); + + SVD::ApplyGivensLeft(W, 0, 1, c, s, 0, 4); + + REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } +} + +// ============================================================================ +// TEST 6: ApplyGivensRight +// ============================================================================ +TEST_CASE("SVD Building Block: ApplyGivensRight", "[Matrix][SVD]") { + // Test: Apply Givens to zero out W[0,1] of a 2×2 matrix + // Input: [[3, 4], [1, 2]] + // Givens on cols 0,1 with x=W[0,0]=3, y=W[0,1]=4 + { + Matrix<5, 5> W{3.0f, 4.0f, 0, 0, 0, 1.0f, 2.0f, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + + float c, s; + SVD::ComputeGivens(3.0f, 4.0f, c, s); + + SVD::ApplyGivensRight(W, 0, 1, c, s, 0, 4); + + // W[0,1] should be ~0 + REQUIRE_THAT(W.Get(0, 1), Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + + // Verify W ≈ W_orig · G + Matrix<5, 5> W_orig{3.0f, 4.0f, 0, 0, 0, 1.0f, 2.0f, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + + // Givens rotation matrix (5×5) + Matrix<5, 5> G{0}; + for (uint8_t i = 0; i < 5; i++) { + G[i][i] = 1.0f; + } + G[0][0] = c; + G[0][1] = -s; + G[1][0] = s; + G[1][1] = c; + + Matrix<5, 5> WG{0}; + W_orig.Mult(G, WG); + float err = frobeniusNorm5(W - WG); + REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + + // Verify G is orthogonal + REQUIRE(isOrthogonal5(G)); + } + + // Test: Apply to larger range of rows + { + Matrix<5, 5> W{3.0f, 4.0f, 0, 0, 0, 1.0f, 2.0f, 0, 0, 0, 5.0f, 6.0f, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + + float c, s; + SVD::ComputeGivens(3.0f, 4.0f, c, s); + + SVD::ApplyGivensRight(W, 0, 1, c, s, 0, 2); + + REQUIRE_THAT(W.Get(0, 1), Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } +} + +// ============================================================================ +// TEST 7: Full Bidiagonalization (composing Householder steps) +// ============================================================================ +TEST_CASE("SVD Building Block: Householder Bidiagonalization", + "[Matrix][SVD]") { + // Test: Bidiagonalize a 3×3 matrix and verify reconstruction + // Input: [[1, 2, 3], [4, 5, 6], [7, 8, 10]] + { + Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0, 0, 4.0f, 5.0f, 6.0f, 0, + 0, 7.0f, 8.0f, 10.0f, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0}; + + // Step 1: Left HH on column 0, rows 1:2 → zero out W[2,0] + { + float x[] = {4.0f, 7.0f}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 2, v, alpha); + SVD::ApplyHouseholderLeft(W, v, 1, 2); + } + + // Step 2: Right HH on row 0, cols 1:2 → zero out W[0,2] + { + float x[] = {W.Get(0, 1), W.Get(0, 2)}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 2, v, alpha); + SVD::ApplyHouseholderRight(W, v, 1, 2); + } + + // Step 3: Left HH on column 1, rows 2:2 → nothing to do (single element) + + // Verify bidiagonal structure: for 3x3, zero elements are A[2][0] (below + // subdiag in col 0) and A[0][2] (above superdiag in row 0) A[2][1] is the + // subdiagonal element of col 1 — valid in bidiagonal form + REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-4f)); + REQUIRE_THAT(W.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-4f)); + } + + // Test: Bidiagonalize a 4×3 matrix + { + Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0, 0, 4.0f, 5.0f, 6.0f, 0, + 0, 7.0f, 8.0f, 9.0f, 0, 0, 10.0f, 11.0f, 12.0f, + 0, 0, 0, 0, 0, 0, 0}; + + // Step 1: Left HH on col 0, rows 1:3 → zero out W[2,0], W[3,0] + { + float x[] = {4.0f, 7.0f, 10.0f}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 3, v, alpha); + SVD::ApplyHouseholderLeft(W, v, 1, 3); + } + + // Step 2: Right HH on row 0, cols 1:2 → zero out W[0,2] + { + float x[] = {W.Get(0, 1), W.Get(0, 2)}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 2, v, alpha); + SVD::ApplyHouseholderRight(W, v, 1, 2); + } + + // Step 3: Left HH on col 1, rows 2:3 → zero out W[3,1] + { + float x[] = {W.Get(2, 1), W.Get(3, 1)}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 2, v, alpha); + SVD::ApplyHouseholderLeft(W, v, 2, 3); + } + + // Verify bidiagonal structure + REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + REQUIRE_THAT(W.Get(3, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + REQUIRE_THAT(W.Get(3, 1), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + } + + // Test: Diagonal matrix (no transformations needed) + { + Matrix<5, 5> W{10.0f, 0, 0, 0, 0, 0, 5.0f, 0, 0, 0, 0, 0, 2.0f, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + + // Householder on zero vector should be identity + float x[] = {0.0f, 0.0f}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 2, v, alpha); + + // Applying identity should not change anything + Matrix<5, 5> W_copy{10.0f, 0, 0, 0, 0, 0, 5.0f, 0, 0, 0, 0, 0, 2.0f, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + SVD::ApplyHouseholderLeft(W_copy, v, 1, 2); + + REQUIRE_THAT(frobeniusNorm5(W - W_copy), + Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } +} + +// ============================================================================ +// // +// ============================================================================ +// TEST 8: Givens QR step on bidiagonal matrix +// =========================================================================== +TEST_CASE("SVD Building Block: Givens QR Step on Bidiagonal", "[Matrix][SVD]") { + // Test: Apply left Givens to zero subdiagonal of a bidiagonal matrix, + // then apply right Givens with restricted row range to restore bidiagonal + // form. + // + // Input: 3x3 bidiagonal [[1, 2, 0], [3, -4, 5], [0, 6, -7]] + // Step 1: Left Givens on rows 0,1 with x=W[0][0]=1, y=W[1][0]=3 -> zero + // W[1][0] Step 2: Right Givens on cols 1,2 with x=W[0][1], y=W[0][2] -> zero + // W[0][2] + // Only applied to row 0 (to not reintroduce subdiagonal non-zeros) + { + Matrix<5, 5> W{1.0f, 2.0f, 0.0f, 0, 0, 3.0f, -4.0f, 5.0f, 0, + 0, 0.0f, 6.0f, -7.0f, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0}; + + float c, s; + SVD::ComputeGivens(W.Get(0, 0), W.Get(1, 0), c, s); + + // Apply from left to zero subdiagonal at W[1][0] + SVD::ApplyGivensLeft(W, 0, 1, c, s, 0, 4); + + REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + + // After left Givens, W[0][2] may have become non-zero (fill-in from row 0) + // Apply right Givens to cols 1,2 with x=W[0][1], y=W[0][2] -> zero W[0][2] + // Only apply to rows 0 (to preserve bidiagonal structure below row 0) + float c2, s2; + SVD::ComputeGivens(W.Get(0, 1), W.Get(0, 2), c2, s2); + SVD::ApplyGivensRight(W, 1, 2, c2, s2, 0, 0); + + // Should be bidiagonal: W[1][0] ~ 0 (from left Givens), W[0][2] ~ 0 (from + // right Givens) + REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + REQUIRE_THAT(W.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + } + + // Test: Verify that a full QR step (left + right Givens) preserves the + // bidiagonal structure when applied correctly with proper row ranges. + { + Matrix<5, 5> W{2.0f, 3.0f, 0, 0, 0, -1.0f, 4.0f, 5.0f, 0, 0, 0, 6.0f, -7.0f, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + + // Left Givens on col 0 (rows 0,1) + float c, s; + SVD::ComputeGivens(W.Get(0, 0), W.Get(1, 0), c, s); + SVD::ApplyGivensLeft(W, 0, 1, c, s, 0, 4); + + REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + + // Right Givens on row 0 (cols 1,2) - only affect row 0 + float c2, s2; + SVD::ComputeGivens(W.Get(0, 1), W.Get(0, 2), c2, s2); + SVD::ApplyGivensRight(W, 1, 2, c2, s2, 0, 0); + + // Bidiagonal structure preserved + REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + REQUIRE_THAT(W.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + } +} + +// ============================================================================TEST +// 9: Orthogonality preservation of Householder transformations +// ============================================================================ +TEST_CASE("SVD Building Block: Householder preserves orthogonality", + "[Matrix][SVD]") { + // Starting with an orthogonal matrix, applying Householder should preserve it + { + // Identity matrix is orthogonal + Matrix<5, 5> M{0}; + for (uint8_t i = 0; i < 5; i++) { + M[i][i] = 1.0f; + } + + // Householder on first 3 elements of column 0 + float x[] = {1.0f, 0.0f, 0.0f}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 3, v, alpha); + + // Apply from left + Matrix<5, 5> M_left = M; + SVD::ApplyHouseholderLeft(M_left, v, 0, 2); + + // M_left should still be orthogonal + REQUIRE(isOrthogonal5(M_left)); + + // Apply from right + Matrix<5, 5> M_right = M; + SVD::ApplyHouseholderRight(M_right, v, 0, 2); + + REQUIRE(isOrthogonal5(M_right)); + } + + // Random orthogonal matrix (rotation) + { + float c = sqrtf(0.5f); + float s = sqrtf(0.5f); + Matrix<5, 5> M{0}; + M[0][0] = c; + M[0][1] = -s; + M[1][0] = s; + M[1][1] = c; + for (uint8_t i = 2; i < 5; i++) { + M[i][i] = 1.0f; + } + + REQUIRE(isOrthogonal5(M)); + + // Apply Householder on rows 0,1 + float x[] = {c, s}; + float v[5] = {0}; + float alpha = 0; + SVD::ComputeHouseholder(x, 2, v, alpha); + + Matrix<5, 5> M_test = M; + SVD::ApplyHouseholderLeft(M_test, v, 0, 1); + + REQUIRE(isOrthogonal5(M_test)); + } +} + +// ============================================================================ +// TEST 10: Orthogonality preservation of Givens transformations +// ============================================================================ +TEST_CASE("SVD Building Block: Givens preserves orthogonality", + "[Matrix][SVD]") { + // Starting with an orthogonal matrix, applying Givens should preserve it + { + Matrix<5, 5> M{0}; + for (uint8_t i = 0; i < 5; i++) { + M[i][i] = 1.0f; + } + + float c, s; + SVD::ComputeGivens(3.0f, 4.0f, c, s); + + // Apply from left + Matrix<5, 5> M_left = M; + SVD::ApplyGivensLeft(M_left, 0, 1, c, s, 0, 4); + + REQUIRE(isOrthogonal5(M_left)); + + // Apply from right + Matrix<5, 5> M_right = M; + SVD::ApplyGivensRight(M_right, 0, 1, c, s, 0, 4); + + REQUIRE(isOrthogonal5(M_right)); + } +} diff --git a/unit-tests/svd-integration-test.cpp b/unit-tests/svd-integration-test.cpp new file mode 100644 index 0000000..aa6051e --- /dev/null +++ b/unit-tests/svd-integration-test.cpp @@ -0,0 +1,269 @@ +#include "Matrix.hpp" +#include "SVD.hpp" +#include +#include +#include + +// Generic helper functions for any matrix size +template +static float frobeniusNorm(const Matrix &M) { + float sum = 0.0f; + for (int i = 0; i < rows; i++) + for (int j = 0; j < columns; j++) { + float v = M.Get(i, j); + sum += v * v; + } + return sqrtf(sum); +} + +template +static bool isOrthogonal(const Matrix &M, float tol = 1e-4f) { + Matrix Mt = M.Transpose(); + Matrix MtM{0}; + Mt.Mult(M, MtM); + for (int i = 0; i < n; i++) + for (int j = 0; j < n; j++) { + float expected = (i == j) ? 1.0f : 0.0f; + if (fabsf(MtM.Get(i, j) - expected) > tol) + return false; + } + return true; +} + +TEST_CASE("SVD Integration: 2x2 [[1,2],[3,4]]", "[Matrix][SVD][Integration]") { + Matrix<2, 2> A{1, 2, 3, 4}; + Matrix<2, 2> U{0}; + Matrix<2, 1> sigma{0}; + Matrix<2, 2> Vt{0}; + + SVD::SVD(A, U, sigma, Vt); + + // Reference singular values from scipy: [5.464985704219, 0.365966190626] + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.4649857f, 1e-3f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(0.3659662f, 1e-3f)); + + // Check orthogonality of U and Vt (first 2x2 blocks) + REQUIRE(isOrthogonal<2>(U)); + REQUIRE(isOrthogonal<2>(Vt)); + + // Check reconstruction: A ≈ U · diag(sigma) · Vt + Matrix<2, 2> recon{0}; + Matrix<2, 2> Usig{0}; + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + Usig[i][j] = U.Get(i, j) * sigma.Get(j, 0); + + Usig.Mult(Vt, recon); + + float err = 0.0f; + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) { + float diff = recon.Get(i, j) - A.Get(i, j); + err += diff * diff; + } + err = sqrtf(err); + REQUIRE_THAT(err, Catch::Matchers::WithinRel(0.0f, 1e-3f)); + + std::cout << "SVD 2x2 [[1,2],[3,4]]:\n"; + std::cout << "Sigma: [" << sigma.Get(0, 0) << ", " << sigma.Get(1, 0) + << "]\n"; +} + +TEST_CASE("SVD Integration: 3x3 diagonal [10,5,2]", + "[Matrix][SVD][Integration]") { + Matrix<3, 3> A{10, 0, 0, 0, 5, 0, 0, 0, 2}; + Matrix<3, 3> U{0}; + Matrix<3, 1> sigma{0}; + Matrix<3, 3> Vt{0}; + + SVD::SVD(A, U, sigma, Vt); + + // Singular values should be [10, 5, 2] (already diagonal) + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(10.0f, 1e-3f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(5.0f, 1e-3f)); + REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(2.0f, 1e-3f)); + + // U and Vt should be identity (or close) for diagonal matrix + float uErr = frobeniusNorm(U - Matrix<3, 3>{1, 0, 0, 0, 1, 0, 0, 0, 1}); + float vtErr = frobeniusNorm(Vt - Matrix<3, 3>{1, 0, 0, 0, 1, 0, 0, 0, 1}); + REQUIRE_THAT(uErr, Catch::Matchers::WithinRel(0.0f, 1e-2f)); + REQUIRE_THAT(vtErr, Catch::Matchers::WithinRel(0.0f, 1e-2f)); +} + +TEST_CASE("SVD Integration: 3x3 rank-deficient [[1,2,3],[4,5,6],[7,8,9]]", + "[Matrix][SVD][Integration]") { + Matrix<3, 3> A{1, 2, 3, 4, 5, 6, 7, 8, 9}; + Matrix<3, 3> U{0}; + Matrix<3, 1> sigma{0}; + Matrix<3, 3> Vt{0}; + + SVD::SVD(A, U, sigma, Vt); + + // Reference: [16.848103352614, 1.068369514555, 0.0] + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(16.8481f, 1e-2f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1.06837f, 1e-2f)); + // Third singular value should be ~0 (rank-deficient) + REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-2f)); + + // Check reconstruction + Matrix<3, 3> recon{0}; + Matrix<3, 3> Usig{0}; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + Usig[i][j] = U.Get(i, j) * sigma.Get(j, 0); + Usig.Mult(Vt, recon); + + float err = 0.0f; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { + float diff = recon.Get(i, j) - A.Get(i, j); + err += diff * diff; + } + err = sqrtf(err); + REQUIRE_THAT(err, Catch::Matchers::WithinRel(0.0f, 1e-2f)); + + std::cout << "SVD 3x3 rank-deficient:\n"; + std::cout << "Sigma: [" << sigma.Get(0, 0) << ", " << sigma.Get(1, 0) << ", " + << sigma.Get(2, 0) << "]\n"; +} + +TEST_CASE("SVD Integration: tall 4x3 matrix", "[Matrix][SVD][Integration]") { + Matrix<4, 3> A{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}; + Matrix<4, 3> U{0}; + Matrix<3, 1> sigma{0}; + Matrix<3, 3> Vt{0}; + + SVD::SVD(A, U, sigma, Vt); + + // Reference: [25.462407436036, 1.290661675761, 0.0] + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(25.4624f, 1e-2f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1.29066f, 1e-2f)); + REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-2f)); + + // Check reconstruction + Matrix<4, 3> recon{0}; + Matrix<4, 3> Usig{0}; + for (int i = 0; i < 4; i++) + for (int j = 0; j < 3; j++) + Usig[i][j] = U.Get(i, j) * sigma.Get(j, 0); + Usig.Mult(Vt, recon); + + float err = 0.0f; + for (int i = 0; i < 4; i++) + for (int j = 0; j < 3; j++) { + float diff = recon.Get(i, j) - A.Get(i, j); + err += diff * diff; + } + err = sqrtf(err); + REQUIRE_THAT(err, Catch::Matchers::WithinRel(0.0f, 1e-2f)); + + std::cout << "SVD tall 4x3:\n"; + std::cout << "Sigma: [" << sigma.Get(0, 0) << ", " << sigma.Get(1, 0) << ", " + << sigma.Get(2, 0) << "]\n"; +} + +TEST_CASE("SVD Integration: wide 3x5 matrix", "[Matrix][SVD][Integration]") { + Matrix<3, 5> A{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}; + Matrix<3, 5> U{0}; + Matrix<5, 1> sigma{0}; // sigma is columns x 1 = 5x1 for wide matrix + Matrix<5, 5> Vt{0}; // Vt is columns x columns = 5x5 + + SVD::SVD(A, U, sigma, Vt); + + // Reference: [35.127223333575, 2.465396696917, 0.0] + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(35.1272f, 1e-2f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(2.46540f, 1e-2f)); + REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-2f)); + + // Check reconstruction: U (3x5) * diag(sigma) (5x3) = 3x3, then * Vt (3x5) = + // 3x5 + Matrix<3, 5> recon{0}; + Matrix<3, 5> Usig{0}; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 5; j++) + Usig[i][j] = U.Get(i, j) * sigma.Get(j, 0); + + // For wide matrix: A = U * Sigma * Vt where U is 3x5, Sigma is 5x5 + // (diagonal), Vt is 5x5 But our implementation returns sigma as 3x1 and Vt as + // 3x5 So we need: recon = Usig (3x5) * Vt (3x5)^T ... no that doesn't work + // either The SVD for wide matrices is: A = U * Sigma * Vt where: + // U is m×m (3×3), Sigma is m×n (3×5), Vt is n×n (5×5) + // But our API returns U as m×n (3×5), sigma as n×1 (3×1), Vt as n×n (3×5) + // So: recon = U (3x5) * diag(sigma) (5x5) * Vt (5x5)^T ... + // Actually, looking at the implementation, for wide matrices we swap roles. + // Let me just check reconstruction using the actual dimensions returned. + + // For wide matrix: A (3x5) = U (3x5) * diag(sigma) (5x5 padded) * Vt (5x5) + // But our API returns Vt as 3x5, not 5x5 + // The implementation stores: U = QR[:,0:p]^T (3x5), sigma (3x1), Vt = + // QL[:,0:p]^T (3x5) Reconstruction: A[i][j] = sum_k U[i][k]*sigma[k]*Vt[j][k] + + float err2 = 0.0f; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 5; j++) { + float recon_val = 0.0f; + for (int k = 0; k < 3; k++) { + recon_val += U.Get(i, k) * sigma.Get(k, 0) * Vt.Get(j, k); + } + float diff = recon_val - A.Get(i, j); + err2 += diff * diff; + } + } + err2 = sqrtf(err2); + REQUIRE_THAT(err2, Catch::Matchers::WithinRel(0.0f, 1e-2f)); + + std::cout << "SVD wide 3x5:\n"; + std::cout << "Sigma: [" << sigma.Get(0, 0) << ", " << sigma.Get(1, 0) << ", " + << sigma.Get(2, 0) << "]\n"; +} + +TEST_CASE("SVD Integration: identity 3x3", "[Matrix][SVD][Integration]") { + Matrix<3, 3> A{1, 0, 0, 0, 1, 0, 0, 0, 1}; + Matrix<3, 3> U{0}; + Matrix<3, 1> sigma{0}; + Matrix<3, 3> Vt{0}; + + SVD::SVD(A, U, sigma, Vt); + + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-3f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1.0f, 1e-3f)); + REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(1.0f, 1e-3f)); + + float err = frobeniusNorm(U - Matrix<3, 3>{1, 0, 0, 0, 1, 0, 0, 0, 1}); + REQUIRE_THAT(err, Catch::Matchers::WithinRel(0.0f, 1e-2f)); +} + +TEST_CASE("SVD Integration: symmetric positive definite 2x2 [[5,3],[3,5]]", + "[Matrix][SVD][Integration]") { + Matrix<2, 2> A{5, 3, 3, 5}; + Matrix<2, 2> U{0}; + Matrix<2, 1> sigma{0}; + Matrix<2, 2> Vt{0}; + + SVD::SVD(A, U, sigma, Vt); + + // For SPD matrix, singular values = eigenvalues: [8, 2] + REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(8.0f, 1e-3f)); + REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(2.0f, 1e-3f)); + + // Check reconstruction + Matrix<2, 2> recon{0}; + Matrix<2, 2> Usig{0}; + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + Usig[i][j] = U.Get(i, j) * sigma.Get(j, 0); + Usig.Mult(Vt, recon); + + float err = 0.0f; + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) { + float diff = recon.Get(i, j) - A.Get(i, j); + err += diff * diff; + } + err = sqrtf(err); + REQUIRE_THAT(err, Catch::Matchers::WithinRel(0.0f, 1e-3f)); + + std::cout << "SVD SPD 2x2 [[5,3],[3,5]]:\n"; + std::cout << "Sigma: [" << sigma.Get(0, 0) << ", " << sigma.Get(1, 0) + << "]\n"; +} diff --git a/unit-tests/svd-reference-values.py b/unit-tests/svd-reference-values.py new file mode 100644 index 0000000..c274210 --- /dev/null +++ b/unit-tests/svd-reference-values.py @@ -0,0 +1,482 @@ +#!/usr/bin/env python3 +""" +Generate reference values for SVD building block unit tests. +Run this to verify/implement the C++ SVD implementation against scipy/numpy. + +Usage: python3 svd-reference-values.py +""" + +import numpy as np +from scipy.linalg import svd, qr as scipy_qr +import json + +def compute_householder(x): + """Compute Householder reflector: H*x = [alpha, 0, 0, ...]^T. + + Returns (v_normalized, alpha) where v is the normalized Householder vector. + H = I - 2*v*v^T / (v^T*v) + """ + x = np.array(x, dtype=np.float64) + norm_x = np.linalg.norm(x) + + if norm_x < 1e-30: + return x.copy(), 0.0 + + alpha = -np.sign(x[0]) * norm_x if x[0] != 0 else -norm_x + v = x.copy() + v[0] -= alpha + v_norm = np.linalg.norm(v) + + if v_norm < 1e-30: + return np.zeros_like(x), alpha + + v /= v_norm + return v, alpha + +def apply_householder_left(A, v, start_row): + """Apply Householder reflection from the left: A = (I - 2vv^T) @ A. + + v is the normalized Householder vector operating on rows [start_row:]. + The length of v must match the number of rows affected. + """ + A = A.copy() + k = len(v) + + for col in range(A.shape[1]): + dot = np.dot(v, A[start_row:start_row+k, col]) + A[start_row:start_row+k, col] -= 2.0 * dot * v + + return A + +def apply_householder_right(A, v, start_col): + """Apply Householder reflection from the right: A = A @ (I - 2vv^T). + + v is the normalized Householder vector operating on columns [start_col:]. + The length of v must match the number of columns affected. + """ + A = A.copy() + k = len(v) + + for row in range(A.shape[0]): + dot = np.dot(A[row, start_col:start_col+k], v) + A[row, start_col:start_col+k] -= 2.0 * dot * v + + return A + +def compute_givens(x, y): + """Compute Givens rotation that zeros out y. + + Returns (c, s) such that [c s; -s c] @ [x; y] = [r; 0]. + """ + r = np.sqrt(x*x + y*y) + if r < 1e-30: + return 1.0, 0.0 + c = x / r + s = y / r + return c, s + +def apply_givens_left(A, i, j, c, s): + """Apply Givens rotation from the left to rows i and j of A. + + [c s] [row_i] + [-s c] @ [row_j] = [new_row_i] + [new_row_j] + """ + A = A.copy() + new_i = c * A[i] + s * A[j] + new_j = -s * A[i] + c * A[j] + A[i] = new_i + A[j] = new_j + return A + +def apply_givens_right(A, i, j, c, s): + """Apply Givens rotation from the right to columns i and j of A. + + [col_i col_j] @ [c -s] = [new_col_i new_col_j] + [s c] + """ + A = A.copy() + new_i = c * A[:, i] + s * A[:, j] + new_j = -s * A[:, i] + c * A[:, j] + A[:, i] = new_i + A[:, j] = new_j + return A + +def householder_bidiagonalization(A): + """Full Householder bidiagonalization: A = Q_L @ B @ Q_R^T. + + Returns (B, Q_L, Q_R) where B is upper bidiagonal. + """ + m, n = A.shape + p = min(m, n) + + QL = np.eye(m, dtype=np.float64) + QR = np.eye(n, dtype=np.float64) + W = A.copy() + + for k in range(p): + # Left HH: zero out W[k+1:, k] + if k < m - 1: + x = W[k+1:, k].copy() + v, alpha = compute_householder(x) + if np.linalg.norm(v) > 1e-30: + W = apply_householder_left(W, v, k + 1) + QL = apply_householder_right(QL, v, k + 1) + + # Right HH: zero out W[k, k+2:] (superdiagonal) + if k < p - 1 and k + 2 <= n: + x = W[k, k+2:].copy() + v, alpha = compute_householder(x) + if np.linalg.norm(v) > 1e-30: + W = apply_householder_right(W, v, k + 2) + QR = apply_householder_right(QR, v, k + 2) + + return W, QL, QR + +def implicit_qr_iteration(B, QR_acc): + """Implicit QR iteration on a bidiagonal matrix. + + Returns (Sigma, QR_acc) where Sigma is diagonal with singular values + and QR_acc contains the accumulated right transformations. + """ + m, n = B.shape + p = min(m, n) + W = B.copy() + + max_iter = 1000 + tol = 1e-10 + + for iteration in range(max_iter): + # Deflate negligible subdiagonal elements + for i in range(p - 1, 0, -1): + if abs(W[i, i-1]) < tol * (abs(W[i-1, i-1]) + abs(W[i, i])): + W[i, i-1] = 0.0 + + # Find smallest unreduced block [start, end] + start = 0 + for i in range(p - 1): + if abs(W[i+1, i]) >= tol * (abs(W[i, i]) + abs(W[i+1, i+1])): + start = i + 1 + + end = p - 1 + for i in range(p - 2, -1, -1): + if abs(W[i+1, i]) >= tol * (abs(W[i, i]) + abs(W[i+1, i+1])): + end = i + break + + if start >= end: + continue + + # Wilkinson shift from bottom 2x2 corner + a, b = W[end-1, end-1], W[end-1, end] + c_val, d = W[end, end-1], W[end, end] + trace = a + d + det = a * d - b * c_val + disc = trace**2 - 4 * det + + if disc >= 0: + sqrt_disc = np.sqrt(disc) + e1, e2 = (trace + sqrt_disc) / 2, (trace - sqrt_disc) / 2 + shift = e1 if abs(e1 - d) < abs(e2 - d) else e2 + else: + shift = d + + # Implicit QR step using Givens rotations + # Process from top to bottom within the block + x = W[start, start] - shift + y = W[start + 1, start] + + for i in range(start, end): + r = np.sqrt(x*x + y*y) + if r < 1e-30: + x = W[i + 1, i] + y = W[i + 1, i + 1] if i + 2 <= end else 0.0 + continue + + c_rot = x / r + s_rot = y / r + + # Apply from left to rows i, i+1 (columns i..n-1) + for j in range(i, n): + t1, t2 = W[i, j], W[i + 1, j] + W[i, j] = c_rot * t1 + s_rot * t2 + W[i + 1, j] = -s_rot * t1 + c_rot * t2 + + # Apply from right to columns i, i+1 (rows 0..i) + if i > start: + for j in range(i + 1): + t1, t2 = W[j, i], W[j, i + 1] + W[j, i] = c_rot * t1 + s_rot * t2 + W[j, i + 1] = -s_rot * t1 + c_rot * t2 + + # Accumulate into QR_acc + for j in range(QR_acc.shape[0]): + t1, t2 = QR_acc[j, i], QR_acc[j, i + 1] + QR_acc[j, i] = c_rot * t1 + s_rot * t2 + QR_acc[j, i + 1] = -s_rot * t1 + c_rot * t2 + + # Prepare for next rotation + x = W[i + 1, i] + y = W[i + 1, i + 1] if i + 2 <= end else 0.0 + + return W, QR_acc + + +def main(): + print("=" * 70) + print("SVB BUILDING BLOCK REFERENCE VALUES") + print("Generated with scipy/numpy for C++ unit test verification") + print("=" * 70) + + # ------------------------------------------------------------------ + # Test 1: Householder Vector Computation + # ------------------------------------------------------------------ + print("\n" + "=" * 70) + print("TEST 1: computeHouseholderVector") + print("=" * 70) + + test_vectors = [ + ("2D [1,3]", [1.0, 3.0]), + ("2D [3,4] (norm=5)", [3.0, 4.0]), + ("3D [1,2,3]", [1.0, 2.0, 3.0]), + ("3D [0,0,1]", [0.0, 0.0, 1.0]), + ("4D [5,-3,2,1]", [5.0, -3.0, 2.0, 1.0]), + ] + + for name, vec in test_vectors: + v, alpha = compute_householder(vec) + x = np.array(vec) + Hx = x - 2 * np.dot(v, x) * v + + print(f"\n{name}:") + print(f" Input: {list(x)}") + print(f" ||x||: {np.linalg.norm(x):.15f}") + print(f" alpha: {alpha:.15f}") + print(f" v (normalized): {[round(float(vi), 12) for vi in v]}") + print(f" H*x = [alpha,0..]: {[round(float(xi), 12) for xi in Hx]}") + print(f" Off-diagonal ~0: {np.allclose(Hx[1:], 0, atol=1e-12)}") + + # ------------------------------------------------------------------ + # Test 2: Householder Apply Left + # ------------------------------------------------------------------ + print("\n" + "=" * 70) + print("TEST 2: applyHouseholderLeft") + print("=" * 70) + + A_test = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], dtype=np.float64) + x_col = A_test[1:, 0].copy() + v_left, _ = compute_householder(x_col) + + print(f"\nInput matrix:\n{A_test}") + print(f"Householder vector (rows 1:3): {[round(float(vi), 12) for vi in v_left]}") + + A_result = apply_householder_left(A_test, v_left, 1) + print(f"\nAfter applyHouseholderLeft:\n{A_result}") + print(f" A[1,0] = {A_result[1,0]:.2e}, A[2,0] = {A_result[2,0]:.2e} (should be ~0)") + + # ------------------------------------------------------------------ + # Test 3: Householder Apply Right + # ------------------------------------------------------------------ + print("\n" + "=" * 70) + print("TEST 3: applyHouseholderRight") + print("=" * 70) + + A_test = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], dtype=np.float64) + x_row = A_test[0, 1:].copy() + v_right, _ = compute_householder(x_row) + + print(f"\nInput matrix:\n{A_test}") + print(f"Householder vector (cols 1:3): {[round(float(vi), 12) for vi in v_right]}") + + A_result = apply_householder_right(A_test, v_right, 1) + print(f"\nAfter applyHouseholderRight:\n{A_result}") + print(f" A[0,1] = {A_result[0,1]:.2e}, A[0,2] = {A_result[0,2]:.2e} (should be ~0)") + + # ------------------------------------------------------------------ + # Test 4: Givens Rotation Computation + # ------------------------------------------------------------------ + print("\n" + "=" * 70) + print("TEST 4: computeGivens") + print("=" * 70) + + givens_tests = [ + ("3-4-5 triangle", 3.0, 4.0), + ("y already zero", 1.0, 0.0), + ("x is zero", 0.0, 5.0), + ("Both negative", -3.0, -4.0), + ("45 degree case", 1.0, -1.0), + ] + + for name, x, y in givens_tests: + c, s = compute_givens(x, y) + result_x = c * x + s * y + result_y = -s * x + c * y + + print(f"\n{name}: x={x}, y={y}") + print(f" r = {np.sqrt(x*x+y*y):.12f}") + print(f" c = {c:.12f}, s = {s:.12f}") + print(f" [c s; -s c] @ [x;y] = [{result_x:.2e}, {result_y:.2e}]") + + # ------------------------------------------------------------------ + # Test 5: Apply Givens Left/Right + # ------------------------------------------------------------------ + print("\n" + "=" * 70) + print("TEST 5: applyGivensLeft / applyGivensRight") + print("=" * 70) + + A_test = np.array([[3.0, 4.0], [1.0, 2.0]], dtype=np.float64) + c, s = compute_givens(3.0, 1.0) + + print(f"\nInput matrix:\n{A_test}") + print(f"Givens rotation (rows 0,1): c={c:.12f}, s={s:.12f}") + + A_left = apply_givens_left(A_test, 0, 1, c, s) + print(f"\nAfter applyGivensLeft:\n{A_left}") + print(f" A[1,0] = {A_left[1,0]:.2e} (should be ~0)") + + A_test = np.array([[3.0, 1.0], [4.0, 2.0]], dtype=np.float64) + c, s = compute_givens(3.0, 4.0) + + print(f"\nInput matrix:\n{A_test}") + print(f"Givens rotation (cols 0,1): c={c:.12f}, s={s:.12f}") + + A_right = apply_givens_right(A_test, 0, 1, c, s) + print(f"\nAfter applyGivensRight:\n{A_right}") + print(f" A[0,1] = {A_right[0,1]:.2e} (should be ~0)") + + # ------------------------------------------------------------------ + # Test 6: Full Bidiagonalization + # ------------------------------------------------------------------ + print("\n" + "=" * 70) + print("TEST 6: householderBidiagonalization") + print("=" * 70) + + bidiag_tests = [ + ("2x2 [[1,2],[3,4]]", np.array([[1.0, 2.0], [3.0, 4.0]])), + ("3x3 SPD [[5,3],[3,5]]", np.array([[5.0, 3.0], [3.0, 5.0]])), + ("3x3 diag [[10,0,0],[0,5,0],[0,0,2]]", + np.array([[10.0, 0, 0], [0, 5.0, 0], [0, 0, 2.0]])), + ("3x3 full [[1,2,3],[4,5,6],[7,8,10]]", + np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 10.0]])), + ("Tall 4x3", np.array([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], dtype=np.float64)), + ] + + for name, A in bidiag_tests: + B, QL, QR = householder_bidiagonalization(A) + m, n = A.shape + p = min(m, n) + + print(f"\n{name}:") + print(f" Original:\n{A}") + print(f"\n Bidiagonal B:\n{B}") + print(f" Diagonal: {[round(float(B[i,i]), 10) for i in range(p)]}") + print(f" Superdiag: {[round(float(B[i,i+1]), 10) for i in range(min(p-1, n-1))]}") + + recon = QL @ B @ QR.T + err = np.linalg.norm(recon - A, 'fro') + print(f" ||QL @ B @ QR^T - A||_F = {err:.2e}") + + # ------------------------------------------------------------------ + # Test 7: Full SVD Reference Values + # ------------------------------------------------------------------ + print("\n" + "=" * 70) + print("TEST 7: Full SVD Reference Values (scipy.linalg.svd)") + print("=" * 70) + + test_matrices = [ + ("Simple 2x2", np.array([[1,2],[3,4]], dtype=np.float64)), + ("SPD 2x2", np.array([[5,3],[3,5]], dtype=np.float64)), + ("Full-rank 3x3", np.array([[1,2,3],[4,5,6],[7,8,10]], dtype=np.float64)), + ("Rank-deficient 3x3", np.array([[1,2,3],[4,5,6],[7,8,9]], dtype=np.float64)), + ("Diagonal 3x3", np.array([[10,0,0],[0,5,0],[0,0,2]], dtype=np.float64)), + ("Tall 4x3", np.array([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], dtype=np.float64)), + ("Wide 3x5", np.array([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15]], dtype=np.float64)), + ("Symmetric tri 5x5", np.array([[2,-1,0,0,0],[-1,2,-1,0,0],[0,-1,2,-1,0],[0,0,-1,2,-1],[0,0,0,-1,2]], dtype=np.float64)), + ("Neg values 2x3", np.array([[0.5,-0.3,0.8],[-0.2,0.7,0.1]], dtype=np.float64)), + ("Near-singular 2x2", np.array([[1,0],[0,1e-6]], dtype=np.float64)), + ("Orthogonal 3x3", np.array([[np.cos(np.pi/4), -np.sin(np.pi/4), 0], + [np.sin(np.pi/4), np.cos(np.pi/4), 0], + [0, 0, 1]], dtype=np.float64)), + ("Identity 3x3", np.eye(3)), + ("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)), + ] + + for name, A in test_matrices: + U, s, Vt = svd(A, full_matrices=False) + + print(f"\n{name}: shape={A.shape}") + print(f" Singular values: {[round(float(x), 12) for x in s]}") + print(f" U:\n{np.array2string(U, precision=6, floatmode='maxprec_equal')}") + print(f" Vt:\n{np.array2string(Vt, precision=6, floatmode='maxprec_equal')}") + recon_err = np.linalg.norm(A - U @ np.diag(s) @ Vt, 'fro') + print(f" Reconstruction error: {recon_err:.2e}") + + # ------------------------------------------------------------------ + # Test 8: Implicit QR Iteration on Bidiagonal + # ------------------------------------------------------------------ + print("\n" + "=" * 70) + print("TEST 8: implicitQRIteration") + print("=" * 70) + + qr_tests = [ + ("2x2 [[1,2],[3,4]]", np.array([[1.0, 2.0], [3.0, 4.0]])), + ("3x3 diag", np.array([[10.0, 0, 0], [0, 5.0, 0], [0, 0, 2.0]])), + ] + + for name, A in qr_tests: + B, QL, QR = householder_bidiagonalization(A) + Sigma, QR_final = implicit_qr_iteration(B.copy(), QR.copy()) + + print(f"\n{name}:") + print(f" Bidiagonal B:\n{B}") + print(f" After QR iteration (Sigma):\n{Sigma}") + print(f" Diagonal entries: {[round(float(Sigma[i,i]), 10) for i in range(min(Sigma.shape))]}") + + # Verify: QL @ Sigma @ QR_final^T ≈ A + recon = QL @ Sigma @ QR_final.T + err = np.linalg.norm(recon - A, 'fro') + print(f" ||QL @ Sigma @ QR^T - A||_F = {err:.2e}") + + # ------------------------------------------------------------------ + # JSON output for easy import into C++ tests + # ------------------------------------------------------------------ + print("\n" + "=" * 70) + print("JSON OUTPUT (for easy C++ integration)") + print("=" * 70) + + json_data = {} + + # Householder test vectors + hh_tests = {} + for name, vec in test_vectors: + v, alpha = compute_householder(vec) + x = np.array(vec) + Hx = x - 2 * np.dot(v, x) * v + hh_tests[name] = { + "input": [float(xi) for xi in x], + "norm": float(np.linalg.norm(x)), + "alpha": float(alpha), + "v_normalized": [round(float(vi), 12) for vi in v], + "Hx": [round(float(xi), 12) for xi in Hx], + } + json_data["householder_vectors"] = hh_tests + + # Full SVD reference values + svd_tests = {} + for name, A in test_matrices: + U, s, Vt = svd(A, full_matrices=False) + svd_tests[name] = { + "shape": list(A.shape), + "singular_values": [round(float(x), 12) for x in s], + "U": [[round(float(U[i,j]), 8) for j in range(U.shape[1])] for i in range(U.shape[0])], + "Vt": [[round(float(Vt[i,j]), 8) for j in range(Vt.shape[1])] for i in range(Vt.shape[0])], + } + json_data["svd_reference"] = svd_tests + + print(json.dumps(json_data, indent=2)) + + +if __name__ == "__main__": + main()