diff --git a/src/SVD.cpp b/src/SVD.cpp index 4c2b2b4..47990ec 100644 --- a/src/SVD.cpp +++ b/src/SVD.cpp @@ -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; } diff --git a/src/SVD.hpp b/src/SVD.hpp index 1477097..7c88362 100644 --- a/src/SVD.hpp +++ b/src/SVD.hpp @@ -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 diff --git a/unit-tests/matrix-tests.cpp b/unit-tests/matrix-tests.cpp index 8390cfa..90c0956 100644 --- a/unit-tests/matrix-tests.cpp +++ b/unit-tests/matrix-tests.cpp @@ -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{}; diff --git a/unit-tests/svd-build-blocks-tests.cpp b/unit-tests/svd-build-blocks-tests.cpp index da19f41..37932d5 100644 --- a/unit-tests/svd-build-blocks-tests.cpp +++ b/unit-tests/svd-build-blocks-tests.cpp @@ -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)); + } +}