Working on an SVD implimentation

CMakeLists.txt

@@ -4,7 +4,7 @@ project(Vector3D)
 add_subdirectory(src)
 add_subdirectory(unit-tests)
 
-set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD 17)
 
 add_compile_options(-Wall -Wextra -Wpedantic)
 add_compile_options (-fdiagnostics-color=always)

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
 )

src/Matrix.cpp

@@ -20,7 +20,8 @@ Matrix<rows, columns>::Matrix(const std::array<float, rows * columns> &array) {
 }
 
 template <uint8_t rows, uint8_t columns>
-template <typename... Args>
+template <typename... Args,
+          std::enable_if_t<(std::is_arithmetic_v<Args> && ...), int>>
 Matrix<rows, columns>::Matrix(Args... args) {
   constexpr uint16_t arraySize{static_cast<uint16_t>(rows) *
                                static_cast<uint16_t>(columns)};
@@ -531,7 +532,7 @@ void Matrix<rows, columns>::QRDecomposition(Matrix<rows, columns> &Q,
   Q.Fill(0);
   R.Fill(0);
   Matrix<rows, 1> a_col, e, u, Q_column_k{};
-  Matrix<1, rows> a_T, e_T{};
+  Matrix<1, rows> e_T{};
 
   for (uint8_t column = 0; column < columns; column++) {
     this->GetColumn(column, a_col);

src/Matrix.hpp

@@ -3,6 +3,7 @@
 #include <array>
 #include <cstdint>
 #include <string>
+#include <type_traits>
 
 // TODO: Add a function to calculate eigenvalues/vectors
 // TODO: Add a function to compute RREF
@@ -29,9 +30,12 @@ public:
   Matrix(const Matrix<rows, columns> &other);
 
   /**
-   * @brief Initialize a matrix directly with any number of arguments
+   * @brief Initialize a matrix directly with scalar values
+   * Uses SFINAE to only accept arithmetic types (int, float, double, etc.)
    */
-  template <typename... Args> Matrix(Args... args);
+  template <typename... Args,
+            std::enable_if_t<(std::is_arithmetic_v<Args> && ...), int> = 0>
+  Matrix(Args... args);
 
   /**
    * @brief Create an identity matrix

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 <cstdint>
+
+// ============================================================================
+// 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<m, n> — first k columns are meaningful
+ *   - sigma: Matrix<n, 1> — first k entries are non-zero singular values
+ *   - Vt: Matrix<n, n> — first k rows are meaningful
+ *
+ * For m < n (wide matrices), we work with Aᵀ and swap roles of U and V.
+ */
+template <uint8_t rows, uint8_t columns>
+void SVD::SVD(Matrix<rows, columns> &matrixToDecompose,
+              Matrix<rows, columns> &U, Matrix<columns, 1> &sigma,
+              Matrix<columns, columns> &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

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 <uint8_t rows, uint8_t columns>
+void SVD(Matrix<rows, columns> &matrixToDecompose, Matrix<rows, columns> &U,
+         Matrix<columns, 1> &sigma, Matrix<columns, columns> &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_

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
 )

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 <array>
@@ -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 <uint8_t rows, uint8_t columns>
+static float frobeniusNorm(const Matrix<rows, columns> &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 <uint8_t rows, uint8_t columns>
+static float svdReconstructionError(const Matrix<rows, columns> &A,
+                                    const Matrix<rows, columns> &U,
+                                    const Matrix<columns, 1> &sigma,
+                                    const Matrix<columns, columns> &Vt) {
+  // Compute U × diag(σ): only first min(rows,cols) columns of U are used
+  constexpr uint8_t k = (rows < columns) ? rows : columns;
+  Matrix<rows, columns> 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<rows, columns> 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<rows, columns> 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 <uint8_t rows, uint8_t columns>
+static float orthogonalityError(const Matrix<rows, columns> &M) {
+  constexpr uint8_t k = (rows < columns) ? rows : columns;
+
+  // Compute Mᵀ × M (should be I_k in top-left)
+  Matrix<columns, rows> Mt = M.Transpose();
+  Matrix<columns, columns> 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 <uint8_t maxCols>
+static bool isSortedDescending(const Matrix<maxCols, 1> &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));
 }

unit-tests/matrix-timing-tests.cpp

@@ -76,7 +76,8 @@ TEST_CASE("Timing Tests", "Matrix") {
 
   SECTION("Determinant") {
     for (uint32_t i{0}; i < 1000000; i++) {
-      float det1 = mat4.Det();
+      float det = mat4.Det();
+      (void)det;
     }
   }
 

unit-tests/svd-build-blocks-tests.cpp

@@ -0,0 +1,785 @@
+// include the unit test framework first
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+// include the module you're going to test next
+#include "Matrix.hpp"
+#include "SVD.hpp"
+
+// any other libraries
+#include <array>
+#include <cmath>
+#include <iostream>
+
+// ============================================================================
+// 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));
+  }
+}

unit-tests/svd-integration-test.cpp

@@ -0,0 +1,269 @@
+#include "Matrix.hpp"
+#include "SVD.hpp"
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <iostream>
+
+// Generic helper functions for any matrix size
+template <uint8_t rows, uint8_t columns>
+static float frobeniusNorm(const Matrix<rows, columns> &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 <uint8_t n>
+static bool isOrthogonal(const Matrix<n, n> &M, float tol = 1e-4f) {
+  Matrix<n, n> Mt = M.Transpose();
+  Matrix<n, n> 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";
+}

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