# Inverse of ill-conditioned symmetric matrix

I've got a matrix K, with dimensions $$(n, n)$$ where each element is computed using the following equation:

$$K_{i, j} = \exp(-\alpha t_i^2 -\gamma(t_i - t_j)^2 - \alpha t_j^2)$$

... where $$t$$ is a sequence of $$(n)$$ numbers equally spaced within the interval $$[-3, 3]$$. This matrix is symmetric, so I'd expect its inverse to be symmetric as well.

Inverting this matrix is difficult due to how quickly the elements tend to zero, but if one adds a small positive number to the diagonal, base R and numpy manage to invert the matrix.

The problem is that this inverse that's computed is not symmetric.

I assume that this might be due to precision issues. Moreover, subtracting the transpose of $$K^{-1}$$ from $$K^{-1}$$ yields some pretty large values (which makes sense - if you have very small values in $$K$$, you'd expect large values in $$K^{-1}$$), but this causes problems. Is there a way to calculate the correct inverse matrix (i.e. one that's symmetric and actually yields something very close to $$K$$ when inverted again) perhaps by using some special library? I also don't mind making minor numerical changes to $$K$$ as long as it remains symmetric.

Edit: Sympy supports inversion of matrices with arbitrary precision, but the vec trick in the answer below and the comment about matlab's inv function are very interesting.

• Why you need the inverse? Taking inverse is usually too expensive and sometimes it won't work for large matrices... Do you want to solve some sort of linear equation by knowing the inverse of $K$? – Alone Programmer Oct 1 '19 at 23:17
• Also these two questions might be interesting for you: math.stackexchange.com/q/16940/599776 and math.stackexchange.com/q/2735/599776 – Alone Programmer Oct 1 '19 at 23:22
• Your inverse is difficult to compute because the matrix is nearly singular- this means that even if you use a method that forces symmetry in the inverse, the inverse will be extremely unstable. It's important that you understand the consequences of this ill-conditioning. – Brian Borchers Oct 2 '19 at 0:01
• If the weighting happens to be separable into a row weight, $r$, and a column weight, $c$, you can compute the weighted sum as $r^T(K^{-1}hh^TK^{-1} - K^{-1})c$, which only requires a couple of linear equation solves. Still, the ill-conditioning problem would remain. – Amit Hochman Oct 3 '19 at 13:52
• Matlab’s inv() function uses the LDL decomposition to compute the inverse, so the inverse of a symmetric matrix comes out symmetric. – Amit Hochman Oct 3 '19 at 14:15

You don't need the inverse even with the goal of finding $$K^{-1} h h^{T} K^{-1} - K^{-1}$$. If you are interested to have this expression, I would explain how you can convert it to a matrix equation and then solve it more efficiently:

Let's define the $$X$$ as:

$$X = K^{-1} h h^{T} K^{-1} - K^{-1}$$

Your objective is to calculate $$X$$ in this equation by assuming that $$K$$, $$h$$, and $$h^{T}$$ are all known:

$$KXK = hh^{T}-K$$

I define:

$$B = hh^{T} - K$$

By using $$\mathrm{vec}$$ operator this equation would be transformed to a standard linear equation as:

$$(K^{T} \otimes K) \mathrm{vec}(X) = \mathrm{vec}(B)$$

Define: $$K^{T} \otimes K = A$$

Finally:

$$A \mathrm{vec}(X) = \mathrm{vec}(B)$$

There are numerous efficient linear solver and as you mentioned in the comment that the dimension of $$K$$ is not that high, it should be fairly efficient.

Update:

This Python code compare the proposed method with direct way to calculate $$X$$:

import numpy as np

d = 10

K = np.zeros((d,d))

alpha = 1
gamma = 10

t_vec = np.linspace(-3,3,d)

for i in range(d):
for j in range(d):
K[i][j] = np.exp(-alpha*(t_vec[i]**2)-alpha*(t_vec[j]**2)-gamma*((t_vec[i]-t_vec[j])**2))

A = np.kron(K.T,K)

h = np.random.rand(d,d)

B = np.matmul(h,h.T) - K

vecB = B.flatten(order='F')

vecX = np.linalg.solve(A,vecB)

X = vecX.reshape((d,d))

print X

Kinv = np.linalg.inv(K)

X_direct = np.matmul(np.matmul(Kinv,np.matmul(h,h.T)),Kinv) - Kinv

print X_direct

print np.abs(X-X_direct)


The outputs are:

X:


[[1.59909954e+16 1.82624715e+12 4.24902134e+10 4.53728217e+08
2.05584992e+08 2.44120203e+08 5.98264287e+08 3.43515584e+10
3.12705250e+12 1.31655477e+16]
[1.82624715e+12 5.09588697e+09 4.78071076e+06 2.16534890e+05
9.59540287e+04 2.92820805e+04 2.33641869e+05 1.22037867e+07
2.91201539e+09 3.74298362e+12]
[4.24902134e+10 4.78071076e+06 1.65187655e+05 2.15188820e+03
8.22687923e+02 8.22528897e+02 2.14605722e+03 1.18784215e+05
1.27880035e+07 3.05868993e+10]
[4.53728217e+08 2.16534890e+05 2.15188820e+03 3.48340497e+01
1.33228000e+01 1.04249016e+01 3.45272712e+01 2.08247428e+03
2.62962350e+05 3.37679580e+08]
[2.05584992e+08 9.59540287e+04 8.22687923e+02 1.33228000e+01
4.33655487e+00 4.31336472e+00 1.42180232e+01 7.12777365e+02
9.01300211e+04 2.01206353e+08]
[2.44120203e+08 2.92820805e+04 8.22528897e+02 1.04249016e+01
4.31336472e+00 4.81918386e+00 1.18107131e+01 7.69453357e+02
7.36892858e+04 2.10911516e+08]
[5.98264287e+08 2.33641869e+05 2.14605722e+03 3.45272712e+01
1.42180232e+01 1.18107131e+01 4.00277702e+01 1.87646704e+03
1.95000703e+05 5.95359066e+08]
[3.43515584e+10 1.22037867e+07 1.18784215e+05 2.08247428e+03
7.12777365e+02 7.69453357e+02 1.87646704e+03 1.41229229e+05
1.68594796e+07 2.65035020e+10]
[3.12705250e+12 2.91201539e+09 1.27880035e+07 2.62962350e+05
9.01300211e+04 7.36892858e+04 1.95000703e+05 1.68594796e+07
3.77226961e+09 2.08638514e+12]
[1.31655477e+16 3.74298362e+12 3.05868993e+10 3.37679580e+08
2.01206353e+08 2.10911516e+08 5.95359066e+08 2.65035020e+10
2.08638514e+12 1.53843211e+16]]


X_direct:

[[1.59909954e+16 1.82624715e+12 4.24902134e+10 4.53728217e+08
2.05584992e+08 2.44120203e+08 5.98264287e+08 3.43515584e+10
3.12705250e+12 1.31655477e+16]
[1.82624715e+12 5.09588697e+09 4.78071076e+06 2.16534890e+05
9.59540287e+04 2.92820805e+04 2.33641869e+05 1.22037867e+07
2.91201539e+09 3.74298362e+12]
[4.24902134e+10 4.78071076e+06 1.65187655e+05 2.15188820e+03
8.22687923e+02 8.22528897e+02 2.14605722e+03 1.18784215e+05
1.27880035e+07 3.05868993e+10]
[4.53728217e+08 2.16534890e+05 2.15188820e+03 3.48340497e+01
1.33228000e+01 1.04249016e+01 3.45272712e+01 2.08247428e+03
2.62962350e+05 3.37679580e+08]
[2.05584992e+08 9.59540287e+04 8.22687923e+02 1.33228000e+01
4.33655487e+00 4.31336472e+00 1.42180232e+01 7.12777365e+02
9.01300211e+04 2.01206353e+08]
[2.44120203e+08 2.92820805e+04 8.22528897e+02 1.04249016e+01
4.31336472e+00 4.81918386e+00 1.18107131e+01 7.69453357e+02
7.36892858e+04 2.10911516e+08]
[5.98264287e+08 2.33641869e+05 2.14605722e+03 3.45272712e+01
1.42180232e+01 1.18107131e+01 4.00277702e+01 1.87646704e+03
1.95000703e+05 5.95359066e+08]
[3.43515584e+10 1.22037867e+07 1.18784215e+05 2.08247428e+03
7.12777365e+02 7.69453357e+02 1.87646704e+03 1.41229229e+05
1.68594796e+07 2.65035020e+10]
[3.12705250e+12 2.91201539e+09 1.27880035e+07 2.62962350e+05
9.01300211e+04 7.36892858e+04 1.95000703e+05 1.68594796e+07
3.77226961e+09 2.08638514e+12]
[1.31655477e+16 3.74298362e+12 3.05868993e+10 3.37679580e+08
2.01206353e+08 2.10911516e+08 5.95359066e+08 2.65035020e+10
2.08638514e+12 1.53843211e+16]]


