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 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.

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:

brew install cmake eigen

Use this project structure:

svd-learning/
├── CMakeLists.txt
└── src/
    └── main.cpp

Create CMakeLists.txt:

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:

cmake -S . -B build \
  -DCMAKE_PREFIX_PATH="$(brew --prefix eigen)"

cmake --build build

Run it:

./build/svd_learning

How the Mathematics Maps to Eigen

Mathematical step	Eigen expression
Define A	`Eigen::Matrix<double, 3, 2> 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`

#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <Eigen/Geometry>
#include <Eigen/SVD>

#include <cmath>
#include <iostream>

int main()
{
    // A is a 3 x 2 matrix.
    Eigen::Matrix<double, 3, 2> 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::Matrix2d> 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<double, 3, 2> Sigma =
        Eigen::Matrix<double, 3, 2>::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<double, 3, 2> 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<Eigen::Matrix<double, 3, 2>> 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<double, 3, 2> Sigma_eigen =
        Eigen::Matrix<double, 3, 2>::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<double, 3, 2> 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

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! 😊

Introduction

Finding the Singular Value Decomposition

Form a Square Symmetric Matrix

Find the Eigenvalues

Calculate the Singular Values

Find the Right Singular Vectors

Construct the Sigma Matrix

Find the Left Singular Vectors

Assemble and Verify the Decomposition

Worked Example

Calculate AᵀA

Find the Eigenvalues

Find the Right Singular Vectors

Calculate the Singular Values and Construct Σ

Find the Left Singular Vectors

Write the Final SVD

What the Two Patterns Mean

Why Software May Display Different Signs

What If the Matrix Has More Columns Than Rows?

Full SVD and Compact SVD

Short Formula Sheet

C++ Implementation with Eigen

Project Setup on macOS

How the Mathematics Maps to Eigen

Complete main.cpp

Expected Key Results

Conclusion

Complete `main.cpp`