cisst-saw
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes
nmrSVD.h File Reference

Declaration of nmrSVD. More...

#include <cisstCommon/cmnThrow.h>
#include <cisstVector/vctFixedSizeMatrix.h>
#include <cisstVector/vctDynamicMatrix.h>
#include <cisstNumerical/nmrNetlib.h>
#include <cisstNumerical/nmrExport.h>

Go to the source code of this file.

Classes

class  nmrSVDDynamicData
 Data for SVD problem (Dynamic). More...
 
class  nmrSVDDynamicData::Friend
 
class  nmrSVDFixedSizeData< _rows, _cols, _storageOrder >
 Data of SVD problem (Fixed size). More...
 
class  nmrSVDFixedSizeData< _rows, _cols, _storageOrder >::Friend
 

Functions

Algorithm SVD: Singular Value Decomposition

The nmrSVD functions compute the singular value decomposition (SVD) of a real $ M \times N $ matrix A, optionally computing the left and/or right singular vectors. The SVD problem is written:

$ A = U * \Sigma * V^{T} $

where $ \Sigma $ is a $ M \times N $ matrix which is zero except for its min(m,n) diagonal elements, U is a $ M \times M $ orthogonal matrix, and V is a $ N \times N $ orthogonal matrix. The diagonal elements of $ \Sigma $ are the singular values of A; they are real and non-negative, and are returned in descending order. The first $ \mbox{min}(m,n) $ columns of U and V are the left and right singular vectors of A.

These functions are wrappers around the LAPACK routine dgesvd, therefore they share some features with the LAPACK routine:

  1. They return $ V^{T} $, not $ V $.
  2. On exit, the content of A is altered.
  3. The vectors and matrices must be compact, i.e. use a contiguous block of memory.

The nmrSVD functions add the following features:

  1. A simplified interface to the cisstVector matrices, either vctDynamicMatrix or vctFixedSizeMatrix.
  2. The possibility to use any storage order, i.e. VCT_ROW_MAJOR or VCT_COL_MAJOR. In the specific case of SVD, the storage order can be treated as a transpose and the properties of the decomposition are preserved (nmrSVD will swap U and Vt pointers). This is performed at no extra cost, i.e. no copy back and forth of the elements nor extra memory allocation.
  3. Input validation checks are performed, i.e. an std::runtime_error exception will be thrown if the sizes or storage order don't match or if the containers are not compact.
  4. Helper classes to allocate memory for the output and workspace: nmrSVDFixedSizeData and nmrSVDDynamicData.

There are different ways to call this function to compute the SVD of the matrix A. These correspond to different overloaded nmrSVD functions:

  1. Using a preallocated data object.

    The user creates the input matrix A:

    vctDynamicMatrix<double> A(12, 24 , VCT_ROW_MAJOR); // 12 x 24 matrix
    vctRandom(A, -10.0, 10.0);

    The user allocates a data object which could be of type nmrSVDFixedSizeData or nmrSVDDynamicData. corresponding to fixed size or dynamic matrix A:

    nmrSVDDynamicData data(A);

    Call the nmrSVD function:

    nmrSVD(A, data);

    The content of input matrix A is modified by this routine. The matrices U, Vt and vector S are available through the following methods

    std::cout << "U: " << data.U() << std::endl
    << "S: " << data.S() << std::endl
    << "V: " << data.Vt().Transpose() << std::endl;

  2. The user provides the matrices U, Vt, and vector S.

    The User allocates memory for these matrices and vector:

    vctDynamicMatrix<double> A(5, 4);
    vctRandom(A, -10.0, 10.0);
    vctDynamicMatrix<double> U(5, 5);
    vctDynamicMatrix<double> Vt(4, 4);
    vctDynamicVector<double> S(4);

    Call the SVD routine:

    nmrSVD(A, U, S, Vt);

    The SVD function verifies that the size of the data objects matches the input, and allocates workspace memory, which is deallocated when the function ends.

  3. The user provides the matrices U, Vt, and vector S along with workspace required by SVD routine of LAPACK.

    Create matrices and vector:

    vctDynamicMatrix<double> A(5, 4);
    vctRandom(A, -10.0, 10.0);
    vctDynamicMatrix<double> U(5, 5);
    vctDynamicMatrix<double> Vt(4, 4);
    vctDynamicVector<double> S(4);

    One also needs to allocate memory the for workspace. This approach is particularly useful when the user is using more than one numerical method from the library and is willing or need to share the workspace between them. In such as case, the user can allocate a block of memory greater than the minimum required by different methods. To help the user determine the minimum workspace required the library provides the helper function nmrSVDDynamicData::WorkspaceSize().

    vctDynamicVector<double> workspace(nmrSVDDynamicData::WorkspaceSize(A));

    Call the SVD function:

    nmrSVD(A, U, S, Vt, workspace);

    For fixed size containers, the example above would be:

    vctFixedSizeMatrix<double, 5, 4> A;
    vctRandom(A, -10.0, 10.0);
    vctFixedSizeMatrix<double, 5, 5> U;
    vctFixedSizeMatrix<double, 4, 4> Vt;
    vctFixedVector<double, 4> S;
    nmrSVDFixedSizeData<4, 3, VCT_COL_MAJOR>::VectorTypeWorkspace Workspace;
    nmrSVD(A, U, S, Vt, Workspace);

