Solve/Pseudo-Inverse with Singular Value Decomposition
DESCRIPTION:
Given the singular-value decomposition of a matrix, either solves (in the least-squares sense) a system of linear equations with that matrix as coefficient matrix or else computes the pseudo-inverse of the matrix.
USAGE: 
solve.svd.Matrix(a, b, tol=0, transpose=F)
REQUIRED ARGUMENTS:
a:
An object of class svd.Matrix, representing the singular-value decomposition of a matrix.
OPTIONAL ARGUMENTS:
b:
A matrix or vector. If transpose=T The number of rows of b must equal the number of rows of the matrix underlying a, while if transpose=F the number of rows of b must equal the number of columns of the matrix underlying a.
tol:
Tolerance for reciprocal condition number. Singular values are considered nonzero only if their ratio with the largest singular value exceeds tol. By default, tol = 0.
transpose:
A logical variable indicating whether or not the transpose (conjugate transpose if complex) of a is to be used for the solve or inverse operation. The default is to use the untransposed matrix.
VALUE:
If A is the matrix whose singular-value decomposition is represented by a, an object of class "Matrix" is returned that is the minimum l2 norm least-squares solution x to the system of equations A %*% x = b if b is present; otherwise the pseudo-inverse of A is returned. Attributes include a copy of the call to solve, the rank of the matrix assumed in order to obtain the solution, and the infinity norm reciprocal condition estimate if tol is nonnegative.
DETAILS:
Can be used for matrices that are rank-deficient, i.e., whose rank is less than their minimum dimension.
REFERENCES:
Golub, G. H., and C. F. Van Loan (1989), Matrix Computations, 2nd edition, Johns Hopkins, Baltimore.
SEE ALSO:
rcond.svd.Matrix , solve , solve.Matrix , solve.qr.Matrix .
EXAMPLES: 
m <- sample(1:9, 1); n <- sample(1:9, 1)
a <- Matrix( sample(-9:9, m*n, replace = T), nrow = m, ncol = n)
b <- rnorm(m)
z <- svd(a)                               # singular-value decomp of a
t(a) %*% (a %*% solve(a,b) - b)           # residual for normal equations
(solve(a) %*% b) - solve(a,b)