Wednesday, July 29, 2009

Compute SVD and Moore-Penrose pseudoinverse in SAS

Let the SVD of X be:


use the following SAS code:

ods output EigenValues=E;
ods output EigenVectors=VT;

PROC PRINCOMP NOINT COV data=X noprint;
var_1....var_k;
RUN;

in order to get the diagonal elements of singular value matrix S of SVD output, use the following formula:



\int_{0}^{1}\frac{x^{4}\left(1-x\right)^{4}}{1+x^{2}}dx
=\frac{22}{7}-\pi

where
E_{ii} 
is the diagonal elements of EigenValue matrix E and N is the number of observations in the data set X;

The U matrix can be calculated as:


To get the Moore-Penrose pseudo-inverse, simply apply:


where $X^{+}$ is the pseudoinverse of X and E^{+} is the pseudoinverse of E by taking the reciprocal of each non-zero element on the diagonal, and leaving the zeros in place. In numerical computation, only elements larger than some small tolerance are taken to be nonzero, and the others are replaced by zeros. In MATLAB, the tolerance is taken to be t = ε•max(m,n)•max(Σ), where ε is the machine epsilon. Here we can apply the same rule, using CONSTANT('MACEPS') function to get machine epsilon in SAS.

The result following this method is numerically solid.

A SAS macro for SVD computation is provided @ Here





No comments: