# Singular Value Decomposition: From the Math to C++ with Eigen > Introduction Imagine a cinema recording how many tickets it sells for each film every day. The rows of its table are days, the columns are films, and each number is a ticket count. Suppose Film A sell - Published: 2026-09-30T07:36:28.361Z - Reading time: 15 minutes - Tags: #svd, linear algebra , C++, eigen, Mathematics, Machine Learning - Original publication: [Hashnode](https://abulasar.com/singular-value-decomposition-cpp-eigen) - Canonical portfolio page: [abulasar.dev](https://abulasar.dev/notes/singular-value-decomposition-cpp-eigen/) ## Introduction Imagine a cinema recording how many tickets it sells for each film every day. The rows of its table are days, the columns are films, and each number is a ticket count. Suppose Film A sells 100 tickets on Friday and 200 on Saturday. Film B sells 200 on Friday and 400 on Saturday. The table contains four numbers, but one simple pattern explains them: Saturday’s sales are twice Friday’s, and Film B’s sales are twice Film A’s. In a table with thousands of days and films, patterns like this would be harder to spot. **Singular Value Decomposition (SVD)** helps us find them. It describes a matrix as a collection of patterns, ordered by how strongly each one contributes. We can use all the patterns to rebuild the original matrix, or keep only the strongest ones to make an approximation. For a matrix **A**, SVD is written as $$A = U\Sigma V^T$$ Think of the three parts as answering three questions: * **U:** How much does each row take part in a pattern? In our example, the rows are days. * **Vᵀ:** How much does each column take part in the matching pattern? Here, the columns are films. * **Σ:** How strong is that pattern? The first column of **U**, the first singular value in **Σ**, and the first row of **Vᵀ** work together to describe one pattern. The second column, value, and row describe another. Adding all of them gives us **A**. This does not mean SVD knows what a “day” or a “film” is. It works with the numbers in the matrix. The ticket table is simply a way to understand what its three parts do. Eigenvalue decomposition helps us study certain directions of a square matrix. SVD gives us a related way to study **any real matrix**, including a rectangular table with a different number of rows and columns. If $$A \in \mathbb{R}^{m \times n},$$ then the matrices in the **full SVD** have these dimensions: $$U \in \mathbb{R}^{m \times m}, \qquad \Sigma \in \mathbb{R}^{m \times n}, \qquad V \in \mathbb{R}^{n \times n}.$$ The columns of **U** and **V** are called *left* and *right singular vectors*. The nonnegative numbers on the diagonal of **Σ** are called *singular values*. Both **U** and **V** are orthogonal matrices: $$U^TU = I_m, \qquad V^TV = I_n.$$ In the rest of this article, we will calculate these three matrices step by step, then implement the same calculation in C++ using Eigen. ## Finding the Singular Value Decomposition Let $$A \in \mathbb{R}^{m \times n}.$$ We will begin with the case **m ≥ n**, where **A** has at least as many rows as columns. ### Form a Square Symmetric Matrix Calculate $$M = A^TA.$$ Because **A** is an **m × n** matrix, $$A^TA \in \mathbb{R}^{n \times n}.$$ The matrix **AᵀA** is square, symmetric, and positive semidefinite. To see why its eigenvalues cannot be negative, take any vector **x**: $$x^T(A^TA)x=(Ax)^T(Ax)=\left\|Ax\right\|_2^2\geq0$$ The eigenvectors of **AᵀA** will give us the columns of **V**, while its eigenvalues will give us the squared singular values. *Note: We find the eigenvalues of AᵀA, not the eigenvalues of A itself.* ### Find the Eigenvalues Solve the characteristic equation $$\det(A^TA - \lambda I) = 0.$$ For **m ≥ n**, this gives **n** eigenvalues. Arrange them from largest to smallest: $$\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_n \ge 0.$$ The eigenvalues may be repeated, and some of them may be zero. ### Calculate the Singular Values The singular values are the square roots of the eigenvalues: $$\sigma_i = \sqrt{\lambda_i}.$$ Therefore, $$\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_n \ge 0.$$ The number of positive singular values is equal to the rank of **A**. *Note: Singular values are conventionally written in decreasing order, not increasing order.* ### Find the Right Singular Vectors For each eigenvalue **λᵢ**, solve $$(A^TA - \lambda_i I)v_i = 0.$$ Normalize the resulting eigenvector: $$v_i\leftarrow\frac{v_i}{\left\|v_i\right\|_2}$$ Place the normalized eigenvectors into the columns of **V**: $$V = \begin{bmatrix}v_1 & v_2 & \cdots & v_n\end{bmatrix}.$$ The order must remain consistent: **vᵢ** must correspond to **λᵢ** and **σᵢ**. The columns of **V** must be orthonormal: $$V^TV = I_n.$$ If an eigenvalue is repeated, choose an orthonormal basis for its eigenspace. ### Construct the Sigma Matrix The matrix **Σ** has the same dimensions as **A**: $$\Sigma \in \mathbb{R}^{m \times n}.$$ Place the singular values on its main diagonal in decreasing order and set every other entry to zero. For a **3 × 2** matrix, $$\Sigma = \begin{bmatrix} \sigma_1 & 0 \\\\ 0 & \sigma_2 \\\\ 0 & 0 \end{bmatrix}.$$ ### Find the Left Singular Vectors For each nonzero singular value, calculate $$u_i = \frac{Av_i}{\sigma_i}.$$ Because **σᵢ = √λᵢ**, the same formula can be written as $$u_i = \frac{Av_i}{\sqrt{\lambda_i}}.$$ Each **uᵢ** must use the matching **vᵢ** and **σᵢ**. For the full SVD, **U** must contain **m** orthonormal columns. If the formula above does not produce all **m** columns, complete **U** with additional orthonormal vectors. For a **3 × 2** full-rank matrix, the first two columns are obtained from the formula, and a third column can be found using $$u_3 = u_1 \times u_2.$$ Normalize **u₃** if necessary. The cross product works only in three-dimensional space. In higher dimensions, use a method such as Gram–Schmidt orthogonalization or find an orthonormal basis for the null space of **Aᵀ**. Finally, $$U = \begin{bmatrix}u_1 & u_2 & \cdots & u_m\end{bmatrix},$$ with $$U^TU = I_m.$$ *Important: If σᵢ = 0, do not calculate Avᵢ/σᵢ because that would require division by zero. Complete the missing columns of U using an orthonormal basis instead.* ### Assemble and Verify the Decomposition The final decomposition is $$A = U\Sigma V^T.$$ Notice that the final factor is **Vᵀ**, not **V**. Verify all three conditions: $$U^TU = I_m,$$ $$V^TV = I_n,$$ and $$U\Sigma V^T = A.$$ Small differences are normal when floating-point values or rounded decimal values are used. ## Worked Example We will calculate the SVD of the same matrix that we later use in the C++ program: $$A = \begin{bmatrix} 3 & 1 \\\\ 2 & 2 \\\\ 1 & 3 \end{bmatrix}.$$ Here, **m = 3** and **n = 2**, so the full SVD has the dimensions $$U_{3 \times 3}, \qquad \Sigma_{3 \times 2}, \qquad V_{2 \times 2}.$$ ### Calculate AᵀA First, $$A^T = \begin{bmatrix} 3 & 2 & 1 \\\\ 1 & 2 & 3 \end{bmatrix}.$$ Therefore, $$A^TA = \begin{bmatrix} 3 & 2 & 1 \\\\ 1 & 2 & 3 \end{bmatrix} \begin{bmatrix} 3 & 1 \\\\ 2 & 2 \\\\ 1 & 3 \end{bmatrix} = \begin{bmatrix} 14 & 10 \\\\ 10 & 14 \end{bmatrix}.$$ ### Find the Eigenvalues Solve $$\det(A^TA - \lambda I) = 0.$$ For this matrix, $$\begin{aligned} \det \begin{bmatrix} 14-\lambda & 10 \\\\ 10 & 14-\lambda \end{bmatrix} &= (14-\lambda)^2 - 100 \\\\ &= \lambda^2 - 28\lambda + 96 \\\\ &= (\lambda-24)(\lambda-4). \end{aligned}$$ Thus, in decreasing order, $$\lambda_1 = 24, \qquad \lambda_2 = 4.$$ ### Find the Right Singular Vectors For **λ₁ = 24**, solving $$(A^TA-24I)v_1=0$$ gives **x = y**. After normalization, choose $$v_1 = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\\\ 1 \end{bmatrix}.$$ For **λ₂ = 4**, solving $$(A^TA-4I)v_2=0$$ gives **x = −y**. After normalization, choose $$v_2 = \frac{1}{\sqrt{2}} \begin{bmatrix} -1 \\\\ 1 \end{bmatrix}.$$ Therefore, $$V = \begin{bmatrix}v_1 & v_2\end{bmatrix} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & -1 \\\\ 1 & 1 \end{bmatrix}.$$ ### Calculate the Singular Values and Construct Σ The singular values are $$\sigma_1 = \sqrt{24} = 2\sqrt{6} \approx 4.89898,$$ $$\sigma_2 = \sqrt{4} = 2.$$ Thus, $$\Sigma = \begin{bmatrix} 2\sqrt{6} & 0 \\\\ 0 & 2 \\\\ 0 & 0 \end{bmatrix}.$$ ### Find the Left Singular Vectors Calculate **u₁** using the matching pair **v₁** and **σ₁**: $$\begin{aligned} u_1 &= \frac{Av_1}{\sigma_1} \\\\ &= \frac{1}{2\sqrt{6}} \begin{bmatrix} 3 & 1 \\\\ 2 & 2 \\\\ 1 & 3 \end{bmatrix} \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\\\ 1 \end{bmatrix} \\\\ &= \frac{1}{\sqrt{3}} \begin{bmatrix} 1 \\\\ 1 \\\\ 1 \end{bmatrix}. \end{aligned}$$ Next, $$\begin{aligned} u_2 &= \frac{Av_2}{\sigma_2} \\\\ &= \frac{1}{2} \begin{bmatrix} 3 & 1 \\\\ 2 & 2 \\\\ 1 & 3 \end{bmatrix} \frac{1}{\sqrt{2}} \begin{bmatrix} -1 \\\\ 1 \end{bmatrix} \\\\ &= \frac{1}{\sqrt{2}} \begin{bmatrix} -1 \\\\ 0 \\\\ 1 \end{bmatrix}. \end{aligned}$$ To complete the full **3 × 3** matrix **U**, calculate $$u_3 = u_1 \times u_2 = \frac{1}{\sqrt{6}} \begin{bmatrix} 1 \\\\ -2 \\\\ 1 \end{bmatrix}.$$ Therefore, $$U = \begin{bmatrix} \frac{1}{\sqrt{3}} & -\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{6}} \\\\ \frac{1}{\sqrt{3}} & 0 & -\frac{2}{\sqrt{6}} \\\\ \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{6}} \end{bmatrix}.$$ ### Write the Final SVD The decomposition is $$A = U\Sigma V^T,$$ where $$V^T = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\\\ -1 & 1 \end{bmatrix}.$$ Multiplying the factors reconstructs the original matrix: $$U\Sigma V^T = \begin{bmatrix} 3 & 1 \\\\ 2 & 2 \\\\ 1 & 3 \end{bmatrix} = A.$$ ### What the Two Patterns Mean We can also write the SVD as the sum of two simpler matrices: $$A=\sigma_1u_1v_1^T+\sigma_2u_2v_2^T =\begin{bmatrix}2&2\\\\2&2\\\\2&2\end{bmatrix} +\begin{bmatrix}1&-1\\\\0&0\\\\-1&1\end{bmatrix}.$$ Imagine the rows are three days and the columns are two films, with ticket counts measured in hundreds. The first part gives us a simple starting point: **200 tickets for each film on each day**. The second part adjusts those numbers. On the first day, Film A gets 100 more tickets and Film B gets 100 fewer. On the second day, nothing changes. On the third day, the adjustment is reversed. Adding both parts gives us the exact original table. If we keep only the first part, we get a simpler approximation in which every value is 2. This is what it means to keep the strongest SVD pattern and leave out the smaller one. ### Why Software May Display Different Signs An eigenvector can be multiplied by **−1** and still be an eigenvector. Consequently, singular vectors are not unique in sign. For any singular-vector pair, replacing both vectors by their negatives does not change the decomposition: $$(-u_i)\sigma_i(-v_i)^T = u_i\sigma_i v_i^T.$$ This is why our manual result and Eigen's `JacobiSVD` result may show opposite signs in corresponding columns of **U** and **V**. Both results are correct if they reconstruct **A**. ## What If the Matrix Has More Columns Than Rows? When **m < n**, the matrix **AAᵀ** is smaller than **AᵀA**, so it is often more convenient to start with $$M = AA^T.$$ Its eigenvectors give the left singular vectors directly: $$AA^T u_i = \lambda_i u_i.$$ The singular values are still $$\sigma_i = \sqrt{\lambda_i}.$$ For each nonzero singular value, obtain the matching right singular vector from $$v_i = \frac{A^T u_i}{\sigma_i}.$$ Then complete **V** with additional orthonormal vectors if the full SVD is required. In summary: * Starting with **AᵀA** gives **V** first; then calculate **U**. * Starting with **AAᵀ** gives **U** first; then calculate **V**. Both approaches are mathematically valid. Choosing the smaller square matrix merely reduces the amount of work. The nonzero eigenvalues of **AᵀA** and **AAᵀ** are the same. ## Full SVD and Compact SVD The full SVD is $$A_{m \times n} = U_{m \times m} \Sigma_{m \times n} V^T_{n \times n}.$$ If **A** has rank **r**, keeping only the singular vectors associated with its positive singular values gives the compact SVD: $$A = U_r\Sigma_rV_r^T,$$ where $$U_r \in \mathbb{R}^{m \times r}, \qquad \Sigma_r \in \mathbb{R}^{r \times r}, \qquad V_r \in \mathbb{R}^{n \times r}.$$ The compact SVD reconstructs **A** exactly while omitting directions associated with zero singular values. ## Short Formula Sheet When starting with **AᵀA**: $$\boxed{A^TA v_i = \sigma_i^2 v_i}$$ $$\boxed{u_i = \frac{Av_i}{\sigma_i}}$$ When starting with **AAᵀ**: $$\boxed{AA^T u_i = \sigma_i^2 u_i}$$ $$\boxed{v_i = \frac{A^T u_i}{\sigma_i}}$$ In both cases: $$\boxed{\sigma_i = \sqrt{\lambda_i}}$$ Finally: $$\boxed{A = U\Sigma V^T}$$ ## C++ Implementation with Eigen The C++ program below deliberately performs the mathematical construction first. It then calls Eigen's built-in `JacobiSVD` implementation so that we can compare the results. ### Project Setup on macOS Install CMake and Eigen with Homebrew: ```bash brew install cmake eigen ``` Use this project structure: ```text svd-learning/ ├── CMakeLists.txt └── src/ └── main.cpp ``` Create `CMakeLists.txt`: ```cmake cmake_minimum_required(VERSION 3.20) project(svd_learning LANGUAGES CXX) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) find_package(Eigen3 3.3 REQUIRED NO_MODULE) add_executable(svd_learning src/main.cpp) target_link_libraries(svd_learning PRIVATE Eigen3::Eigen) ``` Configure and build the project: ```bash cmake -S . -B build \ -DCMAKE_PREFIX_PATH="$(brew --prefix eigen)" cmake --build build ``` Run it: ```bash ./build/svd_learning ``` ### How the Mathematics Maps to Eigen | Mathematical step | Eigen expression | | --- | --- | | Define **A** | `Eigen::Matrix A` | | Calculate **AᵀA** | `A.transpose() * A` | | Find eigenpairs of **AᵀA** | `Eigen::SelfAdjointEigenSolver` | | Build **V** | Place the ordered eigenvectors into its columns | | Calculate **σᵢ** | `std::sqrt(lambda(i))` | | Calculate **uᵢ** | `(A * V.col(i)) / sigma(i)` | | Complete **U** in **R³** | `u1.cross(u2)` | | Verify the factorization | `(A - U * Sigma * V.transpose()).norm()` | | Compute SVD directly | `Eigen::JacobiSVD` | `SelfAdjointEigenSolver` returns eigenvalues in increasing order, so the program reverses the two eigenpairs before constructing **V** and **Σ**. ### Complete `main.cpp` ```cpp #include #include #include #include #include #include int main() { // A is a 3 x 2 matrix. Eigen::Matrix A; A << 3, 1, 2, 2, 1, 3; std::cout << "Matrix A:\n" << A << '\n'; // Step 1: Form A^T A. const Eigen::Matrix2d A_transpose_A = A.transpose() * A; std::cout << "\nA^T * A:\n" << A_transpose_A << '\n'; // Step 2: Find the eigenvalues and normalized eigenvectors. const Eigen::SelfAdjointEigenSolver eigen_solver( A_transpose_A ); if (eigen_solver.info() != Eigen::Success) { std::cerr << "Eigenvalue decomposition failed.\n"; return 1; } const Eigen::Vector2d eigenvalues_ascending = eigen_solver.eigenvalues(); const Eigen::Matrix2d eigenvectors_ascending = eigen_solver.eigenvectors(); std::cout << "\nEigenvalues returned by Eigen:\n" << eigenvalues_ascending.transpose() << '\n'; std::cout << "\nCorresponding normalized eigenvectors:\n" << eigenvectors_ascending << '\n'; // Eigen returns the eigenvalues in increasing order. // Reorder both eigenvalues and matching eigenvectors for SVD. Eigen::Vector2d lambda; lambda << eigenvalues_ascending(1), eigenvalues_ascending(0); Eigen::Matrix2d V; V.col(0) = eigenvectors_ascending.col(1); V.col(1) = eigenvectors_ascending.col(0); std::cout << "\nEigenvalues in decreasing SVD order:\n" << lambda.transpose() << '\n'; std::cout << "\nMatrix V:\n" << V << '\n'; std::cout << "\nV^T * V:\n" << V.transpose() * V << '\n'; // Step 3: Singular values are square roots of eigenvalues. Eigen::Vector2d sigma; sigma << std::sqrt(lambda(0)), std::sqrt(lambda(1)); std::cout << "\nSingular values:\n" << sigma.transpose() << '\n'; // Step 4: Construct the 3 x 2 Sigma matrix. Eigen::Matrix Sigma = Eigen::Matrix::Zero(); Sigma(0, 0) = sigma(0); Sigma(1, 1) = sigma(1); std::cout << "\nMatrix Sigma:\n" << Sigma << '\n'; // Step 5: Calculate the first two columns of U. const Eigen::Vector3d u1 = (A * V.col(0)) / sigma(0); const Eigen::Vector3d u2 = (A * V.col(1)) / sigma(1); std::cout << "\nVector u1:\n" << u1 << '\n'; std::cout << "\nVector u2:\n" << u2 << '\n'; std::cout << "\nNorm of u1: " << u1.norm() << '\n'; std::cout << "Norm of u2: " << u2.norm() << '\n'; std::cout << "Dot product u1^T * u2: " << u1.dot(u2) << '\n'; // Complete U. The cross product is specific to three dimensions. Eigen::Vector3d u3 = u1.cross(u2); u3.normalize(); Eigen::Matrix3d U; U.col(0) = u1; U.col(1) = u2; U.col(2) = u3; std::cout << "\nVector u3:\n" << u3 << '\n'; std::cout << "\nMatrix U:\n" << U << '\n'; std::cout << "\nU^T * U:\n" << U.transpose() * U << '\n'; // Step 6: Reconstruct A and measure the numerical error. const Eigen::Matrix A_reconstructed = U * Sigma * V.transpose(); std::cout << "\nV^T:\n" << V.transpose() << '\n'; std::cout << "\nReconstructed A = U * Sigma * V^T:\n" << A_reconstructed << '\n'; const double reconstruction_error = (A - A_reconstructed).norm(); constexpr double tolerance = 1e-10; std::cout << "\nReconstruction error: " << reconstruction_error << '\n'; std::cout << "Reconstruction valid: " << std::boolalpha << (reconstruction_error < tolerance) << '\n'; // Step 7: Compare with Eigen's direct SVD implementation. const Eigen::JacobiSVD> jacobi_svd( A, Eigen::ComputeFullU | Eigen::ComputeFullV ); if (jacobi_svd.info() != Eigen::Success) { std::cerr << "JacobiSVD computation failed.\n"; return 1; } const Eigen::Matrix3d U_eigen = jacobi_svd.matrixU(); const Eigen::Vector2d sigma_eigen = jacobi_svd.singularValues(); const Eigen::Matrix2d V_eigen = jacobi_svd.matrixV(); Eigen::Matrix Sigma_eigen = Eigen::Matrix::Zero(); Sigma_eigen.diagonal() = sigma_eigen; std::cout << "\n--- Eigen JacobiSVD result ---\n"; std::cout << "\nEigen singular values:\n" << sigma_eigen.transpose() << '\n'; std::cout << "\nEigen matrix U:\n" << U_eigen << '\n'; std::cout << "\nEigen matrix Sigma:\n" << Sigma_eigen << '\n'; std::cout << "\nEigen matrix V:\n" << V_eigen << '\n'; const Eigen::Matrix A_reconstructed_eigen = U_eigen * Sigma_eigen * V_eigen.transpose(); const double eigen_reconstruction_error = (A - A_reconstructed_eigen).norm(); std::cout << "\nA reconstructed using JacobiSVD:\n" << A_reconstructed_eigen << '\n'; std::cout << "\nJacobiSVD reconstruction error: " << eigen_reconstruction_error << '\n'; std::cout << "Singular values match: " << sigma.isApprox(sigma_eigen, tolerance) << '\n'; return 0; } ``` ### Expected Key Results The program should report ```text Eigenvalues in decreasing SVD order: 24 4 Singular values: 4.89898 2 Reconstructed A = U * Sigma * V^T: 3 1 2 2 1 3 Reconstruction valid: true Singular values match: true ``` ## Conclusion The central relationships behind SVD are $$A^TA v_i = \sigma_i^2 v_i,$$ $$Av_i = \sigma_i u_i,$$ and $$A = U\Sigma V^T.$$ In this example, we derived every matrix by hand, reconstructed the original matrix in C++, and confirmed the result using Eigen's built-in `JacobiSVD` implementation. If you have a question or notice something that could be explained more clearly, please leave a comment. Thanks for reading! 😊 ## Continue - [All technical notes](https://abulasar.dev/notes/) - [About Abulasar](https://abulasar.dev/about/) - [Contact Abulasar](https://abulasar.dev/contact/)