The functions nmrSVD return the status code as defined by LAPACK, i.e.:

  1. 0: Successful exit.

  2. Greater than 0: If SBDSQR did not converge, return value specifies how many superdiagonals of an intermediate bidiagonal form B did not converge to zero. See the description of workspace above for details.

    In this case Workspace(2:MIN(M,N)) contains the unconverged superdiagonal elements of an upper bidiagonal matrix $B$ whose diagonal is in $S$ (not necessarily sorted). $B$ satisfies $A = U * B * V^{T}$, so it has the same singular values as $A$, and singular vectors related by $U$ and $V^T$.

Note
The SVD functions make use of LAPACK routines. To activate this code, set the CISST_HAS_CISSTNETLIB flag to ON during the configuration of cisst with CMake.
The general rule for numerical functions which depend on LAPACK is that column major matrices should be used everywhere, and that all matrices should be compact. In this case, both row major and column major are allowed but they must not be mixed in a data object.
template<class _matrixOwnerType >
CISSTNETLIB_INTEGER nmrSVD (vctDynamicMatrixBase< _matrixOwnerType, CISSTNETLIB_DOUBLE > &A, nmrSVDDynamicData &data) throw (std::runtime_error)
 
template<class _matrixOwnerTypeA , class _matrixOwnerTypeU , class _vectorOwnerTypeS , class _matrixOwnerTypeVt , class _vectorOwnerTypeWorkspace >
CISSTNETLIB_INTEGER nmrSVD (vctDynamicMatrixBase< _matrixOwnerTypeA, CISSTNETLIB_DOUBLE > &A, vctDynamicMatrixBase< _matrixOwnerTypeU, CISSTNETLIB_DOUBLE > &U, vctDynamicVectorBase< _vectorOwnerTypeS, CISSTNETLIB_DOUBLE > &S, vctDynamicMatrixBase< _matrixOwnerTypeVt, CISSTNETLIB_DOUBLE > &Vt, vctDynamicVectorBase< _vectorOwnerTypeWorkspace, CISSTNETLIB_DOUBLE > &Workspace)
 
template<class _matrixOwnerTypeA , class _matrixOwnerTypeU , class _vectorOwnerTypeS , class _matrixOwnerTypeVt >
CISSTNETLIB_INTEGER nmrSVD (vctDynamicMatrixBase< _matrixOwnerTypeA, CISSTNETLIB_DOUBLE > &A, vctDynamicMatrixBase< _matrixOwnerTypeU, CISSTNETLIB_DOUBLE > &U, vctDynamicVectorBase< _vectorOwnerTypeS, CISSTNETLIB_DOUBLE > &S, vctDynamicMatrixBase< _matrixOwnerTypeVt, CISSTNETLIB_DOUBLE > &Vt)
 
template<vct::size_type _rows, vct::size_type _cols, vct::size_type _minmn, vct::size_type _work, bool _storageOrder>
CISSTNETLIB_INTEGER nmrSVD (vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _rows, _cols, _storageOrder > &A, vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _rows, _rows, _storageOrder > &U, vctFixedSizeVector< CISSTNETLIB_DOUBLE, _minmn > &S, vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _cols, _cols, _storageOrder > &Vt, vctFixedSizeVector< CISSTNETLIB_DOUBLE, _work > &workspace)
 
template<vct::size_type _rows, vct::size_type _cols, vct::size_type _minmn, bool _storageOrder>
CISSTNETLIB_INTEGER nmrSVD (vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _rows, _cols, _storageOrder > &A, vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _rows, _rows, _storageOrder > &U, vctFixedSizeVector< CISSTNETLIB_DOUBLE, _minmn > &S, vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _cols, _cols, _storageOrder > &Vt)
 
template<vct::size_type _rows, vct::size_type _cols, bool _storageOrder>
CISSTNETLIB_INTEGER nmrSVD (vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _rows, _cols, _storageOrder > &A, nmrSVDFixedSizeData< _rows, _cols, _storageOrder > &data)
 

Detailed Description

Declaration of nmrSVD.

Function Documentation

template<class _matrixOwnerType >
CISSTNETLIB_INTEGER nmrSVD ( vctDynamicMatrixBase< _matrixOwnerType, CISSTNETLIB_DOUBLE > &  A,
nmrSVDDynamicData &  data 
)
throw (std::runtime_error
)
inline

This function solves the SVD problem for a dynamic matrix using an nmrSVDDynamicData.

This function checks for valid input (size, storage order and compact) and calls the LAPACK function. If the input doesn't match the data, an exception is thrown (std::runtime_error).

This function modifies the input matrix A and stores the results in the data. Each component of the result can be obtained via the const methods nmrSVDDynamicData::U(), nmrSVDDynamicData::S() and nmrSVDDynamicData::Vt().

Parameters
AA matrix of size MxN, either vctDynamicMatrix or vctDynamicMatrixRef.
dataA data object corresponding to the input matrix.
Test:
nmrSVDTest::TestDynamicColumnMajor nmrSVDTest::TestDynamicRowMajor nmrSVDTest::TestDynamicColumnMajorUserAlloc nmrSVDTest::TestDynamicRowMajorUserAlloc
template<class _matrixOwnerTypeA , class _matrixOwnerTypeU , class _vectorOwnerTypeS , class _matrixOwnerTypeVt , class _vectorOwnerTypeWorkspace >
CISSTNETLIB_INTEGER nmrSVD ( vctDynamicMatrixBase< _matrixOwnerTypeA, CISSTNETLIB_DOUBLE > &  A,
vctDynamicMatrixBase< _matrixOwnerTypeU, CISSTNETLIB_DOUBLE > &  U,
vctDynamicVectorBase< _vectorOwnerTypeS, CISSTNETLIB_DOUBLE > &  S,
vctDynamicMatrixBase< _matrixOwnerTypeVt, CISSTNETLIB_DOUBLE > &  Vt,
vctDynamicVectorBase< _vectorOwnerTypeWorkspace, CISSTNETLIB_DOUBLE > &  Workspace 
)
inline

This function solves the SVD problem for a dynamic matrix using the storage provided by the user for both the output (U, S, Vt) and the workspace.

Internally, a data is created using the storage provided by the user (see nmrSVDDynamicData::SetRef). While the data is being build, the consistency of the output and workspace is checked. Then, the nmrSVD(A, data) function can be used safely.

Parameters
Ais a reference to a dynamic matrix of size MxN
U,S,VtThe output matrices and vector for SVD
WorkspaceThe workspace for LAPACK.
Test:
nmrSVDTest::TestDynamicColumnMajorUserAlloc nmrSVDTest::TestDynamicRowMajorUserAlloc
template<class _matrixOwnerTypeA , class _matrixOwnerTypeU , class _vectorOwnerTypeS , class _matrixOwnerTypeVt >
CISSTNETLIB_INTEGER nmrSVD ( vctDynamicMatrixBase< _matrixOwnerTypeA, CISSTNETLIB_DOUBLE > &  A,
vctDynamicMatrixBase< _matrixOwnerTypeU, CISSTNETLIB_DOUBLE > &  U,
vctDynamicVectorBase< _vectorOwnerTypeS, CISSTNETLIB_DOUBLE > &  S,
vctDynamicMatrixBase< _matrixOwnerTypeVt, CISSTNETLIB_DOUBLE > &  Vt 
)
inline

This function solves the SVD problem for a dynamic matrix using the storage provided by the user for the output (U, S, Vt). A workspace will be dynamically allocated.

Internally, a data is created using the storage provided by the user (see nmrSVDDynamicData::SetRefOutput). While the data is being build, the consistency of the output is checked. Then, the nmrSVD(A, data) function can be used safely.

Parameters
Ais a reference to a dynamic matrix of size MxN
U,S,VtThe output matrices and vector for SVD
Test:
nmrSVDTest::TestDynamicColumnMajorUserAlloc nmrSVDTest::TestDynamicRowMajorUserAlloc
Warning
Again, this method will dynamically allocate a workspace at each call. This is somewhat inefficient if the method is to be called many times. For a real-time task, the dynamic allocation might even break your application.
template<vct::size_type _rows, vct::size_type _cols, vct::size_type _minmn, vct::size_type _work, bool _storageOrder>
CISSTNETLIB_INTEGER nmrSVD ( vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _rows, _cols, _storageOrder > &  A,
vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _rows, _rows, _storageOrder > &  U,
vctFixedSizeVector< CISSTNETLIB_DOUBLE, _minmn > &  S,
vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _cols, _cols, _storageOrder > &  Vt,
vctFixedSizeVector< CISSTNETLIB_DOUBLE, _work > &  workspace 
)
inline

This function solves the SVD problem for a fixed size matrix using the storage provided by the user for both the output (U, S, Vt) and the workspace.

The storage order of the matrices and their sizes must match at compilation time. This is enforced by the template parameters and matching problems will lead to compilation errors. Since there is no easy way to enforece the size of the vectors S and workspace with template parameters, a runtime check is performed. The test uses CMN_ASSERT to determine what to do if the sizes don't match. By default CMN_ASSERT calls abort() but it can be configured to be ignored or to throw an exception (see CMN_ASSERT for details).

This function modifies the input matrix A and the workspace. It stores the results in U, S and Vt.

Parameters
Ais a reference to a dynamic matrix of size MxN
U,S,VtThe output matrices and vector for SVD
workspaceThe workspace for LAPACK.
Test:
nmrSVDTest::TestFixedSizeColumnMajorMLeqN_T2 nmrSVDTest::TestFixedSizeRowMajorMLeqN_T2 nmrSVDTest::TestFixedSizeColumnMajorMGeqN_T2 nmrSVDTest::TestFixedSizeRowMajorMGeqN_T2
template<vct::size_type _rows, vct::size_type _cols, vct::size_type _minmn, bool _storageOrder>
CISSTNETLIB_INTEGER nmrSVD ( vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _rows, _cols, _storageOrder > &  A,
vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _rows, _rows, _storageOrder > &  U,
vctFixedSizeVector< CISSTNETLIB_DOUBLE, _minmn > &  S,
vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _cols, _cols, _storageOrder > &  Vt 
)
inline

This function solves the SVD problem for a fixed size matrix using the storage provided by the user for the output (U, S, Vt). The workspace will be automatically allocated (stack allocation of a fixed size vector).

This method calls nmrSVD(A, U, S, Vt, workspace).

Parameters
Ais a reference to a dynamic matrix of size MxN
U,S,VtThe output matrices and vector for SVD
Test:
nmrSVDTest::TestFixedSizeColumnMajorMLeqN_T2 nmrSVDTest::TestFixedSizeRowMajorMLeqN_T2 nmrSVDTest::TestFixedSizeColumnMajorMGeqN_T2 nmrSVDTest::TestFixedSizeRowMajorMGeqN_T2
template<vct::size_type _rows, vct::size_type _cols, bool _storageOrder>
CISSTNETLIB_INTEGER nmrSVD ( vctFixedSizeMatrix< CISSTNETLIB_DOUBLE, _rows, _cols, _storageOrder > &  A,
nmrSVDFixedSizeData< _rows, _cols, _storageOrder > &  data 
)
inline

This function solves the SVD problem for a fixed size matrix using the storage provided by the user for both the output (U, S, Vt) and the workspace. This function allows to use the very convenient nmrSVDFixedSizeData to allocate the memory required for U, S, Vt and the workspace:

vctFixedSizeMatrix<double, 12, 7> A;
vctRandom(A, -10.0, 10.0);
nmrSVDFixedSizeData<12, 7> data;
nmrSVD(A, data);
std::cout << "U: " << data.U() << std::endl
<< "S: " << data.S() << std::endl
<< "V: " << data.Vt().Transpose() << std::endl;

This method calls nmrSVD(A, U, S, Vt, workspace).

Parameters
AA fixed size matrix of size MxN.
dataA data object.
Test:
nmrSVDTest::TestFixedSizeColumnMajorMLeqN nmrSVDTest::TestFixedSizeRowMajorMLeqN nmrSVDTest::TestFixedSizeColumnMajorMGeqN nmrSVDTest::TestFixedSizeRowMajorMGeqN
Generated by   doxygen 1.8.6