Calculate inverse of a non-square matrix in R

Question

I'm pretty new to the R language and trying to find out how you can calculate the inverse of a matrix that isn't square. (non-square? irregular? I'm unsure of the correct terminology).

From my book and a quick google search, (see source), I've found you can use solve(a) to find the inverse of a matrix if a is square.

The matrix I have created is, and from what I understand, not square:

  > matX<-matrix(c(rep(1, 8),2,3,4,0,6,4,3,7,-2,-4,3,5,7,8,9,11), nrow=8, ncol=3);
  > matX

       [,1] [,2] [,3]
  [1,]    1    2   -2
  [2,]    1    3   -4
  [3,]    1    4    3
  [4,]    1    0    5
  [5,]    1    6    7
  [6,]    1    4    8
  [7,]    1    3    9
  [8,]    1    7   11
  >

Is there a function for solving a matrix of this size or will I have to do something to each element? as the solve() function gives this error:

  Error in solve.default(matX) : 'a' (8 x 3) must be square

The calculation I'm trying to achieve from the above matrix is: (matX'matX)^-1

Thanks in advance.

Referring to your last line, try: `solve(t(matX) %*% matX)` – Mark Heckmann Jan 26 '14 at 13:53 — Mark Heckmann, Jan 26 '14 at 13:53

G. Grothendieck · Accepted Answer · 2019-11-15T14:18:47.157

ginv ginv in the MASS package will give the generalized inverse of a matrix. Premultiplying the original matrix by it will give the identity:

library(MASS)
inv <- ginv(matX)

# test it out
inv %*% matX
##               [,1]         [,2]          [,3]
## [1,]  1.000000e+00 6.661338e-16  4.440892e-15
## [2,] -8.326673e-17 1.000000e+00 -1.110223e-15
## [3,]  6.938894e-17 8.326673e-17  1.000000e+00

As suggested in the comments this can be displayed in a nicer way using zapsmall:

zapsmall(inv %*% matX)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

The inverse of matX'matX is now:

tcrossprod(inv)
##              [,1]         [,2]         [,3]
## [1,]  0.513763935 -0.104219636 -0.002371406
## [2,] -0.104219636  0.038700372 -0.007798748
## [3,] -0.002371406 -0.007798748  0.006625269

solve however, if you aim is to calculate the inverse of matX'matX you don't need it in the first place. This does not involve MASS:

solve(crossprod(matX))
##              [,1]         [,2]         [,3]
## [1,]  0.513763935 -0.104219636 -0.002371406
## [2,] -0.104219636  0.038700372 -0.007798748
## [3,] -0.002371406 -0.007798748  0.006625269

svd The svd could also be used and similarly does not require MASS:

with(svd(matX), v %*% diag(1/d^2) %*% t(v))
##              [,1]         [,2]         [,3]
## [1,]  0.513763935 -0.104219636 -0.002371406
## [2,] -0.104219636  0.038700372 -0.007798748
## [3,] -0.002371406 -0.007798748  0.006625269

ADDED some additional info.

`zapsmall` might be nice for showing the results of `inv %*% matX` more cleanly ... — Ben Bolker, Jan 26 '14 at 16:26
Oh cool :D Thanks for taking the time to answer :D That's great. — Reanimation, Jan 26 '14 at 18:36
For the svd solution, why square the diagonal matrix? diag(1/d^2) — Katelynn ruan, Feb 21 '20 at 07:35
@user1143669, Please read the second last line of the question. — G. Grothendieck, Feb 21 '20 at 13:36

Marc in the box · Answer 2 · 2014-01-26T16:29:54.540

2

You can do what's called a "Moore–Penrose pseudoinverse". Here's a function exp.matthat will do this for you. There is also an example outlining it's use here.

`exp.mat()`:

#The exp.mat function performs can calculate the pseudoinverse of a matrix (EXP=-1)
#and other exponents of matrices, such as square roots (EXP=0.5) or square root of 
#its inverse (EXP=-0.5). 
#The function arguments are a matrix (MAT), an exponent (EXP), and a tolerance
#level for non-zero singular values.
exp.mat<-function(MAT, EXP, tol=NULL){
    MAT <- as.matrix(MAT)
    matdim <- dim(MAT)
    if(is.null(tol)){
        tol=min(1e-7, .Machine$double.eps*max(matdim)*max(MAT))
    }
    if(matdim[1]>=matdim[2]){ 
        svd1 <- svd(MAT)
        keep <- which(svd1$d > tol)
        res <- t(svd1$u[,keep]%*%diag(svd1$d[keep]^EXP, nrow=length(keep))%*%t(svd1$v[,keep]))
    }
    if(matdim[1]<matdim[2]){ 
        svd1 <- svd(t(MAT))
        keep <- which(svd1$d > tol)
        res <- svd1$u[,keep]%*%diag(svd1$d[keep]^EXP, nrow=length(keep))%*%t(svd1$v[,keep])
    }
    return(res)
}

example of use:

source("exp.mat.R")
X <- matrix(c(rep(1, 8),2,3,4,0,6,4,3,7,-2,-4,3,5,7,8,9,11), nrow=8, ncol=3)
iX <- exp.mat(X, -1)
zapsmall(iX %*% X) # results in identity matrix
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

edited Jan 26 '14 at 16:29

answered Jan 26 '14 at 13:54

Marc in the box

11,769
4
47
97

Hello. Thanks for your time. Is there a specific package I need to use to use the function you state above? – Reanimation Jan 26 '14 at 14:09
Oh wait, I guess you need to create that function yourself in your program? – Reanimation Jan 26 '14 at 14:14
@Reanimation - Your welcome. No, you can just load the function as is (i.e. pass the function code to your R environment). If you prefer to load a package rather than `source`this function, a similar function (`pseudoinverse`) can also be found in the package `corpcor`. The advantage to the `exp.mat` function is that it can also take other exponential values (e.g. `EXP=-0.5`). – Marc in the box Jan 26 '14 at 14:15
1

`MASS::ginv` also implements a pseudo-inverse, and comes with a standard R installation ... – Ben Bolker Jan 26 '14 at 16:26
@BenBolker - Cheers. I can't believe I didn't know about `ginv` - thanks for pointing it out. – Marc in the box Jan 26 '14 at 16:32

Calculate inverse of a non-square matrix in R

2 Answers2

exp.mat():

example of use:

`exp.mat()`: