NAME

svd.h - singular value decomposition and friends

SYNOPSIS

Included here is a basic SVD call and some basic applications of the SVD, such as pseudo matrix inversion and principal component analysis.

#include <svd.h>

void svd(double *U,
  double *S, double *V, unsigned nrow,
  unsigned ncol);

void pinv(double **A, double **Ainv, unsigned n);

void spinv(double **A, double **Ainv, unsigned n);

int pca(DATASET *data,
  unsigned m, double ***D, double **S,
  double **M);

DESCRIPTION

The heart of this module is the SVD routine. The supplied SVD function is very small and simple but is by no means optimal. If this module is compiled with either a -DSL_EXTENSIONS or a -DLAPACK command line option and linked with either "-llapack -lblas -lf2c" or "-lsle", then the LAPACK SVD will be used which is vastly superior to the enclosed version. Also, at runtime, one can switch back to the original SVD routine by setting svd_original to a nonzero value. The proper way to declare svd_original is:

extern int svd_original;

Function Definitions

The following function prototypes are given in the header file svd.h.

void svd(double *U,
  double *S, double *V, unsigned nrow,
  unsigned ncol);

Computes the singular value decomposition of the matrix, U. The contents of U are replaced such that A = U*S*V' where A represents the initial value of U. S is understood to have only room for ncol elements. The matrix V may be NULL, in which case, no values are returned for V.

void pinv(double **A, double **Ainv, unsigned n);

Matrix pseudo-inversion, based on a singular value decomposition. The two matrices may point to the same memory.

void spinv(double **A, double **Ainv, unsigned n);

Matrix pseudo-inversion for symmetric positive definite matrices based on a singular value decomposition. The two matrices may point to the same memory.

int pca(DATASET *data,
  unsigned m, double ***D, double **S,
  double **M);

Computes the m principal components of the x portion of the data contained in data. The last three arguments are pointers to arrays of doubles that will be filled with useful values. You do not need to allocate space for these arrays, as it is done for you. After the call, D will point to an (m x n) matrix that contains the m largest principle components, sorted from largest to smallest, S will point to a vector that contains the associated eigenvalues, and M will point to a vector with the mean of the data (which you will need because PCA is properly defined on zero mean data.) The return value is the number of nonzero components actually returned, so if all goes well, it should be equal to m.

It is an error for m to be larger than the dimensionality of the input data. The space returned in D, S, and M must be deallocated with deallocate_array().

BUGS

The original SVD has a tolerance defined in the source code file svd.c. You may wish to change this.

AUTHOR

Gary William Flake (gary.flake@usa.net).