pseudoinverse: Pseudoinverse of a Matrix

In corpcor: Efficient Estimation of Covariance and (Partial) Correlation

Description

The standard definition for the inverse of a matrix fails if the matrix is not square or singular. However, one can generalize the inverse using singular value decomposition. Any rectangular real matrix M can be decomposed as

M = U D V',

where U and V are orthogonal, V' means V transposed, and D is a diagonal matrix containing only the positive singular values (as determined by tol, see also fast.svd).

The pseudoinverse, also known as Moore-Penrose or generalized inverse is then obtained as

iM = V D^(-1) U' .

Usage

pseudoinverse(m, tol)

Arguments

m	matrix
tol	tolerance - singular values larger than tol are considered non-zero (default value: tol = max(dim(m))*max(D)*.Machine$double.eps)

Details

The pseudoinverse has the property that the sum of the squares of all the entries in iM %*% M - I, where I is an appropriate identity matrix, is minimized. For non-singular matrices the pseudoinverse is equivalent to the standard inverse.

Value

A matrix (the pseudoinverse of m).

Author(s)

Korbinian Strimmer (https://strimmerlab.github.io).

See Also

solve, fast.svd

Examples

# load corpcor library
library("corpcor")

# a singular matrix
m = rbind(
c(1,2),
c(1,2)
)

# not possible to invert exactly
try(solve(m))

# pseudoinverse
p = pseudoinverse(m)
p

# characteristics of the pseudoinverse
zapsmall( m %*% p %*% m )  ==  zapsmall( m )
zapsmall( p %*% m %*% p )  ==  zapsmall( p )
zapsmall( p %*% m )  ==  zapsmall( t(p %*% m ) )
zapsmall( m %*% p )  ==  zapsmall( t(m %*% p ) )


# example with an invertable matrix
m2 = rbind(
c(1,1),
c(1,0)
)
zapsmall( solve(m2) ) == zapsmall( pseudoinverse(m2) )

Example output

Error in solve.default(m) : 
  Lapack routine dgesv: system is exactly singular: U[2,2] = 0
     [,1] [,2]
[1,]  0.1  0.1
[2,]  0.2  0.2
     [,1] [,2]
[1,] TRUE TRUE
[2,] TRUE TRUE
     [,1] [,2]
[1,] TRUE TRUE
[2,] TRUE TRUE
     [,1] [,2]
[1,] TRUE TRUE
[2,] TRUE TRUE
     [,1] [,2]
[1,] TRUE TRUE
[2,] TRUE TRUE
     [,1] [,2]
[1,] TRUE TRUE
[2,] TRUE TRUE

corpcor documentation built on Sept. 16, 2021, 5:12 p.m.