Singular Value Decomposition (SVD) in C++: A Comprehensive Guide

Introduction

Singular Value Decomposition (SVD) is a fundamental matrix factorization method that decomposes a matrix into three components: \( U \), \( \Sigma \), and \( V^T \). This decomposition has numerous applications in data compression, dimensionality reduction, and solving linear systems. In this guide, we’ll explore how to implement and use SVD effectively in C++.

Mathematical Background

Singular Value Decomposition (SVD) is a mathematical method for breaking a matrix into three key components that reveal its structure and properties. This decomposition is widely used because it helps simplify complex problems and provides insights into the matrix’s behavior. Let’s take a closer look at the formula and its components.

The SVD Formula

For any matrix \( A \in \mathbb{R}^{m \times n} \), SVD represents \( A \) as:

Here’s what each component means:

• \( U \in \mathbb{R}^{m \times m} \): An orthogonal matrix whose columns are called the left singular vectors of \( A \).
• \( \Sigma \in \mathbb{R}^{m \times n} \): A diagonal matrix containing the singular values of \( A \), which are non-negative and sorted in descending order.
• \( V^T \in \mathbb{R}^{n \times n} \): The transpose of an orthogonal matrix, where the columns of \( V \) are the right singular vectors of \( A \).

Understanding the Components

Each part of the decomposition plays a critical role:

1. Singular Values (\( \Sigma \))

• The singular values (\( \sigma_i \)) are the diagonal entries of \( \Sigma \) and represent the “strength” or importance of each dimension in \( A \).
• They are ordered as \( \sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_r \geq 0 \), where \( r \) is the rank of \( A \).

2. Orthogonal Matrices (\( U \) and \( V \))

• \( U \): Columns of \( U \) are eigenvectors of \( A A^T \) and are orthogonal (\( U^T U = I \)).
• \( V \): Columns of \( V \) are eigenvectors of \( A^T A \) and are also orthogonal (\( V^T V = I \)).

Computing SVD

The process of computing SVD involves the following steps:

1. Form \( A A^T \) and \( A^T A \), which are square matrices used to compute eigenvalues and eigenvectors.
2. Find the eigenvalues of \( A A^T \) and \( A^T A \). The square roots of these eigenvalues are the singular values.
3. Determine the eigenvectors of \( A A^T \) to form \( U \), and the eigenvectors of \( A^T A \) to form \( V \).
4. Construct \( \Sigma \) by placing the singular values on its diagonal, filling the rest with zeros.

Example with a 2×2 Matrix

To understand how SVD works, let’s decompose a simple \( 2 \times 2 \) matrix:

Using Singular Value Decomposition, we express \( A \) as:

where:

• \( U \): contains the left singular vectors of \( A \)
• \( \Sigma \): a diagonal matrix with the singular values of \( A \)
• \( V^T \): the transpose of a matrix whose columns are the right singular vectors of \( A \)

The computed decomposition gives:

To verify this decomposition, we can reconstruct \( A \) by multiplying \( U \), \( \Sigma \), and \( V^T \):

This confirms that SVD accurately decomposes \( A \) into its constituent components.

Truncated SVD

To simplify large matrices, we can use truncated SVD, which keeps only the top \( k \) singular values and their corresponding vectors. The low-rank approximation is:

Truncated SVD is widely used in data compression, noise reduction, and dimensionality reduction.

Key Properties

• Uniqueness: Singular values are unique, though \( U \) and \( V \) may vary when singular values are repeated.
• Rank: The rank of \( A \) equals the number of non-zero singular values.
• Condition Number: The ratio \( \sigma_1 / \sigma_r \) measures the sensitivity of \( A \) to numerical errors.
• Frobenius Norm: \( \|A\|_F = \sqrt{\sum \sigma_i^2} \), representing the “size” of the matrix.

Implementation

This section demonstrates how to implement Singular Value Decomposition (SVD) in C++ using a custom matrix class. We’ll use the power iteration method to compute singular values and vectors, offering a simple approach for educational purposes.

Matrix Class Implementation

First, we need a basic Matrix class to handle matrix operations such as multiplication, transposition, and norm calculation. Here’s the code:

#include <vector>
#include <iostream>
#include <cmath>
#include <numeric>
#include <stdexcept>

