diff --git a/unit-tests/matrix-tests.cpp b/unit-tests/matrix-tests.cpp
index 90c0956..f32035a 100644
--- a/unit-tests/matrix-tests.cpp
+++ b/unit-tests/matrix-tests.cpp
@@ -916,6 +916,132 @@ TEST_CASE("SVD: 5×5 Symmetric Tridiagonal", "Matrix") {
   REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 }
 
+TEST_CASE("SVD: 5×5 Full-Rank Random", "Matrix") {
+  // Reference: scipy.linalg.svd, np.random.default_rng(7).standard_normal((5,5))
+  // cond ≈ 11.5
+  // σ = [3.04651784, 2.22681732, 1.84290662, 1.02101969, 0.264826749]
+  Matrix<5, 5> A{0.00123015f, 0.298746f, -0.274138f, -0.890592f, -0.454671f,
+                 -0.991647f, 0.0601436f, 1.34022f, -0.492207f, -0.620475f,
+                 0.489842f, 0.356887f, 0.105414f, -0.930468f, -0.0292518f,
+                 0.695303f, -1.34421f, -0.457616f, -1.90122f, -1.28954f,
+                 -1.84174f, -0.235091f, -1.26745f, 0.271264f, 0.156751f};
+  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.04651784f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(2.22681732f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(1.84290662f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(3, 0),
+               Catch::Matchers::WithinRel(1.02101969f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(4, 0), Catch::Matchers::WithinRel(0.264826749f, 1e-4f));
+
+  REQUIRE(isSortedDescending(sigma, 5));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+}
+
+TEST_CASE("SVD: 5×5 Rank-Deficient (rank 3)", "Matrix") {
+  // Reference: scipy.linalg.svd of rng.standard_normal((5,3)) @
+  // rng.standard_normal((3,5)) — exactly rank 3
+  // σ = [7.29829771, 2.76487864, 1.57392325, ~1e-16, ~1e-16]
+  Matrix<5, 5> A{3.46252f, -1.52873f, -0.111526f, 1.28954f, -5.18688f,
+                 -1.34702f, 1.92936f, -0.0410797f, -0.958791f, 0.449623f,
+                 0.844592f, 0.0986352f, 0.408213f, 0.124867f, -2.45393f,
+                 1.34177f, -0.587312f, -1.39847f, 0.580032f, -0.167692f,
+                 0.915462f, -0.311165f, -1.16141f, 0.377283f, 0.055373f};
+  Matrix<5, 5> U{}, Vt{};
+  Matrix<5, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(7.29829771f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(2.76487864f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(1.57392325f, 1e-4f));
+  // The two rank-deficient singular values must be at noise level
+  REQUIRE(sigma.Get(3, 0) < 1e-3f);
+  REQUIRE(sigma.Get(4, 0) < 1e-3f);
+
+  REQUIRE(isSortedDescending(sigma, 5));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  // Rank-3 matrix: the top-3 SVD terms must reproduce A
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 5e-3f));
+}
+
+TEST_CASE("SVD: 5×5 Wide Dynamic Range (cond ≈ 9000)", "Matrix") {
+  // Symmetric banded, diagonal decays 50 → 1e-3, off-diagonals 3 → 0.01
+  // Reference: scipy.linalg.svd
+  // σ = [50.1496395, 10.1719501, 1.02846061, 0.0207253626, 0.00557535757]
+  Matrix<5, 5> A{50.0f, 3.0f, 0.0f, 0.0f, 0.0f,
+                 -3.0f, 10.0f, 0.5f, 0.0f, 0.0f,
+                 0.0f, -0.5f, 1.0f, 0.08f, 0.0f,
+                 0.0f, 0.0f, -0.08f, 0.01f, 0.01f,
+                 0.0f, 0.0f, 0.0f, -0.01f, 0.001f};
+  Matrix<5, 5> U{}, Vt{};
+  Matrix<5, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0),
+               Catch::Matchers::WithinRel(50.1496395f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0),
+               Catch::Matchers::WithinRel(10.1719501f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0),
+               Catch::Matchers::WithinRel(1.02846061f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(3, 0),
+               Catch::Matchers::WithinRel(0.0207253626f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(4, 0),
+               Catch::Matchers::WithinRel(0.00557535757f, 1e-3f));
+
+  REQUIRE(isSortedDescending(sigma, 5));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  // Relative to the largest entry (‖A‖F ≈ 50.1)
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 5e-2f));
+}
+
+TEST_CASE("SVD: 5×5 Symmetric Indefinite", "Matrix") {
+  // Symmetric with negative eigenvalues — σ must equal |eigenvalues|
+  // Reference: scipy.linalg.svd
+  // σ = [4.70141723, 3.76392521, 2.73426097, 1.15773083, 0.642665756]
+  Matrix<5, 5> A{2.0f, -1.0f, 0.0f, 0.0f, 0.5f,
+                 -1.0f, 2.0f, -1.0f, 0.0f, 0.0f,
+                 0.0f, -1.0f, 3.0f, -1.0f, 0.0f,
+                 0.0f, 0.0f, -1.0f, 2.0f, -1.0f,
+                 0.5f, 0.0f, 0.0f, -1.0f, 4.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(4.70141723f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0),
+               Catch::Matchers::WithinRel(3.76392521f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0),
+               Catch::Matchers::WithinRel(2.73426097f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(3, 0),
+               Catch::Matchers::WithinRel(1.15773083f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(4, 0),
+               Catch::Matchers::WithinRel(0.642665756f, 1e-4f));
+
+  REQUIRE(isSortedDescending(sigma, 5));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+}
+
 TEST_CASE("SVD: 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]