diff --git a/src/Matrix.cpp b/src/Matrix.cpp
index 5e087be..b1b6621 100644
--- a/src/Matrix.cpp
+++ b/src/Matrix.cpp
@@ -20,6 +20,20 @@ void EigenQR(Matrix<N, N> &matrixToDecompose, Matrix<N, N> &eigenVectors,
 #include "QR.hpp"
 #endif
 
+// Forward-declare SVD::SVD so the Matrix::SVD implementation below can call
+// it even when Matrix.cpp is pulled in through SVD.hpp's own include chain
+// (SVD.hpp -> Matrix.hpp -> Matrix.cpp), where the SVD namespace has not
+// been declared yet at this point. If we are not already inside that chain,
+// pull in the full SVD library so its template definition is available.
+namespace SVD {
+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);
+}
+#ifndef SVD_H_
+#include "SVD.hpp"
+#endif
+
 #ifdef MATRIX_H_ // since the .cpp file has to be included by the .hpp file this
                  // will evaluate to true
 #include "Matrix.hpp"
@@ -592,4 +606,18 @@ void Matrix<rows, columns>::EigenQR(Matrix<rows, rows> &eigenVectors,
   QR::EigenQR(A, eigenVectors, eigenValues, maxIterations, tolerance);
 }
 
+template <uint8_t rows, uint8_t columns>
+void Matrix<rows, columns>::SVD(Matrix<rows, columns> &U,
+                                Matrix<columns, 1> &sigma,
+                                Matrix<columns, columns> &Vt) const {
+  // Delegate to the SVD library (see src/SVD.hpp for the algorithm and
+  // conventions). NB: the fully-qualified ::SVD is required here — inside
+  // this member the unqualified name SVD refers to this method, which
+  // would shadow the namespace in a qualified lookup. SVD::SVD takes its
+  // input by non-const reference but does not modify it; pass a copy so
+  // the const-ness of *this is preserved.
+  Matrix<rows, columns> A = *this;
+  ::SVD::SVD<rows, columns>(A, U, sigma, Vt);
+}
+
 #endif // MATRIX_H_
\ No newline at end of file
diff --git a/src/Matrix.hpp b/src/Matrix.hpp
index 5422898..aa0e909 100644
--- a/src/Matrix.hpp
+++ b/src/Matrix.hpp
@@ -6,7 +6,6 @@
 #include <type_traits>
 
 // TODO: Add a function to compute RREF
-// TODO: Add a function for SVD decomposition
 // TODO: Add a function for LQ decomposition
 
 template <uint8_t rows, uint8_t columns> class Matrix {
@@ -253,6 +252,24 @@ public:
   void EigenQR(Matrix<rows, rows> &eigenVectors, Matrix<rows, 1> &eigenValues,
                uint32_t maxIterations = 1000, float tolerance = 1e-6f) const;
 
+  /**
+   * @brief Compute the Singular Value Decomposition (SVD) of this matrix.
+   *
+   * Wrapper around SVD::SVD (see SVD.hpp for the full algorithm
+   * description, output storage conventions, and stack-usage notes).
+   * Decomposes A = U · Σ · Vᵀ where U is rows×columns, Σ is the vector
+   * of singular values (columns×1, sorted descending), and Vᵀ is
+   * columns×columns. Works for any shape (wide matrices are handled
+   * internally by computing SVD(Aᵀ) and swapping the factors back).
+   * This matrix is not modified.
+   *
+   * @param U Output: left singular vectors (rows×columns)
+   * @param sigma Output: singular values in descending order (columns×1)
+   * @param Vt Output: right singular vectors, transposed (columns×columns)
+   */
+  void SVD(Matrix<rows, columns> &U, Matrix<columns, 1> &sigma,
+           Matrix<columns, columns> &Vt) const;
+
 protected:
   std::array<float, rows * columns> matrix;
 
diff --git a/unit-tests/svd-integration-test.cpp b/unit-tests/svd-integration-test.cpp
index d05f81c..48559e3 100644
--- a/unit-tests/svd-integration-test.cpp
+++ b/unit-tests/svd-integration-test.cpp
@@ -250,3 +250,114 @@ TEST_CASE("SVD Integration: symmetric positive definite 2x2 [[5,3],[3,5]]",
   std::cout << "Sigma: [" << sigma.Get(0, 0) << ", " << sigma.Get(1, 0)
             << "]\n";
 }
+
+// ----------------------------------------------------------------------------
+// Matrix::SVD member wrapper (delegates to SVD::SVD)
+// ----------------------------------------------------------------------------
+
+/**
+ * Reconstruction error ‖U·diag(sigma)·Vᵀ − A‖_F. Zero-padded entries of
+ * U/sigma/Vt (wide/tall cases) are zero by the output conventions, so the
+ * full product equals U[:, :k]·diag(sigma[:k])·Vt[:k, :].
+ */
+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) {
+  Matrix<rows, columns> recon{0};
+  Matrix<rows, columns> Usig{0};
+  for (int i = 0; i < rows; i++)
+    for (int j = 0; j < columns; 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 < rows; i++)
+    for (int j = 0; j < columns; j++) {
+      float diff = recon.Get(i, j) - A.Get(i, j);
+      err += diff * diff;
+    }
+  return sqrtf(err);
+}
+
+/**
+ * Orthonormality of the first k columns of M: the k×k leading block of
+ * MᵀM must equal I_k. (For a tall SVD, U has k = min(rows, cols)
+ * meaningful columns and this is the full UᵀU.)
+ */
+template <uint8_t r, uint8_t c>
+static bool leadingColumnsOrthonormal(const Matrix<r, c> &M, uint8_t k,
+                                      float tol = 1e-4f) {
+  Matrix<c, r> Mt = M.Transpose();
+  Matrix<c, c> MtM{0};
+  Mt.Mult(M, MtM);
+  for (int i = 0; i < k; i++)
+    for (int j = 0; j < k; j++) {
+      float expected = (i == j) ? 1.0f : 0.0f;
+      if (fabsf(MtM.Get(i, j) - expected) > tol)
+        return false;
+    }
+  return true;
+}
+
+/**
+ * Orthonormality of the first k rows of M: the k×k leading block of
+ * M·Mᵀ must equal I_k. (Vᵀ may have zero-padded trailing rows in the
+ * wide case, so check only the meaningful leading block.)
+ */
+template <uint8_t r, uint8_t c>
+static bool leadingRowsOrthonormal(const Matrix<r, c> &M, uint8_t k,
+                                   float tol = 1e-4f) {
+  Matrix<c, r> Mt = M.Transpose();
+  Matrix<r, r> MMt{0};
+  M.Mult(Mt, MMt);
+  for (int i = 0; i < k; i++)
+    for (int j = 0; j < k; j++) {
+      float expected = (i == j) ? 1.0f : 0.0f;
+      if (fabsf(MMt.Get(i, j) - expected) > tol)
+        return false;
+    }
+  return true;
+}
+
+TEST_CASE("Matrix::SVD wrapper: 3x2 tall [[1,2],[3,4],[5,6]]",
+          "[Matrix][SVD][Wrapper]") {
+  Matrix<3, 2> A{1, 2, 3, 4, 5, 6};
+  Matrix<3, 2> U{0};
+  Matrix<2, 1> sigma{0};
+  Matrix<2, 2> Vt{0};
+
+  A.SVD(U, sigma, Vt);
+
+  // Reference singular values from numpy: [9.52552, 0.514301]
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(9.52552f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(0.514301f, 1e-3f));
+
+  REQUIRE(leadingColumnsOrthonormal(U, 2));
+  REQUIRE(leadingRowsOrthonormal(Vt, 2));
+
+  float err = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+}
+
+TEST_CASE("Matrix::SVD wrapper: 2x3 wide [[1,2,3],[4,5,6]]",
+          "[Matrix][SVD][Wrapper]") {
+  Matrix<2, 3> A{1, 2, 3, 4, 5, 6};
+  Matrix<2, 3> U{0};
+  Matrix<3, 1> sigma{0};
+  Matrix<3, 3> Vt{0};
+
+  A.SVD(U, sigma, Vt);
+
+  // Reference singular values from numpy: [9.50803, 0.77287]; the third
+  // entry (wide-matrix padding) must be zero.
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(9.50803f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(0.77287f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+
+  REQUIRE(leadingColumnsOrthonormal(U, 2));
+  REQUIRE(leadingRowsOrthonormal(Vt, 2));
+
+  float err = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+}