Search code examples
rmatrixmatrix-inverse

Calculate inverse of a non-square matrix in R


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.


Solution

  • 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.