Fix SVD pipeline bugs; extract and unit-test bidiagonal block solvers
- Snapshot bidiagonal block before solving (W overwrite corrupted Ublock) - Wide-matrix QL row extent fix (rowsQL = n when transposed) - Vt output copy bound j<columns (was OOB write clobbering caller sigma) - Rank-deficient blocks: residual singular values sigma_i = ||B*v_i|| instead of sqrt(eigs of B*B) (condition-number squaring); Gram-Schmidt orthonormalization of U columns with rank-deficiency completion - JacobiEigenSymmetric: return signed eigenvalues (removed fabsf) - Fix 28 WithinRel(0.0f,.) matcher misuse -> WithinAbs in matrix-tests - New unit tests: SolveBidiagonalBlock2x2 (7 scipy-referenced cases), JacobiEigenSymmetric (3 cases), DeflateBidiagonal/BidiagonalIsDiagonal
This commit is contained in:
+137
-18
@@ -421,7 +421,10 @@ void SVD::JacobiEigenSymmetric(float T[5][5], uint8_t n, float evals[5],
|
||||
}
|
||||
|
||||
for (uint8_t i = 0; i < n; i++) {
|
||||
evals[i] = fabsf(T[i][i]);
|
||||
// NOTE: no fabsf — the eigenvalues keep their sign (this is a
|
||||
// general symmetric eigen solver, not just for PSD matrices like
|
||||
// T = BᵀB). Callers needing magnitudes take them.
|
||||
evals[i] = T[i][i];
|
||||
}
|
||||
}
|
||||
|
||||
@@ -441,43 +444,61 @@ void SVD::ApplyBlockFactorsToAccumulators(uint8_t blockStart, uint8_t blockSize,
|
||||
// rowsQL / rowsQR are the meaningful row extents of the accumulators:
|
||||
// for a transposed (wide) problem W = Aᵀ has n rows, so QL carries n
|
||||
// meaningful rows while in the normal case it carries m.
|
||||
//
|
||||
// BUG FIX: Use temporary buffers to avoid in-place corruption.
|
||||
// The old code updated QL[j][blockStart+i] while still reading from
|
||||
// QL[j][blockStart+k] for later i values, corrupting subsequent columns.
|
||||
|
||||
float newQL[5][5] = {{0}};
|
||||
for (uint8_t j = 0; j < rowsQL; j++) {
|
||||
for (uint8_t i = 0; i < blockSize; i++) {
|
||||
float sum = 0.0f;
|
||||
for (uint8_t k = 0; k < blockSize; k++) {
|
||||
sum += QL[j][blockStart + k] * Ublock[k][i];
|
||||
}
|
||||
QL[j][blockStart + i] = sum;
|
||||
newQL[j][blockStart + i] = sum;
|
||||
}
|
||||
}
|
||||
for (uint8_t j = 0; j < rowsQL; j++)
|
||||
for (uint8_t i = 0; i < blockSize; i++)
|
||||
QL[j][blockStart + i] = newQL[j][blockStart + i];
|
||||
|
||||
float newQR[5][5] = {{0}};
|
||||
for (uint8_t j = 0; j < rowsQR; j++) {
|
||||
for (uint8_t i = 0; i < blockSize; i++) {
|
||||
float sum = 0.0f;
|
||||
for (uint8_t k = 0; k < blockSize; k++) {
|
||||
sum += QR[j][blockStart + k] * Vblock[k][i];
|
||||
}
|
||||
QR[j][blockStart + i] = sum;
|
||||
newQR[j][blockStart + i] = sum;
|
||||
}
|
||||
}
|
||||
for (uint8_t j = 0; j < rowsQR; j++)
|
||||
for (uint8_t i = 0; i < blockSize; i++)
|
||||
QR[j][blockStart + i] = newQR[j][blockStart + i];
|
||||
}
|
||||
|
||||
void SVD::SolveBidiagonalBlockJacobi(Matrix<5, 5> &W, uint8_t blockStart,
|
||||
uint8_t blockSize, uint8_t rowsQL,
|
||||
uint8_t rowsQR, Matrix<5, 5> &QL,
|
||||
Matrix<5, 5> &QR, float tol) {
|
||||
// Full SVD of an unreduced upper-bidiagonal block of size > 2 via
|
||||
// eigen-decomposition of the tridiagonal T = BᵀB:
|
||||
// Full SVD of an unreduced upper-bidiagonal block of size > 2, computed
|
||||
// as the eigen-decomposition of the symmetric tridiagonal T = BᵀB:
|
||||
//
|
||||
// 1. Snapshot the ORIGINAL block diagonal/superdiagonal from W
|
||||
// 2. Form T = BᵀB (tridiagonal symmetric)
|
||||
// 3. JacobiEigenSymmetric → eigenvalues + eigenvector matrix V
|
||||
// 4. Sort eigenvalues descending, reordering V
|
||||
// 5. Ublock = B_orig · V · Σ⁻¹ (computed from the SNAPSHOT so that
|
||||
// overwriting W's diagonal does not corrupt it)
|
||||
// 5. Singular values are the RESIDUAL norms σᵢ = ‖B·vᵢ‖ rather than
|
||||
// sqrt(evals[i]): forming BᵀB squares the condition number, so
|
||||
// float noise in T swamps the smallest eigenvalues of
|
||||
// rank-deficient / near-deficient blocks (σ error ~1e-3 instead of
|
||||
// ~1e-6). ‖B·vᵢ‖ stays accurate to ~eps·‖B‖. It also makes
|
||||
// uᵢ = B·vᵢ/σᵢ unit-norm by construction; when σᵢ ≈ 0 (true rank
|
||||
// deficiency) uᵢ is replaced by a Gram–Schmidt orthogonal
|
||||
// completion so the accumulated QL/QR stay orthogonal.
|
||||
// 6. Fold Ublock/Vblock into QL/QR via ApplyBlockFactorsToAccumulators
|
||||
// 7. Only now write sqrt(eigenvalues) into W's diagonal and zero the
|
||||
// block's superdiagonals
|
||||
// 7. Only now write σ onto W's diagonal and zero the superdiagonals
|
||||
(void)tol; // Jacobi convergence tolerance is internal
|
||||
|
||||
// Step 1: snapshot original block values (diagonal d[i], superdiag e[i])
|
||||
@@ -506,12 +527,10 @@ void SVD::SolveBidiagonalBlockJacobi(Matrix<5, 5> &W, uint8_t blockStart,
|
||||
}
|
||||
}
|
||||
|
||||
// Step 3: Jacobi eigenvalue algorithm
|
||||
// Steps 3-4: Jacobi eigen-decomposition, then sort descending
|
||||
float evals[5] = {0};
|
||||
float V[5][5] = {{0}};
|
||||
SVD::JacobiEigenSymmetric(T, blockSize, evals, V);
|
||||
|
||||
// Step 4: sort eigenvalues descending, reordering eigenvector columns
|
||||
for (uint8_t i = 0; i < blockSize - 1; i++) {
|
||||
for (uint8_t j = i + 1; j < blockSize; j++) {
|
||||
if (evals[j] > evals[i]) {
|
||||
@@ -527,17 +546,117 @@ void SVD::SolveBidiagonalBlockJacobi(Matrix<5, 5> &W, uint8_t blockStart,
|
||||
}
|
||||
}
|
||||
|
||||
// Step 5: Ublock = B_orig · V · Σ⁻¹, from the SNAPSHOT values.
|
||||
// Column i of Ublock is u_i = (B_orig · v_i) / σ_i.
|
||||
float Ublock[5][5] = {{0}};
|
||||
// Step 5: residual singular values σᵢ = ‖B·vᵢ‖ and unit-norm uᵢ
|
||||
float Bv[5][5];
|
||||
for (uint8_t i = 0; i < blockSize; i++) {
|
||||
float sigmaI = sqrtf(evals[i]);
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
float result = d[r] * V[r][i];
|
||||
if (r + 1 < blockSize) {
|
||||
result += e[r] * V[r + 1][i];
|
||||
}
|
||||
Ublock[r][i] = (sigmaI > 1e-30f) ? result / sigmaI : 0.0f;
|
||||
Bv[r][i] = result;
|
||||
}
|
||||
}
|
||||
float sigma[5];
|
||||
for (uint8_t i = 0; i < blockSize; i++) {
|
||||
float n = 0.0f;
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
n += Bv[r][i] * Bv[r][i];
|
||||
}
|
||||
sigma[i] = sqrtf(n);
|
||||
}
|
||||
// Re-sort sigma descending, swapping V AND Bv columns consistently
|
||||
// (residual norms can differ slightly in ordering from sqrt(evals))
|
||||
for (uint8_t i = 0; i < blockSize - 1; i++) {
|
||||
for (uint8_t j = i + 1; j < blockSize; j++) {
|
||||
if (sigma[j] > sigma[i]) {
|
||||
float tmpS = sigma[i];
|
||||
sigma[i] = sigma[j];
|
||||
sigma[j] = tmpS;
|
||||
for (uint8_t k = 0; k < blockSize; k++) {
|
||||
float tmpV = V[k][i];
|
||||
V[k][i] = V[k][j];
|
||||
V[k][j] = tmpV;
|
||||
float tmpB = Bv[k][i];
|
||||
Bv[k][i] = Bv[k][j];
|
||||
Bv[k][j] = tmpB;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
float Ublock[5][5] = {{0}};
|
||||
// Pass 1: raw u-columns. For non-degenerate σ, uᵢ = B·vᵢ/σᵢ.
|
||||
for (uint8_t i = 0; i < blockSize; i++) {
|
||||
if (sigma[i] > 1e-30f) {
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
Ublock[r][i] = Bv[r][i] / sigma[i];
|
||||
}
|
||||
}
|
||||
// Degenerate (sigma[i] ~ 0): leave zero, completed in pass 2
|
||||
}
|
||||
// Pass 2: enforce a full orthonormal set. Even when B·vᵢ ≠ 0, the
|
||||
// smallest-σ columns come from eigenvectors of the noise-dominated
|
||||
// tail of T = BᵀB, so their u-directions are near-random and NOT
|
||||
// mutually orthogonal. Gram–Schmidt against the previous columns and
|
||||
// re-normalize (no-op for the well-conditioned columns); if the
|
||||
// residual is ~0 (true rank deficiency) pick a basis-vector seed and
|
||||
// complete orthonally instead.
|
||||
for (uint8_t i = 0; i < blockSize; i++) {
|
||||
for (uint8_t p = 0; p < i; p++) {
|
||||
float dot = 0.0f;
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
dot += Ublock[r][i] * Ublock[r][p];
|
||||
}
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
Ublock[r][i] -= dot * Ublock[r][p];
|
||||
}
|
||||
}
|
||||
float nn = 0.0f;
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
nn += Ublock[r][i] * Ublock[r][i];
|
||||
}
|
||||
if (nn > 1e-12f) {
|
||||
float inv = 1.0f / sqrtf(nn);
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
Ublock[r][i] *= inv;
|
||||
}
|
||||
} else {
|
||||
// True rank deficiency: direction arbitrary. Seed with a basis
|
||||
// vector, Gram–Schmidt against previous columns, normalize.
|
||||
bool found = false;
|
||||
for (uint8_t seed = 0; seed < blockSize && !found; seed++) {
|
||||
float g[5];
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
g[r] = (r == seed) ? 1.0f : 0.0f;
|
||||
}
|
||||
for (uint8_t p = 0; p < i; p++) {
|
||||
float dot = 0.0f;
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
dot += g[r] * Ublock[r][p];
|
||||
}
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
g[r] -= dot * Ublock[r][p];
|
||||
}
|
||||
}
|
||||
float gg = 0.0f;
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
gg += g[r] * g[r];
|
||||
}
|
||||
if (gg > 1e-12f) {
|
||||
float inv = 1.0f / sqrtf(gg);
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
Ublock[r][i] = g[r] * inv;
|
||||
}
|
||||
found = true;
|
||||
}
|
||||
}
|
||||
if (!found) {
|
||||
// Degenerate fallback (cannot happen for i < blockSize): e_0
|
||||
for (uint8_t r = 0; r < blockSize; r++) {
|
||||
Ublock[r][i] = (r == 0) ? 1.0f : 0.0f;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -548,7 +667,7 @@ void SVD::SolveBidiagonalBlockJacobi(Matrix<5, 5> &W, uint8_t blockStart,
|
||||
// Step 7: W last — write singular values onto the diagonal and zero
|
||||
// the block's superdiagonals
|
||||
for (uint8_t i = 0; i < blockSize; i++) {
|
||||
W[blockStart + i][blockStart + i] = sqrtf(evals[i]);
|
||||
W[blockStart + i][blockStart + i] = sigma[i];
|
||||
if (i < blockSize - 1) {
|
||||
W[blockStart + i][blockStart + i + 1] = 0.0f;
|
||||
}
|
||||
|
||||
+13
-7
@@ -224,13 +224,19 @@ static void ApplyBlockFactorsToAccumulators(uint8_t blockStart,
|
||||
* W[blockStart..blockStart+blockSize−1] via eigendecomposition of the
|
||||
* tridiagonal T = BᵀB:
|
||||
* 1. Snapshot the ORIGINAL block diagonal/superdiagonal from W
|
||||
* 2. JacobiEigenSymmetric on T → eigenvalues (unsorted) + V
|
||||
* 3. Sort eigenvalues descending, reordering V
|
||||
* 4. Ublock = B_orig · V · Σ⁻¹ (from the snapshot, so W is not
|
||||
* overwritten before Ublock is computed)
|
||||
* 5. Fold Ublock/Vblock into QL/QR via ApplyBlockFactorsToAccumulators
|
||||
* 6. Only then write sqrt(eigenvalues) into W's diagonal and zero the
|
||||
* block's superdiagonals
|
||||
* 2. Form T = BᵀB (tridiagonal symmetric)
|
||||
* 3. JacobiEigenSymmetric on T → eigenvalues (unsorted) + V
|
||||
* 4. Sort eigenvalues descending, reordering V columns
|
||||
* 5. Compute RESIDUAL singular values: σᵢ = ‖B_orig · vᵢ‖
|
||||
* (NOT sqrt(eigenvalue) — forming BᵀB squares the condition number,
|
||||
* causing float noise to swamp true tiny eigenvalues for
|
||||
* rank-deficient blocks)
|
||||
* 6. Re-sort σ descending, keeping V consistent
|
||||
* 7. Build Ublock: uᵢ = B_orig · vᵢ / σᵢ (unit norm); for σᵢ ≈ 0,
|
||||
* use Gram-Schmidt orthogonal completion against prior U columns
|
||||
* 8. Fold Ublock/Vblock into QL/QR via ApplyBlockFactorsToAccumulators
|
||||
* 9. Write residual norms into W's diagonal and zero the block's
|
||||
* superdiagonals
|
||||
*
|
||||
* @param W Input/output: bidiagonal matrix (5×5 working array); the block's
|
||||
* diagonal holds the singular values and its superdiagonals are
|
||||
|
||||
@@ -826,7 +826,8 @@ TEST_CASE("SVD: Rank-Deficient 3x3 Matrix", "Matrix") {
|
||||
|
||||
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};
|
||||
// Row-major init: [10,0,0, 0,5,0, 0,0,2] = diag(10,5,2)
|
||||
Matrix<3, 3> A{10.0f, 0.0f, 0.0f, 0.0f, 5.0f, 0.0f, 0.0f, 0.0f, 2.0f};
|
||||
Matrix<3, 3> U{}, Vt{};
|
||||
Matrix<3, 1> sigma{};
|
||||
|
||||
|
||||
@@ -1406,3 +1406,232 @@ TEST_CASE("SVD Phase 1: Bidiagonalize reconstruction property", "[Matrix][SVD]")
|
||||
REQUIRE_THAT(err, Catch::Matchers::WithinAbs(1e-3f, 1e-3f));
|
||||
}
|
||||
}
|
||||
|
||||
// ============================================================================
|
||||
// TEST: SolveBidiagonalBlock2x2 — 2×2 upper-bidiagonal block SVD
|
||||
// ============================================================================
|
||||
// Reference singular values generated with scipy.linalg.svd for
|
||||
// B = [[a, b], [0, d]].
|
||||
TEST_CASE("SVD Building Block: SolveBidiagonalBlock2x2", "[Matrix][SVD]") {
|
||||
struct Case2x2 {
|
||||
float a, b, d;
|
||||
float refSigma[2];
|
||||
};
|
||||
const Case2x2 cases[] = {
|
||||
{2.5f, -1.3f, 0.8f, {2.84346151f, 0.70336806f}},
|
||||
{3.0f, 0.0f, 1.0f, {3.0f, 1.0f}},
|
||||
{1.0f, 2.0f, 0.0f, {2.23606798f, 0.0f}},
|
||||
{-1.5f, 0.7f, -2.2f, {2.37779179f, 1.38784229f}},
|
||||
{1.0f, 1e-4f, 0.0f, {1.0f, 0.0f}},
|
||||
{-1.770486f, 0.281880f, 0.208573f, {1.79308863f, 0.20594385f}},
|
||||
{0.866025f, 1.0f, 0.5f, {1.37890797f, 0.31402567f}},
|
||||
};
|
||||
|
||||
for (const auto &tc : cases) {
|
||||
float Ublock[2][2] = {{0}}, Vblock[2][2] = {{0}}, sigma[2] = {0};
|
||||
SVD::SolveBidiagonalBlock2x2(tc.a, tc.b, tc.d, Ublock, Vblock, sigma);
|
||||
|
||||
// 1. Singular values match scipy
|
||||
REQUIRE_THAT(sigma[0],
|
||||
Catch::Matchers::WithinRel(tc.refSigma[0], 1e-3f));
|
||||
if (tc.refSigma[1] > 0.0f) {
|
||||
REQUIRE_THAT(sigma[1],
|
||||
Catch::Matchers::WithinRel(tc.refSigma[1], 1e-3f));
|
||||
} else {
|
||||
REQUIRE(sigma[1] < 1e-3f);
|
||||
}
|
||||
REQUIRE(sigma[0] >= sigma[1]);
|
||||
|
||||
// 2. Ublock and Vblock are orthogonal (MᵀM = I)
|
||||
for (int i = 0; i < 2; i++) {
|
||||
for (int j = i; j < 2; j++) {
|
||||
float dotU = Ublock[0][i] * Ublock[0][j] + Ublock[1][i] * Ublock[1][j];
|
||||
float dotV = Vblock[0][i] * Vblock[0][j] + Vblock[1][i] * Vblock[1][j];
|
||||
float expected = (i == j) ? 1.0f : 0.0f;
|
||||
REQUIRE_THAT(dotU, Catch::Matchers::WithinAbs(expected, 1e-3f));
|
||||
REQUIRE_THAT(dotV, Catch::Matchers::WithinAbs(expected, 1e-3f));
|
||||
}
|
||||
}
|
||||
|
||||
// 3. Ublock · diag(sigma) · Vblockᵀ reproduces B = [[a,b],[0,d]]
|
||||
// (C[i][j] = sum_k U[i][k] * sigma[k] * V[j][k])
|
||||
float C[2][2] = {{0}, {0}};
|
||||
for (int i = 0; i < 2; i++)
|
||||
for (int j = 0; j < 2; j++)
|
||||
for (int k = 0; k < 2; k++)
|
||||
C[i][j] += Ublock[i][k] * sigma[k] * Vblock[j][k];
|
||||
REQUIRE_THAT(C[0][0], Catch::Matchers::WithinAbs(tc.a, 1e-2f));
|
||||
REQUIRE_THAT(C[0][1], Catch::Matchers::WithinAbs(tc.b, 1e-2f));
|
||||
REQUIRE_THAT(C[1][0], Catch::Matchers::WithinAbs(0.0f, 1e-2f));
|
||||
REQUIRE_THAT(C[1][1], Catch::Matchers::WithinAbs(tc.d, 1e-2f));
|
||||
}
|
||||
}
|
||||
|
||||
// ============================================================================
|
||||
// TEST: JacobiEigenSymmetric — cyclic Jacobi eigenvalue decomposition
|
||||
// ============================================================================
|
||||
// Reference eigenvalues generated with scipy.linalg.eigvalsh (desc).
|
||||
TEST_CASE("SVD Building Block: JacobiEigenSymmetric", "[Matrix][SVD]") {
|
||||
struct CaseJac {
|
||||
float S[5][5];
|
||||
uint8_t n;
|
||||
float refEig[5];
|
||||
};
|
||||
|
||||
// (i) T = BᵀB from a real bidiagonalization (3×3)
|
||||
float T3[5][5] = {
|
||||
{65.999993f, -124.470864f, 0.0f, 0.0f, 0.0f},
|
||||
{-124.470864f, 237.877008f, -0.499065f, 0.0f, 0.0f},
|
||||
{0.0f, -0.499065f, 0.122959f, 0.0f, 0.0f},
|
||||
{0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
|
||||
{0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
|
||||
};
|
||||
// (ii) random-looking 3×3 symmetric (seed 42)
|
||||
float S3[5][5] = {
|
||||
{0.304717f, -0.04971f, 0.439146f, 0.0f, 0.0f},
|
||||
{-0.04971f, -1.951035f, -0.809211f, 0.0f, 0.0f},
|
||||
{0.439146f, -0.809211f, -0.016801f, 0.0f, 0.0f},
|
||||
{0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
|
||||
{0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
|
||||
};
|
||||
// (iii) random-looking 4×4 symmetric (seed 42)
|
||||
float S4[5][5] = {
|
||||
{-0.853044f, 1.00332f, -0.090545f, -0.307449f, 0.0f},
|
||||
{1.00332f, 0.467509f, 0.009579f, 0.795646f, 0.0f},
|
||||
{-0.090545f, 0.009579f, -0.049926f, -0.169696f, 0.0f},
|
||||
{-0.307449f, 0.795646f, -0.169696f, -0.428328f, 0.0f},
|
||||
{0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
|
||||
};
|
||||
|
||||
float refs[3][5] = {
|
||||
{303.195295f, 0.765908223f, 0.0387564408f, 0, 0},
|
||||
{0.7227162f, -0.13661881f, -2.24921639f, 0, 0},
|
||||
{1.22127596f, -0.01555681f, -0.31307273f, -1.75643542f, 0},
|
||||
};
|
||||
uint8_t ns[3] = {3, 3, 4};
|
||||
float (*mats[3])[5] = {T3, S3, S4};
|
||||
float maxAbs[3] = {237.877008f, 1.951035f, 1.00332f};
|
||||
|
||||
for (int c = 0; c < 3; c++) {
|
||||
float T[5][5];
|
||||
for (int i = 0; i < 5; i++)
|
||||
for (int j = 0; j < 5; j++)
|
||||
T[i][j] = mats[c][i][j];
|
||||
float S_orig[5][5];
|
||||
for (int i = 0; i < 5; i++)
|
||||
for (int j = 0; j < 5; j++)
|
||||
S_orig[i][j] = mats[c][i][j];
|
||||
|
||||
float evals[5] = {0};
|
||||
float V[5][5] = {{0}};
|
||||
SVD::JacobiEigenSymmetric(T, ns[c], evals, V);
|
||||
|
||||
// 1. Sorted eigenvalues match scipy
|
||||
float sorted[5] = {0};
|
||||
for (int i = 0; i < ns[c]; i++) sorted[i] = evals[i];
|
||||
// Sort descending to match the scipy reference order
|
||||
for (int i = 0; i < ns[c] - 1; i++) {
|
||||
int maxIdx = i;
|
||||
for (int j = i + 1; j < ns[c]; j++)
|
||||
if (sorted[j] > sorted[maxIdx])
|
||||
maxIdx = j;
|
||||
if (maxIdx != i) {
|
||||
float t = sorted[i];
|
||||
sorted[i] = sorted[maxIdx];
|
||||
sorted[maxIdx] = t;
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < ns[c]; i++) {
|
||||
if (fabsf(refs[c][i]) > 0.01f) {
|
||||
REQUIRE_THAT(sorted[i],
|
||||
Catch::Matchers::WithinRel(refs[c][i], 1e-3f));
|
||||
} else {
|
||||
REQUIRE_THAT(sorted[i], Catch::Matchers::WithinAbs(refs[c][i], 1e-3f));
|
||||
}
|
||||
}
|
||||
|
||||
// 2. V is orthogonal (VᵀV = I on the n×n part)
|
||||
for (int i = 0; i < ns[c]; i++) {
|
||||
for (int j = i; j < ns[c]; j++) {
|
||||
float dot = 0.0f;
|
||||
for (int k = 0; k < ns[c]; k++) dot += V[k][i] * V[k][j];
|
||||
float expected = (i == j) ? 1.0f : 0.0f;
|
||||
REQUIRE_THAT(dot, Catch::Matchers::WithinAbs(expected, 1e-3f));
|
||||
}
|
||||
}
|
||||
|
||||
// 3. Residual ‖S_orig·V − V·diag(evals)‖ small
|
||||
// (col i of S_orig·V must equal evals_i · col i of V)
|
||||
float residual = 0.0f;
|
||||
for (int i = 0; i < ns[c]; i++) {
|
||||
for (int r = 0; r < ns[c]; r++) {
|
||||
float Sv = 0.0f;
|
||||
for (int k = 0; k < ns[c]; k++) Sv += S_orig[r][k] * V[k][i];
|
||||
float diff = Sv - evals[i] * V[r][i];
|
||||
residual += diff * diff;
|
||||
}
|
||||
}
|
||||
residual = sqrtf(residual);
|
||||
REQUIRE_THAT(residual,
|
||||
Catch::Matchers::WithinAbs(0.0f,
|
||||
1e-2f * maxAbs[c]));
|
||||
}
|
||||
}
|
||||
|
||||
// ============================================================================
|
||||
// TEST: DeflateBidiagonal / BidiagonalIsDiagonal
|
||||
// ============================================================================
|
||||
TEST_CASE("SVD Building Block: DeflateBidiagonal and BidiagonalIsDiagonal",
|
||||
"[Matrix][SVD]") {
|
||||
float tol = 1e-8f;
|
||||
|
||||
// IsDiagonal: true on a diagonal matrix
|
||||
{
|
||||
Matrix<5, 5> W{10.0f, 0.0f, 0.0f, 0.0f, 0.0f,
|
||||
0.0f, 5.0f, 0.0f, 0.0f, 0.0f,
|
||||
0.0f, 0.0f, 2.0f, 0.0f, 0.0f,
|
||||
0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
|
||||
0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
|
||||
REQUIRE(SVD::BidiagonalIsDiagonal(W, 5, tol));
|
||||
}
|
||||
|
||||
// IsDiagonal: false when a superdiagonal is significant
|
||||
{
|
||||
Matrix<5, 5> W{10.0f, 1e-3f, 0.0f, 0.0f, 0.0f,
|
||||
0.0f, 5.0f, 0.0f, 0.0f, 0.0f,
|
||||
0.0f, 0.0f, 2.0f, 0.0f, 0.0f,
|
||||
0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
|
||||
0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
|
||||
REQUIRE_FALSE(SVD::BidiagonalIsDiagonal(W, 5, tol));
|
||||
}
|
||||
|
||||
// Deflate: small superdiagonals zeroed, significant ones kept
|
||||
{
|
||||
Matrix<5, 5> W{1.0f, 0.5f, 0.0f, 0.0f, 0.0f,
|
||||
0.0f, 2.0f, 1e-9f, 0.0f, 0.0f,
|
||||
0.0f, 0.0f, 3.0f, 0.3f, 0.0f,
|
||||
0.0f, 0.0f, 0.0f, 4.0f, 1e-12f,
|
||||
0.0f, 0.0f, 0.0f, 0.0f, 5.0f};
|
||||
SVD::DeflateBidiagonal(W, 5, tol);
|
||||
REQUIRE_THAT(W.Get(0, 1), Catch::Matchers::WithinAbs(0.5f, 1e-6f));
|
||||
REQUIRE(W.Get(1, 2) == 0.0f);
|
||||
REQUIRE_THAT(W.Get(2, 3), Catch::Matchers::WithinAbs(0.3f, 1e-6f));
|
||||
REQUIRE(W.Get(3, 4) == 0.0f);
|
||||
// Diagonal untouched
|
||||
REQUIRE_THAT(W.Get(0, 0), Catch::Matchers::WithinAbs(1.0f, 1e-6f));
|
||||
REQUIRE_THAT(W.Get(4, 4), Catch::Matchers::WithinAbs(5.0f, 1e-6f));
|
||||
// NOT fully diagonal: significant superdiagonals (0.5, 0.3) remain
|
||||
REQUIRE_FALSE(SVD::BidiagonalIsDiagonal(W, 5, tol));
|
||||
}
|
||||
|
||||
// Deflate on an already-diagonal-ish matrix makes IsDiagonal true
|
||||
{
|
||||
Matrix<5, 5> W{1.0f, 1e-9f, 0.0f, 0.0f, 0.0f,
|
||||
0.0f, 2.0f, 1e-11f, 0.0f, 0.0f,
|
||||
0.0f, 0.0f, 3.0f, 0.0f, 0.0f,
|
||||
0.0f, 0.0f, 0.0f, 4.0f, 1e-10f,
|
||||
0.0f, 0.0f, 0.0f, 0.0f, 5.0f};
|
||||
SVD::DeflateBidiagonal(W, 5, tol);
|
||||
REQUIRE(SVD::BidiagonalIsDiagonal(W, 5, tol));
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user