np.abs(X - X_direct):

[[4.00000000e+00 4.88281250e-04 1.52587891e-05 0.00000000e+00
2.98023224e-08 0.00000000e+00 0.00000000e+00 1.14440918e-05
4.88281250e-04 6.00000000e+00]
[1.46484375e-03 3.81469727e-06 1.02445483e-08 8.73114914e-11
2.91038305e-11 2.54658516e-11 2.91038305e-11 5.58793545e-09
9.53674316e-07 9.76562500e-04]
[0.00000000e+00 3.72529030e-09 5.82076609e-11 0.00000000e+00
3.41060513e-13 3.41060513e-13 9.09494702e-13 5.82076609e-11
1.86264515e-09 3.81469727e-06]
[5.96046448e-08 1.45519152e-10 4.54747351e-13 2.13162821e-14
7.10542736e-15 1.77635684e-15 1.42108547e-14 0.00000000e+00
1.16415322e-10 5.96046448e-08]
[2.98023224e-08 2.91038305e-11 1.13686838e-13 7.10542736e-15
1.77635684e-15 3.55271368e-15 3.55271368e-15 2.27373675e-13
4.36557457e-11 2.98023224e-08]
[1.19209290e-07 2.91038305e-11 2.27373675e-13 1.77635684e-15
4.44089210e-15 1.77635684e-15 1.77635684e-15 1.13686838e-13
4.36557457e-11 2.98023224e-08]
[0.00000000e+00 2.91038305e-11 4.54747351e-13 0.00000000e+00
3.55271368e-15 1.77635684e-15 7.10542736e-15 0.00000000e+00
5.82076609e-11 2.38418579e-07]
[7.62939453e-06 5.58793545e-09 1.45519152e-11 9.09494702e-13
1.13686838e-13 2.27373675e-13 4.54747351e-13 5.82076609e-11
3.72529030e-09 0.00000000e+00]
[0.00000000e+00 9.53674316e-07 3.72529030e-09 1.16415322e-10
1.45519152e-11 4.36557457e-11 8.73114914e-11 3.72529030e-09
1.90734863e-06 2.44140625e-04]
[2.00000000e+00 4.88281250e-04 7.62939453e-06 1.19209290e-07
5.96046448e-08 2.98023224e-08 2.38418579e-07 3.81469727e-06
1.70898438e-03 2.00000000e+00]]


Which you can see the difference is pretty small and shows that the proposed method actually works.

• Thank you for this! I'll try to implement it. – InfProbSciX Oct 1 '19 at 23:47
• @InfProbSciX There was a typo and now everything is correct. – Alone Programmer Oct 1 '19 at 23:55
• I think cond(A) = cond(K)^2 so I’m a bit surprised this approach works. – Amit Hochman Oct 2 '19 at 14:30
• I have strong doubts about the efficiency of this method. Isn't this $O(n^6)$ instead of $O(n^3)$? As soon as you start ramping up the dimensions it's going to be very slow. – Federico Poloni Oct 5 '19 at 9:36
• The complexity is $O(d^6)$ (using your notation) as it is written because $A$ is $d^2\times d^2$ and np.linalg.solve has cubic cost. – Federico Poloni Oct 5 '19 at 15:02