class Matrix {
private:
    std::vector<std::vector<double>> data;
    size_t rows, cols;

public:
    // Constructor initializes matrix with zeros
    Matrix(size_t r, size_t c) : rows(r), cols(c) {
        data.resize(r, std::vector<double>(c, 0.0));
    }

    // Access operators for reading and modifying elements
    double& operator()(size_t i, size_t j) { return data[i][j]; }
    const double& operator()(size_t i, size_t j) const { return data[i][j]; }

    // Get dimensions
    size_t numRows() const { return rows; }
    size_t numCols() const { return cols; }

    // Matrix multiplication
    Matrix operator*(const Matrix& other) const {
        if (cols != other.rows) {
            throw std::invalid_argument("Matrix dimensions don't match for multiplication");
        }
        Matrix result(rows, other.cols);
        for (size_t i = 0; i < rows; ++i) {
            for (size_t j = 0; j < other.cols; ++j) {
                double sum = 0.0;
                for (size_t k = 0; k < cols; ++k) {
                    sum += data[i][k] * other.data[k][j];
                }
                result(i, j) = sum;
            }
        }
        return result;
    }

    // Transpose operation
    Matrix transpose() const {
        Matrix result(cols, rows);
        for (size_t i = 0; i < rows; ++i) {
            for (size_t j = 0; j < cols; ++j) {
                result(j, i) = data[i][j];
            }
        }
        return result;
    }

    // Frobenius norm (square root of sum of squared elements)
    double norm() const {
        double sum = 0.0;
        for (size_t i = 0; i < rows; ++i) {
            for (size_t j = 0; j < cols; ++j) {
                sum += data[i][j] * data[i][j];
            }
        }
        return std::sqrt(sum);
    }
};

Explanation of the Code

This class provides the essential matrix operations needed for SVD. Here’s what each method does:

• operator(): Allows accessing or modifying matrix elements using row and column indices.
• operator*: Implements matrix multiplication using nested loops. Throws an exception if the dimensions are incompatible.
• transpose(): Returns the transposed matrix by flipping rows and columns.
• norm(): Computes the Frobenius norm, which is the square root of the sum of squared elements.

SVD Implementation

Now, let’s define the SVD algorithm using the power iteration method, which repeatedly refines approximations for the largest singular values and their corresponding singular vectors.

struct SVDResult {
    Matrix U;
    Matrix S;
    Matrix V;

    SVDResult(size_t m, size_t n) :
        U(m, m), S(m, n), V(n, n) {}
};

class SVD {
private:
    // Power iteration to find the largest singular value and vectors
    static void powerIteration(const Matrix& A,
                             std::vector<double>& v,
                             double& sigma,
                             std::vector<double>& u,
                             size_t maxIter = 100,
                             double tol = 1e-10) {
        size_t m = A.numRows();
        size_t n = A.numCols();

        // Step 1: Initialize vector v with random values
        for (auto& val : v) {
            val = static_cast<double>(rand()) / RAND_MAX;
        }

        // Step 2: Normalize v
        double norm_v = std::sqrt(std::inner_product(v.begin(), v.end(), v.begin(), 0.0));
        for (auto& val : v) {
            val /= norm_v;
        }

        std::vector<double> prev_v = v;

        // Step 3: Iterate until convergence or maximum iterations
        for (size_t iter = 0; iter < maxIter; ++iter) {
            // Multiply A by v to get u
            u.assign(m, 0.0);
            for (size_t i = 0; i < m; ++i) {
                for (size_t j = 0; j < n; ++j) {
                    u[i] += A(i, j) * v[j];
                }
            }

            // Compute sigma (norm of u)
            sigma = std::sqrt(std::inner_product(u.begin(), u.end(), u.begin(), 0.0));

            // Normalize u
            for (auto& val : u) {
                val /= sigma;
            }

            // Multiply A^T by u to get new v
            v.assign(n, 0.0);
            for (size_t i = 0; i < n; ++i) {
                for (size_t j = 0; j < m; ++j) {
                    v[i] += A(j, i) * u[j];
                }
            }

            // Normalize v
            norm_v = std::sqrt(std::inner_product(v.begin(), v.end(), v.begin(), 0.0));
            for (auto& val : v) {
                val /= norm_v;
            }

            // Check convergence
            double diff = 0.0;
            for (size_t i = 0; i < n; ++i) {
                diff += std::abs(v[i] - prev_v[i]);
            }
            if (diff < tol) {
                break;
            }
            prev_v = v;
        }
    }

public:
    static SVDResult compute(const Matrix& A, size_t k) {
        size_t m = A.numRows();
        size_t n = A.numCols();
        SVDResult result(m, n);

        // Make a copy of A to work on
        Matrix B = A;

        // Step 4: Compute k singular values and vectors
        for (size_t i = 0; i < k; ++i) {
            std::vector<double> v(n);
            std::vector<double> u(m);
            double sigma;

            powerIteration(B, v, sigma, u);

            // Store results in U, S, and V matrices
            for (size_t j = 0; j < m; ++j) {
                result.U(j, i) = u[j];
            }
            result.S(i, i) = sigma;
            for (size_t j = 0; j < n; ++j) {
                result.V(j, i) = v[j];
            }

            // Deflate B by subtracting the outer product of u, v scaled by sigma
            for (size_t r = 0; r < m; ++r) {
                for (size_t c = 0; c < n; ++c) {
                    B(r, c) -= sigma * u[r] * v[c];
                }
            }
        }
        return result;
    }
};

