Exercise 4: Singular Value Decomposition

In this exercise we turn our attention to the Singular Value Decomposition (SVD). Remember that in real symmetric (or complex Hermitian) matrices, singular values coincide with eigenvalues, but in general this is not the case. The SVD is defined for any matrix, even rectangular ones. Singular values are always non-negative real values.

This example works also by reading a matrix from a file. In particular, the matrix to be used is related to a 2D reaction-diffusion model. More details about this problem can be found at Matrix Market.

Compiling

Copy the file ex14.c [plain text] to your directory and add these lines to the makefile

ex14: ex14.o
        -${CLINKER} -o ex14 ex14.o ${SLEPC_SVD_LIB}
        ${RM} ex14.o

Note: In the above text, the blank space in the 2nd and 3rd lines represents a tab.

Build the executable with the command

$ make ex14

Running the Program

In order to run the program, type the following

$ ./ex14 -file $SLEPC_DIR/share/slepc/datafiles/matrices/rdb200.petsc

You should get an output similar to this

Singular value problem stored in file.

 Reading REAL matrix from a binary file...
 Number of iterations of the method: 3
 Solution method: cross

 Number of requested singular values: 1
 Stopping condition: tol=1e-08, maxit=100
 SVD solve converged (1 singular triplet) due to CONVERGED_TOL; iterations 3
 ---------------------- --------------------
          sigma            relative error
 ---------------------- --------------------
        35.007519           1.40809e-10
 ---------------------- --------------------

Source Code Details

The way in which the SVD object works is very similar to that of EPS. However, some important differences exist. Examine the source code of the example program and pay attention to the differences with respect to EPS. After loading the matrix, the problem is solved by the following sequence of function calls:

SVDCreate(MPI_Comm comm,SVD *svd);
SVDSetOperators(SVD svd,Mat A,Mat B);
SVDSetFromOptions(SVD svd);
SVDSolve(SVD svd);
SVDGetConverged(SVD svd, int *nconv);
SVDGetSingularTriplet(SVD svd,int i,PetscReal *sigma,Vec u,Vec v);
SVDDestroy(SVD svd);

First, the singular value solver (SVD) context is created and the matrix associated with the problem is specified. Then various options are set for customized solution. After that, the program solves the problem, retrieves the solution, and finally destroys the SVD context.

Note that the singular value, sigma, is defined as a PetscReal, and that the singular vectors are simple Vec's.

SVD Options

Most of the options available in the EPS object have their equivalent in SVD. A full list of command-line options can be obtained by running the example with the option -help.

To show information about the SVD solver, add the -svd_view option:

$ ./ex14 -file $SLEPC_DIR/share/slepc/datafiles/matrices/rdb200.petsc -svd_view

Note: All the command-line options related to the SVD object have the -svd_ prefix.

Your output will include something like this

SVD Object: 1 MPI processes
  type: cross
    implicit matrix
    EPS Object: (svd_cross_) 1 MPI processes
      type: krylovschur
        50% of basis vectors kept after restart
        using the locking variant
      problem type: symmetric eigenvalue problem
      selected portion of the spectrum: largest real parts
      number of eigenvalues (nev): 1
      number of column vectors (ncv): 16
      maximum dimension of projected problem (mpd): 16
      maximum number of iterations: 100
      tolerance: 1e-09
      convergence test: relative to the eigenvalue
    BV Object: (svd_cross_) 1 MPI processes
      type: svec
      17 columns of global length 200
      vector orthogonalization method: classical Gram-Schmidt
      orthogonalization refinement: if needed (eta: 0.7071)
      block orthogonalization method: GS
      doing matmult as a single matrix-matrix product
    DS Object: (svd_cross_) 1 MPI processes
      type: hep
      solving the problem with: Implicit QR method (_steqr)
    ST Object: (svd_cross_) 1 MPI processes
      type: shift
      shift: 0.
      number of matrices: 1
      using a shell matrix
  transpose mode: explicit
  selected portion of the spectrum: largest
  number of singular values (nsv): 1
  number of column vectors (ncv): 16
  maximum dimension of projected problem (mpd): 16
  maximum number of iterations: 100
  tolerance: 1e-08
  convergence test: relative to the singular value

The output shows all the options that are susceptible of being changed, either from the command line or from the source code of the program: the method, the portion of the spectrum (largest or smallest singular values), the number of singular values (nsv), etc.

Try to change some of the values, for instance:

$ ./ex14 -file $SLEPC_DIR/share/slepc/datafiles/matrices/rdb200.petsc
         -svd_nsv 10 -svd_ncv 40 -svd_smallest

The "transpose mode" refers to whether the transpose of matrix A is being built explicitly or not (see SVDSetImplicitTranspose for an explanation).

Note that in the sample output above, the SVD object contains an EPS object. This only happens in some SVD solver types, as detailed below.

Changing the Singular Value Solver

SLEPc provides several solvers for computing the SVD, which can be selected in the source code with the function SVDSetType, or at run time:

$ ./ex14 -file $SLEPC_DIR/share/slepc/datafiles/matrices/rdb200.petsc
         -svd_type trlanczos

The following table shows the list of SVD solvers available in SLEPc.

Solver	Command-line Name	Parameter
Cross Product	cross	SVDCROSS
Cyclic Matrix	cyclic	SVDCYCLIC
Lanczos with Explicit Restart	lanczos	SVDLANCZOS
Lanczos with Thick Restart	trlanczos	SVDTRLANCZOS
Lapack	lapack	SVDLAPACK

Note: The Lapack solver is not really a full-featured singular value solver but simply an interface to some LAPACK routines. These routines operate sequentially in dense mode and therefore are suitable only for small size problems. This solver should be used only for debugging purposes.

Note: The default solver is cross.

The first two solvers, cross and cyclic, are not real methods implemented in the SVD module, but are two convenient ways of solving the SVD problem by making use of the eigensolvers available in the EPS module. In those two cases, the SVD object manages an EPS object internally, whose parameters can be set as desired (typically only the method). For example:

$ ./ex14 -file $SLEPC_DIR/share/slepc/datafiles/matrices/rdb200.petsc
         -svd_type cyclic -svd_cyclic_eps_type lanczos -svd_cyclic_eps_lanczos_reorthog local