Based on Federico's suggestions and ideas, more straight forward formulation of extracting $$K^{-1}hh^{T}K^{-1}-K^{-1}$$ would be:

$$X = K^{-1}hh^{T}K^{-1}-K^{-1}$$

$$KXK = hh^{T}-K$$

$$Z = XK$$

Solve for $$Z$$:

$$KZ = hh^{T} - K$$

and then find $$X^{T}$$ from:

$$K^{T} X^{T} = Z^{T}$$

and finally $$X$$:

$$X = (X^{T})^{T}$$

Let's define the error between direct inversion of matrices and my initial proposed method and Federico's method as:

$$\varepsilon = \frac{||X_{direct}-X||_{F}}{||X_{direct}||_{F}}$$

Where $$||\cdot||_{F}$$ is Frobenius norm.

I wrote this code based on numpy:

import numpy as np
import time

d = 100

K = np.zeros((d,d))

alpha = 1
gamma = 10

t_vec = np.linspace(-3,3,d)

for i in range(d):
for j in range(d):
K[i][j] = np.exp(-alpha*(t_vec[i]**2)-alpha*(t_vec[j]**2)-gamma*((t_vec[i]-t_vec[j])**2))

A = np.kron(K.T,K)

h = np.random.rand(d,d)

B = np.matmul(h,h.T) - K

vecB = B.flatten(order='F')

start = time.time()

vecX = np.linalg.solve(A,vecB)

end = time.time()

print "My method time elapsed: " + str(end - start) + " seconds"

X = vecX.reshape((d,d))

Kinv = np.linalg.inv(K)

X_direct = np.matmul(np.matmul(Kinv,np.matmul(h,h.T)),Kinv) - Kinv

epsilon1 = np.linalg.norm(X-X_direct) / np.linalg.norm(X_direct)

start = time.time()

Z = np.linalg.solve(K,B)

X = np.linalg.solve(K.T,Z.T).T

end = time.time()

epsilon2 = np.linalg.norm(X-X_direct) / np.linalg.norm(X_direct)

print "Federico's method time elapsed: " + str(end - start) + " seconds"

print "My method error: " + str(epsilon1)

print "Federico's method error: " + str(epsilon2)


and these are the results:

My method time elapsed: 19.122369051 seconds
Federico's method time elapsed: 0.000936031341553 seconds
My method error: 0.99999999989649
Federico's method error: 0.6635256191409429


You see that Federico's method is about 5 orders of magnitude faster and also its error is about the half of my proposed method.

• Sorry, but you seem completely off the track to me. Bartels-Stewart solves equations of the form $AX+XB=C$, and can be generalized to $AXB+CXD=E$. This problem is much simpler, $AXB=C$, with only one summand. You don't need anything fancy to solve it, just use the formula $X=A^{-1}CB^{-1}$: compute the two inverses, and carry out the products inv(A) @ C @ inv(B); this takes $O(d^3)$. Or, slightly better, use two linsolves instead of the inverses: Z = (np.linalg.linsolve(B.T, C.T)).T; X =np.linalg.linsolve(A, Z) (or something like that; I don't use numpy a lot so the syntax may be off). – Federico Poloni Oct 5 '19 at 21:05
• @FedericoPoloni I think you don't get the primary concept here: Avoid the matrix inverse computation. Your first method: two inverses are completely useless... Second proposition: There is no np.linalg.linsolve class and I don't get your formula, please elaborate it. And finally: Yes, we could do it with Bartels-Stewart algorithm and it's not that fancy and it's completely relevant here when: $C = 0$ and $D = 0$ and would be really more efficient than my initial proposition. I'm not sure what do you want to accomplish here cause you don't propose your method. Just guessing... – Alone Programmer Oct 5 '19 at 21:12
• Sorry, the function is called np.linalg.solve(A, B). It computes $A^{-1}B$ using a LU decomposition rather than an explicit inverse. Since I need $Z=CB^{-1}$, with an inverse on the right side, and I'm not sure if that is available in Numpy, I rewrote it as $Z=((B^T)^{-1}C^T)^T$. – Federico Poloni Oct 5 '19 at 21:16
• Anyway, even if Bartels-Stewart would work here, the problem is that it requires Schur decompositions, which do cost $O(d^3)$ but the multiplicative constant hidden in the big-O notation is much larger than the cost of solving two linear systems. – Federico Poloni Oct 5 '19 at 21:22
• @FedericoPoloni If we define the error as: $\varepsilon = \frac{|| X_{direct} - X||_{F}}{||X_{direct}||_{F}}$. The error of your method is at least 50% lower than mine and it is really faster. – Alone Programmer Oct 5 '19 at 22:03

My first attempt would be pulling out diagonal scaling: $$K_{i, j} = \exp(-\alpha t_i^2 -\gamma(t_i - t_j)^2 - \alpha t_j^2) = \exp(-\alpha t_i^2)\exp(-\gamma(t_i - t_j)^2)\exp(-\alpha t_j^2),$$ so $$K = DMD$$, where $$D$$ is diagonal with $$D_{ii} = \exp(-\alpha t_i^2))$$ and $$M_{ij} = \exp(-\gamma(t_i - t_j)^2)$$. Then you can rearrange the sum and reduce to inverting $$M$$ instead of $$K$$. You could try expanding the square and pulling out the terms $$\exp(-\gamma t_i^2)$$ as well, but I'm not sure if the resulting matrix will be easier to invert, because doing as I suggested $$M$$ has ones on the diagonal and rapidly decaying elements outside, which looks like a nice structure.

I haven't tried actually doing experiments, and more importantly you didn't specify the values of $$\alpha$$ and $$\gamma$$ which is important information, but I would guess that this gives a huge improvement.

[EDIT: I have quickly checked the conditioning of the resulting matrix for $$\alpha=1,\gamma=10$$; it decreases from ca. $$10^{21}$$ to ca. $$10^{18}$$, numerically; so it seems like things improve but not by much. (And I wouldn't trust those numbers anyway, since they are close to the inverse of machine precision.) What is more important, though, is that if I am not mistaken $$M$$ is a known matrix, a so-called Gaussian Toeplitz matrix. There is a closed-form expression for its Cholesky factorization, which might help you inverting it.]

Another rearranging that may help is $$X = K^{-1} h h^{T} K^{-1} - K^{-1} = K^{-1}(hh^T-K)K^{-1}$$, which allows you to use linsolve and LDL decompositions rather than explicit inverses (which, as you should know, are almost always a bad idea).

I suggest you combine these two tricks.

• Thank you! At the moment, I'm playing around with $\alpha = 1, \gamma = 10$. – InfProbSciX Oct 5 '19 at 9:53
• I'm not sure what did you add here besides decomposing diagonal and symmetric parts? Also based on your comment above I'm not sure how you reduced the complexity from $O(n^{6})$ to $O(n^{3})$ here just by decomposing diagonal and symmetric parts? I think my answer's complexity is not $O(n^{6})$ or $O(n^{3})$, it is $O(n^{4})$ and based on OP's comment the $K$ is not that big to be concerned about making the problem too expensive by applying vectorization. – Alone Programmer Oct 5 '19 at 14:33
• @AloneProgrammer Yes, my main contribution here is "decomposing diagonal and symmetric part", which I believe could have a big impact numerically, at least for large values of $\alpha$. The expression I suggested can be implemented with np.linalg.solve by working only on $n\times n$ matrices, and does not require vectorization, so it costs $O(n^3)$. I'll coment above about the cost of your solution. – Federico Poloni Oct 5 '19 at 15:05