9.7 Inverse of a Real Symmetric Matrix

SymInv

## function (SymMat) 
## {
##     n = dim(SymMat)[1]
##     L = SymMat
##     Linv = diag(1, n)
##     RetMat = matrix(0, nrow = n, ncol = n)
##     for (i in 1:n) {
##         for (j in 1:i) {
##             if (j > 1) 
##                 for (k in 1:(j - 1)) L[i, j] = L[i, j] - L[i, 
##                   k] * L[j, k] * L[k, k]
##             if (j < i) 
##                 L[i, j] = L[i, j]/L[j, j]
##         }
##     }
##     for (i in 2:n) {
##         for (j in (i - 1):1) {
##             for (k in j:(i - 1)) Linv[i, j] = Linv[i, j] - L[i, 
##                 k] * Linv[k, j]
##         }
##     }
##     for (i in 1:n) {
##         for (j in 1:i) {
##             for (k in i:n) RetMat[i, j] = RetMat[i, j] + Linv[k, 
##                 i]/L[k, k] * Linv[k, j]
##             RetMat[j, i] = RetMat[i, j]
##         }
##     }
##     return(RetMat)
## }
## <bytecode: 0x000001e2b3ccfab0>
## <environment: namespace:math>

SymInv(M)

##            [,1]    [,2]       [,3]    [,4]       [,5]
## [1,]  6.061e-01 -0.5152 -2.860e-17  0.1515 -6.061e-02
## [2,] -5.152e-01  1.0379 -5.000e-01 -0.1288  1.515e-01
## [3,] -2.860e-17 -0.5000  1.000e+00 -0.5000 -2.691e-16
## [4,]  1.515e-01 -0.1288 -5.000e-01  1.0379 -5.152e-01
## [5,] -6.061e-02  0.1515 -2.691e-16 -0.5152  6.061e-01

solve(M)

##            [,1]    [,2] [,3]    [,4]       [,5]
## [1,]  6.061e-01 -0.5152  0.0  0.1515 -6.061e-02
## [2,] -5.152e-01  1.0379 -0.5 -0.1288  1.515e-01
## [3,]  1.619e-17 -0.5000  1.0 -0.5000 -6.476e-17
## [4,]  1.515e-01 -0.1288 -0.5  1.0379 -5.152e-01
## [5,] -6.061e-02  0.1515  0.0 -0.5152  6.061e-01