Singular matrix - 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

do you mean this?Traceback (most recent call last): File "/home/me/PP/Est/est.py", line 68, in E = getE(G, K) File "/home/me/PP/Est/est.py", line 33, in getE KInv = linalg.inv(k) File "/usr/lib/python2.7/dist-packages/numpy/linalg/linalg.py", line 520, in inv ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj) File "/usr/lib/python2.7/dist-packages/numpy/linalg/linalg.py", line 90, in _raise_linalgerror_singular raise LinAlgError("Singular matrix") numpy.linalg.linalg.LinAlgError: Singular matrix — johnhenry, Feb 01 '15 at 23:34
Yes. That tells us that `getE` is throwing the error. Specifically, that the inverse is the offending operation. — andrew, Feb 01 '15 at 23:39
Ok thank you! that was exactly what I was thinking, but do you have any idea of why? And how to fix it? — johnhenry, Feb 01 '15 at 23:44
Another possibility is to use a library that implements pseudoinverses ( http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse ). — , Feb 02 '15 at 00:34

score 7 · Accepted Answer · edited May 23 '17 at 12:26

7

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.

numerically stable inverse of a 2x2 matrix

edited May 23 '17 at 12:26

Community

1
1

answered Feb 01 '15 at 23:49

andrew

3,929
1
25
38

1

Thank you very much! But using your suggestion – johnhenry Feb 02 '15 at 00:00
def getE(g, k): m = 10 ^ (-6) KInv = linalg.inv(k + np.eye(k.size)*m) Ktrans = k.transpose() KtransInv = linalg.inv(Ktrans) e = KtransInv * g * KInv return e – johnhenry Feb 02 '15 at 00:00
I get: ValueError: operands could not be broadcast together with shapes (4,4) (16,16) – johnhenry Feb 02 '15 at 00:00
Ok I think that substituting k.size with 4 makes it work, but I still get, – johnhenry Feb 02 '15 at 00:08
raise LinAlgError("Singular matrix") numpy.linalg.linalg.LinAlgError: Singular matrix – johnhenry Feb 02 '15 at 00:08
even changing the matrix k and making it full of high values (order of magnitude 1000) does not change the situation, I still get the same error – johnhenry Feb 02 '15 at 00:39
But unfortunately it does not change the situation. Still getting the "singular matrix" error – johnhenry Feb 02 '15 at 00:44
Take a look at the new stack trace. I think you will see that the error is now being thrown by the `KtransInv = linalg.inv(Ktrans)` line. Try using the same trick we used on the first matrix inverse in `getE`. – andrew Feb 02 '15 at 00:53

Singular matrix - python

1 Answers1

Linked