Explanation of the Code

The compute method calculates the SVD by iteratively finding the largest singular values and their corresponding vectors using power iteration. Key points include:

• powerIteration: Computes one singular value and its corresponding singular vectors.
• deflation: Removes the influence of the current singular value to find subsequent ones.
• The loop runs k times to compute the top k singular values.

Usage Example

Here’s how to use the SVD implementation to decompose a test matrix:

int main() {
    // Define a 3x2 matrix
    Matrix A(3, 2);
    A(0,0) = 4.0;
    A(0,1) = 0.0;
    A(1,0) = 3.0; A(1,1) = -5.0;
    A(2,0) = 0.0; A(2,1) = 2.0;

    // Compute SVD, requesting 2 singular values
    SVDResult result = SVD::compute(A, 2);

    // Print singular values
    std::cout << "Singular Values:\n";
    for (size_t i = 0; i < 2; ++i) {
        std::cout << result.S(i, i) << " ";
    }
    std::cout << "\n";

    // Verify reconstruction
    Matrix reconstructed = result.U * result.S * result.V.transpose();
    double error = 0.0;
    for (size_t i = 0; i < A.numRows(); ++i) {
        for (size_t j = 0; j < A.numCols(); ++j) {
            error += std::abs(A(i, j) - reconstructed(i, j));
        }
    }
    std::cout << "Reconstruction error: " << error << "\n";

    return 0;
}

Explanation of the Code

This example shows how to use the SVD class to decompose a 3×2 matrix. Key steps:

• Create the matrix A and initialize it with values.
• Call the SVD::compute method to calculate the singular values and vectors. Here, k = 2 specifies the top 2 singular values to compute.
• Print the singular values from the diagonal of the S matrix.
• Reconstruct the original matrix using the formula \( A = U \Sigma V^T \) and calculate the reconstruction error.

Expected Output

Singular Values:
6.49097 3.44489
Reconstruction error: 1.22125e-15

Implementation with Eigen

The Eigen library provides powerful tools for performing Singular Value Decomposition (SVD). This section demonstrates how to use Eigen's SVD classes for practical purposes.

Overview of Eigen’s SVD Classes

Eigen offers three main SVD implementations, each optimized for specific use cases:

• JacobiSVD: Provides the most accurate results but is slower for large matrices.
• BDCSVD: Optimized for large matrices, using bidiagonalization for speed.
• CompleteOrthogonalDecomposition: Useful for rank-revealing decompositions.

Basic SVD with Eigen

Let’s perform SVD using the JacobiSVD class, which computes the singular values and singular vectors. The following code demonstrates this process:

#include <Eigen/Dense>
#include <iostream>

int main() {
    // Step 1: Create a 3x2 matrix
    Eigen::MatrixXd A(3, 2);
    A << 4, 0,
         3, -5,
         0, 2;

    std::cout << "Original matrix A:\n" << A << "\n\n";

    // Step 2: Compute the SVD using JacobiSVD
    const Eigen::JacobiSVD<Eigen::MatrixXd> svd(A,
        Eigen::ComputeFullU | Eigen::ComputeFullV);

    // Step 3: Get the singular values
    std::cout << "Singular values:\n" << svd.singularValues() << "\n\n";

    // Step 4: Extract U, S, and V matrices
    Eigen::MatrixXd U = svd.matrixU();
    Eigen::MatrixXd V = svd.matrixV();
    Eigen::MatrixXd S = Eigen::MatrixXd::Zero(A.rows(), A.cols());
    S.diagonal() = svd.singularValues();

    std::cout << "U matrix:\n" << U << "\n\n";
    std::cout << "S matrix:\n" << S << "\n\n";
    std::cout << "V matrix:\n" << V << "\n\n";

    // Step 5: Verify the decomposition
    Eigen::MatrixXd reconstructed = U * S * V.transpose();
    std::cout << "Reconstruction error: "
              << (A - reconstructed).norm() << "\n";

    return 0;
}

Explanation of the Code

The code above demonstrates how to use Eigen to perform SVD. Here’s a breakdown of the steps:

• Step 1: Create a matrix \(A\) with predefined values. This matrix serves as the input for the SVD.
• Step 2: Use JacobiSVD to compute the singular values and vectors of \(A\). The flags Eigen::ComputeFullU and Eigen::ComputeFullV ensure that the full \(U\) and \(V\) matrices are computed.
• Step 3: Extract the singular values using svd.singularValues(). These values are stored in a vector.
• Step 4: Retrieve the \(U\), \(S\), and \(V\) matrices. The \(S\) matrix is constructed as a diagonal matrix from the singular values.
• Step 5: Verify the decomposition by reconstructing \(A\) using \(U \Sigma V^T\). Calculate the reconstruction error, which should be near zero.

Expected Output

When the code is executed, it produces the following output:

Original matrix A:
 4  0
 3 -5
 0  2

Singular values:
6.49097
3.44489

U matrix:
-0.405933  0.873624 -0.268328
-0.884011 -0.300872  0.357771
 0.231825  0.382436  0.894427

S matrix:
6.49097       0
      0 3.44489
      0       0

V matrix:
-0.658725  0.752384
 0.752384  0.658725

Reconstruction error: 2.47508e-15

Explanation of the Output

The output consists of:

• Singular values: The diagonal elements of the \(S\) matrix, representing the magnitude of each dimension in the transformation.
• U matrix: Contains the left singular vectors as columns, representing the eigenvectors of \(A A^T\).
• S matrix: A diagonal matrix with singular values arranged in descending order.
• V matrix: Contains the right singular vectors as columns, representing the eigenvectors of \(A^T A\).
• Reconstruction error: The Frobenius norm of the difference between \(A\) and the reconstructed matrix \(U \Sigma V^T\), indicating numerical precision.

Setting Up Eigen with CMake

To compile your project with Eigen, use the following CMake configuration:

cmake_minimum_required(VERSION 3.29)
project(lu_decomposition)

set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNDEBUG -DEIGEN_NO_DEBUG")

# Path to the Eigen library
set(EIGEN3_INCLUDE_DIR "/Users/your_username/Downloads/eigen-3.4.0/")

add_executable(lu_decomposition main.cpp)

# Include directories
target_include_directories(${PROJECT_NAME} PRIVATE
        ${EIGEN3_INCLUDE_DIR}
)
        

Compilation Steps

Once the CMake configuration file (CMakeLists.txt) is set up, compile your project with the following commands:

1. Navigate to your project directory:

   cd /path/to/your/project

2. Create a build directory and navigate into it:

   mkdir build && cd build

3. Run CMake to generate the build files:

   cmake ..

4. Compile your project:

   cmake --build .

Performance Considerations

To enhance the performance of SVD for large matrices or real-time applications, consider the following tips:

Quick Example

The following example demonstrates how to efficiently compute the Singular Value Decomposition (SVD) of a matrix while optimizing performance using Eigen's BDCSVD class with thin computation.

#include <Eigen/Dense>

// Efficient SVD computation for large matrices
Eigen::MatrixXd efficientSVD(const Eigen::MatrixXd& matrix, int k) {
    // Use BDCSVD with thin computation for better performance
    Eigen::BDCSVD<Eigen::MatrixXd> svd(matrix,
        Eigen::ComputeThinU | Eigen::ComputeThinV);

    // Only compute top k components
    return svd.matrixU().leftCols(k) *
           svd.singularValues().head(k).asDiagonal() *
           svd.matrixV().leftCols(k).transpose();
}

Explanation of the Code

The function efficientSVD computes the SVD of the input matrix and reconstructs an approximation of the matrix using only the top k singular values and their corresponding singular vectors. Here's a breakdown:

