johnhenry - 10 months ago 83

Python Question

The following code shows a problem of singularity of a matrix, since working in Pycharm I get

`raise LinAlgError("Singular matrix")`

numpy.linalg.linalg.LinAlgError: Singular matrix

I guess the problem is K but I cannot understand exactly how:

`from numpy import zeros`

from numpy.linalg import linalg

import math

def getA(kappa):

matrix = zeros((n, n), float)

for i in range(n):

for j in range(n):

matrix[i][j] = 2*math.cos((2*math.pi/n)*(abs(j-i))*kappa)

return matrix

def getF(csi, a):

csiInv = linalg.inv(csi)

valueF = csiInv * a * csiInv * a

traceF = valueF.trace()

return 0.5 * traceF

def getG(csi, f, a):

csiInv = linalg.inv(csi)

valueG = (csiInv * a * csiInv) / (2 * f)

return valueG

def getE(g, k):

KInv = linalg.inv(k)

Ktrans = linalg.transpose(k)

KtransInv = linalg.inv(Ktrans)

e = KtransInv * g * KInv

return e

file = open('transformed.txt', 'r')

n = 4

transformed = zeros(n)

for counter, line in enumerate(file):

if counter == n:

break

transformed[counter] = float(line)

CSI = zeros((n, n))

for i in range(n):

for j in range(n):

CSI[i][j] = transformed[abs(i-j)]

A = getA(1)

F = getF(CSI, A)

G = getG(CSI, F, A)

K = zeros((n, n), float)

for j in range(n):

K[0][j] = 0.0001

for i in range(1, n):

for j in range(n):

K[i][j] = ((3.0*70.0*70.0*0.3)/(2.0*300000.0*300000.0))*((j*(i-j))/i)*(1.0+(70.0/300000.0)*j)

E = getE(G, K)

print G

print K

Does anyone has any suggestions to fix it? Thank you

Answer Source

Inverting matrices that are very "close" to being singular often causes computation problems. A quick hack is to add a very small value to the diagonal of your matrix before inversion.

```
def getE(g, k):
m = 10^-6
KInv = linalg.inv(k + numpy.eye(k.shape[1])*m)
Ktrans = linalg.transpose(k)
KtransInv = linalg.inv(Ktrans + + numpy.eye(Ktrans.shape[1])*m)
e = KtransInv * g * KInv
return e
```

I think of that as being good enough for homework. But if you want to really deploy something computationally robust, you should look into alternatives to inverting.