Example 10: Singular Value Decomposition (SVD)

This example demonstrates computing the SVD of a matrix \(A = U \Sigma V^H\).

Key Concepts

  1. Singular Values Only: Computing just \(\Sigma\) using svd_vals.

  2. Full SVD: Computing \(\Sigma\), \(U\), and \(V^H\) using svd (gesvd).

  3. Partial Vectors: Computing only \(U\) or only \(V^H\).

C++ Example

Setup (Lines 29-30)

slate::Matrix<scalar_type> A( m, n, nb, ... );
std::vector<real_t> Sigma( min_mn );

  • A: Input matrix m by n.

  • Sigma: Vector to store the singular values. Size must be at least min(m, n).

Singular Values Only (Lines 41-46)

slate::svd_vals( A, Sigma );
// or
slate::svd( A, Sigma );

If you only need the singular values (not the vectors U and V), use svd_vals or svd without matrix arguments. This is significantly faster than computing vectors.

Singular Vectors (Lines 57-70)

slate::Matrix<scalar_type>  U( m, min_mn, nb, ... );
slate::Matrix<scalar_type> VH( min_mn, n, nb, ... );

slate::svd( A, Sigma, U, VH );

To compute vectors:

  • U: Left singular vectors. Dimensions m by min(m, n).

  • VH: Right singular vectors (transposed). Dimensions min(m, n) by n.

  • Note: This example computes the reduced SVD.

Partial Vectors (Lines 75-80)

slate::Matrix<scalar_type> Uempty, Vempty;
slate::svd( A, Sigma, U, Vempty );  // only U
slate::svd( A, Sigma, Uempty, VH ); // only V^H

You can compute just U or just VH by passing an empty matrix placeholder for the unwanted component. This saves computation time.

  1// ex10_svd.cc
  2// Solve Singular Value Decomposition, A = U Sigma V^H
  3
  4/// !!!   Lines between `//---------- begin label`          !!!
  5/// !!!             and `//---------- end label`            !!!
  6/// !!!   are included in the SLATE Users' Guide.           !!!
  7
  8#include <slate/slate.hh>
  9
 10#include "util.hh"
 11
 12int mpi_size = 0;
 13int mpi_rank = 0;
 14int grid_p = 0;
 15int grid_q = 0;
 16
 17//------------------------------------------------------------------------------
 18template <typename scalar_type>
 19void test_svd()
 20{
 21    using real_t = blas::real_type<scalar_type>;
 22
 23    print_func( mpi_rank );
 24
 25    int64_t m=2000, n=1000, nb=256;
 26
 27    //---------- begin svd1
 28    int64_t min_mn = std::min( m, n );
 29    slate::Matrix<scalar_type> A( m, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
 30    std::vector<real_t> Sigma( min_mn );
 31    // ...
 32    //---------- end svd1
 33
 34    A.insertLocalTiles();
 35    random_matrix( A );
 36
 37    //----------------------------------------
 38    //---------- begin svd2
 39
 40    // A = U Sigma V^H, singular values only
 41    slate::svd_vals( A, Sigma );
 42    //---------- end svd2
 43    random_matrix( A );
 44
 45    //---------- begin svd3
 46    slate::svd( A, Sigma );
 47    //---------- end svd3
 48
 49    // todo: full SVD not yet supported?
 50
 51    //----------------------------------------
 52    //---------- begin svd4
 53
 54    // Singular vectors
 55    // U is m x min_mn (reduced SVD) or m x m (full SVD)
 56    // V is min_mn x n (reduced SVD) or n x n (full SVD)
 57    slate::Matrix<scalar_type>  U( m, min_mn, nb, grid_p, grid_q, MPI_COMM_WORLD );
 58    slate::Matrix<scalar_type> VH( min_mn, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
 59    // empty, 0-by-0 matrices as placeholders for U and VH.
 60    slate::Matrix<scalar_type> Uempty, Vempty;
 61    // ...
 62    //---------- end svd4
 63
 64    U.insertLocalTiles();
 65    VH.insertLocalTiles();
 66    random_matrix( A );
 67
 68    //---------- begin svd5
 69
 70    slate::svd( A, Sigma, U, VH );      // both U and V^H
 71    //---------- end svd5
 72    random_matrix( A );
 73
 74    //---------- begin svd6
 75    slate::svd( A, Sigma, U, Vempty );  // only U
 76    //---------- end svd6
 77    random_matrix( A );
 78
 79    //---------- begin svd7
 80    slate::svd( A, Sigma, Uempty, VH ); // only V^H
 81    //---------- end svd7
 82}
 83
 84//------------------------------------------------------------------------------
 85int main( int argc, char** argv )
 86{
 87    try {
 88        // Parse command line to set types for s, d, c, z precisions.
 89        bool types[ 4 ];
 90        parse_args( argc, argv, types );
 91
 92        int provided = 0;
 93        slate_mpi_call(
 94            MPI_Init_thread( &argc, &argv, MPI_THREAD_MULTIPLE, &provided ) );
 95        assert( provided == MPI_THREAD_MULTIPLE );
 96
 97        slate_mpi_call(
 98            MPI_Comm_size( MPI_COMM_WORLD, &mpi_size ) );
 99
100        slate_mpi_call(
101            MPI_Comm_rank( MPI_COMM_WORLD, &mpi_rank ) );
102
103        // Determine p-by-q grid for this MPI size.
104        grid_size( mpi_size, &grid_p, &grid_q );
105        if (mpi_rank == 0) {
106            printf( "mpi_size %d, grid_p %d, grid_q %d\n",
107                    mpi_size, grid_p, grid_q );
108        }
109
110        // so random_matrix is different on different ranks.
111        srand( 100 * mpi_rank );
112
113        if (types[ 0 ]) {
114            test_svd< float >();
115        }
116
117        if (types[ 1 ]) {
118            test_svd< double >();
119        }
120
121        if (types[ 2 ]) {
122            test_svd< std::complex<float> >();
123        }
124
125        if (types[ 3 ]) {
126            test_svd< std::complex<double> >();
127        }
128
129        slate_mpi_call(
130            MPI_Finalize() );
131    }
132    catch (std::exception const& ex) {
133        fprintf( stderr, "%s", ex.what() );
134        return 1;
135    }
136    return 0;
137}