• Eigen::BDCSVD<Eigen::MatrixXd> svd(matrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
  • This initializes the SVD computation using Eigen's BDCSVD class, which is optimized for large matrices through bidiagonalization. It supports efficient computation of singular values and vectors.
  • The flags Eigen::ComputeThinU and Eigen::ComputeThinV specify that the algorithm should compute only the first min(rows, cols) columns of the matrices U and V. This is ideal for performance and memory usage when the full matrices are not required.
• svd.singularValues().head(k).asDiagonal():
  • Extracts the top k singular values from the diagonal matrix \( \Sigma \).
• svd.matrixU().leftCols(k) and svd.matrixV().leftCols(k):
  • Select the corresponding singular vectors (columns of U and V) associated with the top k singular values.
• svd.matrixU().leftCols(k) * svd.singularValues().head(k).asDiagonal() * svd.matrixV().leftCols(k).transpose():
  • Reconstructs an approximation of the original matrix using the selected singular values and vectors, retaining only the most significant features.

What ComputeThinU and ComputeThinV Do

By default, SVD computes the full matrices U (\( m \times m \)) and V (\( n \times n \)), which can be computationally expensive and unnecessary for many applications. The options ComputeThinU and ComputeThinV optimize this process:

• ComputeThinU: Instead of computing the full matrix \( U \) (\( m \times m \)), this flag computes only the first min(m, n) columns of \( U \), reducing memory usage and computation time.
• ComputeThinV: Similarly, this flag computes only the first min(m, n) columns of \( V \), rather than the full matrix \( V \) (\( n \times n \)).

Using these options is particularly advantageous when you are interested only in the most significant singular values and vectors, as in the case of dimensionality reduction or low-rank approximation.

Advanced SVD Analysis

Let’s dive into some advanced features of Eigen’s SVD implementation, including rank determination, pseudo-inverse computation, truncated SVD, and variance analysis.

The following example showcases a complete implementation of SVD analysis using Eigen's BDCSVD, which is optimized for large matrices:

#include <Eigen/Dense>
#include <iostream>
#include <iomanip>
#include <limits>

void analyzeSVD(const Eigen::MatrixXd& A) {
    std::cout << "\nOriginal matrix A:\n" << A << "\n\n";

    // Using BDCSVD for better performance on large matrices
    Eigen::BDCSVD<Eigen::MatrixXd> svd(A,
        Eigen::ComputeFullU | Eigen::ComputeFullV);

    // Get rank using singular values
    double threshold = std::max(A.cols(), A.rows()) *
                      svd.singularValues().maxCoeff() *
                      std::numeric_limits<double>::epsilon();

    int rank = 0;
    for (int i = 0; i < svd.singularValues().size(); ++i) {
        if (svd.singularValues()(i) > threshold) {
            rank++;
        }
    }

    // Compute condition number
    double conditionNumber = svd.singularValues()(0) /
                            svd.singularValues()(svd.singularValues().size() - 1);

    std::cout << "Matrix Analysis:\n";
    std::cout << "Rank: " << rank << "\n";
    std::cout << "Condition number: " << conditionNumber << "\n";

    // Compute pseudo-inverse
    Eigen::MatrixXd pinv = svd.matrixV() *
        (svd.singularValues().array().inverse().matrix().asDiagonal()) *
        svd.matrixU().transpose();

    std::cout << "\nPseudo-inverse:\n" << pinv << "\n";

    // Truncated SVD (keep top k singular values)
    int k = 1; // Number of components to keep
    Eigen::MatrixXd truncated = svd.matrixU().leftCols(k) *
                               svd.singularValues().head(k).asDiagonal() *
                               svd.matrixV().leftCols(k).transpose();

    std::cout << "\nTruncated SVD (k=" << k << "):\n" << truncated << "\n";

    // Variance explained
    double totalVariance = svd.singularValues().array().square().sum();
    Eigen::ArrayXd explained = svd.singularValues().array().square() /
                              totalVariance * 100;

    std::cout << "\nVariance explained by each component:\n";
    for (int i = 0; i < explained.size(); ++i) {
        std::cout << "Component " << i + 1 << ": "
                 << std::fixed << std::setprecision(2)
                 << explained(i) << "%\n";
    }

    // Print singular values for verification
    std::cout << "\nSingular values:\n" << svd.singularValues() << "\n";
}

int main() {
    // Example 1: Well-conditioned matrix
    Eigen::MatrixXd A1(3, 3);
    A1 << 4, 0, 0,
          0, 3, 0,
          0, 0, 2;
    std::cout << "=== Test 1: Well-conditioned diagonal matrix ===";
    analyzeSVD(A1);

    // Example 2: Rank-deficient matrix
    Eigen::MatrixXd A2(3, 3);
    A2 << 1, 2, 2,
          2, 4, 4,
          3, 6, 6;
    std::cout << "\n=== Test 2: Rank-deficient matrix ===";
    analyzeSVD(A2);

    // Example 3: Rectangular matrix
    Eigen::MatrixXd A3(4, 2);
    A3 << 1, 2,
          3, 4,
          5, 6,
          7, 8;
    std::cout << "\n=== Test 3: Rectangular matrix ===";
    analyzeSVD(A3);

    return 0;
}

Explanation of the Code

The code above demonstrates advanced usage of SVD with Eigen. Here’s what each part does:

• Rank Calculation: The rank of a matrix is determined by counting singular values above a calculated threshold. The threshold is based on the size of the matrix, the largest singular value, and machine precision.
• Condition Number: This is the ratio of the largest singular value to the smallest, indicating the numerical stability of the matrix. A high condition number suggests a nearly singular matrix.
• Pseudo-inverse: The pseudo-inverse is computed using the formula \( V \Sigma^{-1} U^T \). It is especially useful for solving systems of linear equations when the matrix is not square or full-rank.
• Truncated SVD: A low-rank approximation of the matrix is generated by keeping only the top \( k \) singular values. This is useful for dimensionality reduction.
• Variance Explained: The proportion of variance explained by each component is computed as the square of singular values divided by the total variance.

Output and Analysis

Below are the complete outputs for each test case, including rank, pseudo-inverse, truncated SVD, variance explained, and singular values.

=== Test 1: Well-conditioned diagonal matrix ===
Original matrix A:
4 0 0
0 3 0
0 0 2

Matrix Analysis:
Rank: 3
Condition number: 2

Pseudo-inverse:
    0.25        0        0
       0 0.333333        0
       0        0      0.5

Truncated SVD (k=1):
4 0 0
0 0 0
0 0 0

Variance explained by each component:
Component 1: 55.17%
Component 2: 31.03%
Component 3: 13.79%

Singular values:
4.00
3.00
2.00

=== Test 2: Rank-deficient matrix ===
Original matrix A:
1 2 2
2 4 4
3 6 6

Matrix Analysis:
Rank: 1
Condition number: inf

Pseudo-inverse:
 inf -inf  nan
-inf  inf  nan
 nan  nan  nan

Truncated SVD (k=1):
1 2 2
2 4 4
3 6 6

Variance explained by each component:
Component 1: 100.00%
Component 2: 0.00%
Component 3: 0.00%

Singular values:
11.22
 0.00
 0.00

=== Test 3: Rectangular matrix ===
Original matrix A:
1 2
3 4
5 6
7 8

Matrix Analysis:
Rank: 2
Condition number: 22.76

Pseudo-inverse:
-1.00 -0.50 -0.00  0.50
 0.85  0.45  0.05 -0.35

Truncated SVD (k=1):
1.40 1.67
3.20 3.83
5.01 5.99
6.82 8.15

Variance explained by each component:
Component 1: 99.81%
Component 2: 0.19%

Singular values:
14.27
 0.63

Practical Applications

SVD has numerous applications across various domains. Here, we demonstrate a practical use case: image compression. By leveraging SVD, we can represent an image with reduced storage while retaining most of its essential features.

Image Compression with SVD

The following code compresses an image by retaining only the top \(k\) singular values, which capture the most significant information about the image:

#include <iostream>
#include <Eigen/Dense>

Eigen::MatrixXd compressImage(const Eigen::MatrixXd& image, int k) {
    Eigen::BDCSVD<Eigen::MatrixXd> svd(image,
        Eigen::ComputeThinU | Eigen::ComputeThinV);

    // Keep only k singular values
    return svd.matrixU().leftCols(k) *
           svd.singularValues().head(k).asDiagonal() *
           svd.matrixV().leftCols(k).transpose();
}

int main() {
    // Create a sample image matrix
    Eigen::MatrixXd image = Eigen::MatrixXd::Random(100, 100);

    // Compress using different numbers of components
    std::vector<int> components = {5, 10, 20};

    for(int k : components) {
        Eigen::MatrixXd compressed = compressImage(image, k);

        double error = (image - compressed).norm() / image.norm();
        std::cout << "Components: " << k
                 << ", Relative error: " << error
                 << ", Compression ratio: "
                 << (double)(100*100)/(k*(100+100)) << "\n";
    }

    return 0;
}

Explanation of the Code

The code above demonstrates how to use SVD for image compression:

• Input: The function compressImage takes an image (represented as a matrix) and a parameter \(k\), which specifies the number of singular values to retain.
• SVD Computation: The BDCSVD class decomposes the image matrix into its singular value components. The ComputeThinU and ComputeThinV options ensure efficient computation of reduced matrices when \(k \ll \text{min}(m, n)\).
• Truncation: The function reconstructs a compressed image by multiplying the truncated \(U\), \(\Sigma\), and \(V^T\) matrices.
• Output: For each specified \(k\), the relative error (normalized reconstruction error) and compression ratio are computed and displayed.

Output and Results

Below is the output for a sample run using a random \(100 \times 100\) image matrix:

Components: 5, Relative error: 0.910224, Compression ratio: 10
Components: 10, Relative error: 0.82944, Compression ratio: 5
Components: 20, Relative error: 0.684721, Compression ratio: 2.5

As seen in the output, increasing the number of components (\(k\)) improves the reconstruction quality because more singular values capture additional details about the original image matrix. Each singular value represents a specific level of significance in the data, with the largest singular values encoding the most critical patterns or features. By retaining more singular values, the compressed image retains more of these features, resulting in a more accurate reconstruction with a lower relative error. However, this also reduces the compression ratio because the storage requirements for the retained \(U\), \(\Sigma\), and \(V^T\) matrices increase. The balance between reconstruction quality and compression ratio highlights the trade-off inherent in image compression using SVD.

Applications of Singular Value Decomposition (SVD)

Singular Value Decomposition (SVD) is far more than a theoretical tool in linear algebra—it’s a practical technique that powers numerous real-world applications. From reducing complex datasets to compressing images and enabling personalized recommendations, SVD’s versatility makes it indispensable. Below, we’ll dive into some key areas where SVD plays a pivotal role.

1. Dimensionality Reduction

In the age of big data, reducing high-dimensional datasets to lower dimensions without losing critical information is essential. SVD is the backbone of techniques like Principal Component Analysis (PCA), helping extract the most significant features while retaining data variance. This is invaluable in tasks ranging from facial recognition to predictive analytics.

2. Image Compression

Storing high-resolution images can be costly in terms of storage and bandwidth. SVD simplifies this by compressing images while preserving visual quality. By keeping only the most significant singular values, SVD enables compact storage and efficient transmission, commonly used in image processing applications.

3. Recommendation Systems

Ever wondered how streaming platforms suggest content tailored to your taste? SVD is at the heart of collaborative filtering methods, where user-item interactions are decomposed into latent factors. This approach predicts user preferences and enhances recommendations, creating personalized experiences.

4. Signal Processing

In signal processing, SVD is a go-to tool for denoising and feature extraction. By isolating noise and reconstructing signals using dominant singular components, SVD improves clarity and retains essential characteristics, finding applications in audio processing, radar systems, and beyond.

5. Data Analysis and Statistics

SVD’s ability to compute matrix properties like rank, condition number, and explained variance makes it indispensable in statistical data analysis. It helps uncover patterns, relationships, and insights in complex datasets, empowering decision-making in fields like finance, healthcare, and social sciences.

Conclusion

Throughout this guide, we've explored Singular Value Decomposition from both theoretical and practical perspectives. We've covered custom implementation, Eigen library usage, and various real-world applications. Here are the key takeaways:

• SVD is a powerful technique for matrix factorization with applications in dimensionality reduction, data compression, and signal processing
• While custom implementations are valuable for learning, production code should leverage optimized libraries like Eigen
• The choice between different SVD algorithms depends on matrix size, required accuracy, and performance constraints
• Modern C++ features and the Eigen library make SVD implementation both efficient and maintainable

Congratulations on reading to the end of this comprehensive guide! If you found this guide valuable for your C++ journey, please consider citing or sharing it with fellow developers. Your support helps us continue creating comprehensive C++ resources for the development community.

Be sure to explore the Further Reading section for additional resources on linear algebra, matrix decompositions, and modern C++ practices.