Inverse of ill-conditioned symmetric matrixLapack symmetric update $B^-1AB^-T$Diagonalization of Dense Ill Conditioned Matricescomputing the determinant of a dense nonsymmetric 100x100 matrix having very big and very small eigenvalueslapack singular matrixUnder what circumstances can two (nearly) identical sparse matrices give different solutions to Mx = b?Accuracy between ill-conditioned matrix-free vs. matrix-based operatorsPractical example of why it is not good to invert a matrixInverting big symmetric and singular matricesHow to directly compute the inverse of an ill-conditioned dense matrix

Pregnant spouse slipped abortion pills unknowingly. What would the legal ramifications be?

How can communicating in human language with an unconscious alien species be treated as an attack?

How to get rid of or reduce the gap of an equation break?

Is there a word/short phrase for "the most" of something (not necessarily the majority)?

I missed an important client meeting and hurt my standing. How can I recover?

Pass on your radiation

Retrieving a list of change events

Description of attributes beyond 5

Has a remote rogue DNS client managed to slip through the first iptables rule of a Linux DNS server?

Upgrade new bike to 11 speed or downgrade trainers wheel to 8 speed

being overqualified as a barrier for getting a job

Are the stars distributed in uniform distribution, on the celestial dome, with respect to brightness?

Taking a Switch on vacation from EU to USA. What region eshop will I be able to access?

How scammy are cashback sites?

Populate field in one layer with corresponding value from another layer

Why there isn't public transport on Christmas Day in the UK

Can a public school in the USA force a 14yr old to create a Twitter account for a passing grade?

Advent calendar

Carlsen beat a high ranking GM with 1 Nh3. Conclusions?

Specific heat for tetraatomic gas molecules

What does a reduction gearbox do in a turbofan engine?

"Can't open display" when using `scp`

Who was Wodehouse’s intended audience for Jeeves and Wooster?

Senate Impeachment rule change?

Inverse of ill-conditioned symmetric matrix

Lapack symmetric update $B^-1AB^-T$Diagonalization of Dense Ill Conditioned Matricescomputing the determinant of a dense nonsymmetric 100x100 matrix having very big and very small eigenvalueslapack singular matrixUnder what circumstances can two (nearly) identical sparse matrices give different solutions to Mx = b?Accuracy between ill-conditioned matrix-free vs. matrix-based operatorsPractical example of why it is not good to invert a matrixInverting big symmetric and singular matricesHow to directly compute the inverse of an ill-conditioned dense matrix

.everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty,.everyoneloves__bot-mid-leaderboard:empty
margin-bottom:0;

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.

edited Oct 14 at 10:52

Federico Poloni

4,69714 silver badges27 bronze badges

asked Oct 1 at 23:11

InfProbSciX

1435 bronze badges

2

$begingroup$
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$?
$endgroup$
– Alone Programmer
Oct 1 at 23:17

$begingroup$
Also these two questions might be interesting for you: math.stackexchange.com/q/16940/599776 and math.stackexchange.com/q/2735/599776
$endgroup$
– Alone Programmer
Oct 1 at 23:22

6

$begingroup$
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.
$endgroup$
– Brian Borchers
Oct 2 at 0:01

1

$begingroup$
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^-1hh^TK^-1 - K^-1)c$, which only requires a couple of linear equation solves. Still, the ill-conditioning problem would remain.
$endgroup$
– Amit Hochman
Oct 3 at 13:52

2

$begingroup$
Matlab’s inv() function uses the LDL decomposition to compute the inverse, so the inverse of a symmetric matrix comes out symmetric.
$endgroup$
– Amit Hochman
Oct 3 at 14:15

|
show 1 more comment

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.

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.

edited Oct 14 at 10:52

Federico Poloni

4,69714 silver badges27 bronze badges

asked Oct 1 at 23:11

InfProbSciX

1435 bronze badges

2

$begingroup$
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$?
$endgroup$
– Alone Programmer
Oct 1 at 23:17

$begingroup$
Also these two questions might be interesting for you: math.stackexchange.com/q/16940/599776 and math.stackexchange.com/q/2735/599776
$endgroup$
– Alone Programmer
Oct 1 at 23:22

6

$begingroup$
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.
$endgroup$
– Brian Borchers
Oct 2 at 0:01

1

$begingroup$
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^-1hh^TK^-1 - K^-1)c$, which only requires a couple of linear equation solves. Still, the ill-conditioning problem would remain.
$endgroup$
– Amit Hochman
Oct 3 at 13:52

2

$begingroup$
Matlab’s inv() function uses the LDL decomposition to compute the inverse, so the inverse of a symmetric matrix comes out symmetric.
$endgroup$
– Amit Hochman
Oct 3 at 14:15

|
show 1 more comment

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.

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.

edited Oct 14 at 10:52

Federico Poloni

4,69714 silver badges27 bronze badges

asked Oct 1 at 23:11

InfProbSciX

1435 bronze badges

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.

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.

matrix precision inverse

edited Oct 14 at 10:52

Federico Poloni

4,69714 silver badges27 bronze badges

asked Oct 1 at 23:11

InfProbSciX

1435 bronze badges

edited Oct 14 at 10:52

Federico Poloni

4,69714 silver badges27 bronze badges

asked Oct 1 at 23:11

InfProbSciX

1435 bronze badges

edited Oct 14 at 10:52

Federico Poloni

4,69714 silver badges27 bronze badges

edited Oct 14 at 10:52

Federico Poloni

4,69714 silver badges27 bronze badges

edited Oct 14 at 10:52

Federico Poloni

4,69714 silver badges27 bronze badges

asked Oct 1 at 23:11

InfProbSciX

1435 bronze badges

asked Oct 1 at 23:11

InfProbSciX

1435 bronze badges

asked Oct 1 at 23:11

InfProbSciX

1435 bronze badges

2

$begingroup$
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$?
$endgroup$
– Alone Programmer
Oct 1 at 23:17

$begingroup$
Also these two questions might be interesting for you: math.stackexchange.com/q/16940/599776 and math.stackexchange.com/q/2735/599776
$endgroup$
– Alone Programmer
Oct 1 at 23:22

6

$begingroup$
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.
$endgroup$
– Brian Borchers
Oct 2 at 0:01

1

$begingroup$
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^-1hh^TK^-1 - K^-1)c$, which only requires a couple of linear equation solves. Still, the ill-conditioning problem would remain.
$endgroup$
– Amit Hochman
Oct 3 at 13:52

2

$begingroup$
Matlab’s inv() function uses the LDL decomposition to compute the inverse, so the inverse of a symmetric matrix comes out symmetric.
$endgroup$
– Amit Hochman
Oct 3 at 14:15

|
show 1 more comment

2

$begingroup$
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$?
$endgroup$
– Alone Programmer
Oct 1 at 23:17

$begingroup$
Also these two questions might be interesting for you: math.stackexchange.com/q/16940/599776 and math.stackexchange.com/q/2735/599776
$endgroup$
– Alone Programmer
Oct 1 at 23:22

6

$begingroup$
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.
$endgroup$
– Brian Borchers
Oct 2 at 0:01

1

$begingroup$
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^-1hh^TK^-1 - K^-1)c$, which only requires a couple of linear equation solves. Still, the ill-conditioning problem would remain.
$endgroup$
– Amit Hochman
Oct 3 at 13:52

2

$begingroup$
Matlab’s inv() function uses the LDL decomposition to compute the inverse, so the inverse of a symmetric matrix comes out symmetric.
$endgroup$
– Amit Hochman
Oct 3 at 14:15

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 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 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 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^-1hh^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 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 at 14:15

|
show 1 more comment

3 Answers
3

active

oldest

votes

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 $ mathrmvec$ operator this equation would be transformed to a standard linear equation as:

$$(K^T otimes K) mathrmvec(X) = mathrmvec(B)$$

Define: $K^T otimes K = A$

Finally:

$$A mathrmvec(X) = mathrmvec(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.

edited Oct 14 at 2:14

answered Oct 1 at 23:41

Alone Programmer

9782 silver badges16 bronze badges

1

$begingroup$
Thank you for this! I'll try to implement it.
$endgroup$
– InfProbSciX
Oct 1 at 23:47

$begingroup$
@InfProbSciX There was a typo and now everything is correct.
$endgroup$
– Alone Programmer
Oct 1 at 23:55

$begingroup$
I think cond(A) = cond(K)^2 so I’m a bit surprised this approach works.
$endgroup$
– Amit Hochman
Oct 2 at 14:30

$begingroup$
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.
$endgroup$
– Federico Poloni
Oct 5 at 9:36

$begingroup$
The complexity is $O(d^6)$ (using your notation) as it is written because $A$ is $d^2times d^2$ and np.linalg.solve has cubic cost.
$endgroup$
– Federico Poloni
Oct 5 at 15:02

add a comment
|

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

$$X = K^-1hh^TK^-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 = fracX_direct-XX_direct$$

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.

edited Oct 7 at 23:23

answered Oct 5 at 20:23

Alone Programmer

9782 silver badges16 bronze badges

$begingroup$
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^-1CB^-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).
$endgroup$
– Federico Poloni
Oct 5 at 21:05

$begingroup$
@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...
$endgroup$
– Alone Programmer
Oct 5 at 21:12

$begingroup$
Sorry, the function is called np.linalg.solve(A, B). It computes $A^-1B$ 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)^-1C^T)^T$.
$endgroup$
– Federico Poloni
Oct 5 at 21:16

$begingroup$
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.
$endgroup$
– Federico Poloni
Oct 5 at 21:22

1

$begingroup$
@FedericoPoloni If we define the error as: $varepsilon = fracX_direct$. The error of your method is at least 50% lower than mine and it is really faster.
$endgroup$
– Alone Programmer
Oct 5 at 22:03

|
show 10 more comments

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.

edited Oct 8 at 11:29

answered Oct 5 at 9:49

Federico Poloni

4,69714 silver badges27 bronze badges

$begingroup$
Thank you! At the moment, I'm playing around with $alpha = 1, gamma = 10$.
$endgroup$
– InfProbSciX
Oct 5 at 9:53

$begingroup$
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.
$endgroup$
– Alone Programmer
Oct 5 at 14:33

$begingroup$
@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 $ntimes n$ matrices, and does not require vectorization, so it costs $O(n^3)$. I'll coment above about the cost of your solution.
$endgroup$
– Federico Poloni
Oct 5 at 15:05

add a comment
|

Your Answer

StackExchange.ready(function()
var channelOptions =
tags: "".split(" "),
id: "363"
;
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function()
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled)
StackExchange.using("snippets", function()
createEditor();
);

else
createEditor();

);

function createEditor()
StackExchange.prepareEditor(
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: false,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
imageUploader:
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/4.0/"u003ecc by-sa 4.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
,
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
);

);

draft saved

draft discarded

StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fscicomp.stackexchange.com%2fquestions%2f33517%2finverse-of-ill-conditioned-symmetric-matrix%23new-answer', 'question_page');

);

Post as a guest

Name

Required, but never shown

3 Answers
3

active

oldest

votes

3 Answers
3

active

oldest

votes

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 $ mathrmvec$ operator this equation would be transformed to a standard linear equation as:

$$(K^T otimes K) mathrmvec(X) = mathrmvec(B)$$

Define: $K^T otimes K = A$

Finally:

$$A mathrmvec(X) = mathrmvec(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.

edited Oct 14 at 2:14

answered Oct 1 at 23:41

Alone Programmer

9782 silver badges16 bronze badges

1

$begingroup$
Thank you for this! I'll try to implement it.
$endgroup$
– InfProbSciX
Oct 1 at 23:47

$begingroup$
@InfProbSciX There was a typo and now everything is correct.
$endgroup$
– Alone Programmer
Oct 1 at 23:55

$begingroup$
I think cond(A) = cond(K)^2 so I’m a bit surprised this approach works.
$endgroup$
– Amit Hochman
Oct 2 at 14:30

$begingroup$
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.
$endgroup$
– Federico Poloni
Oct 5 at 9:36

$begingroup$
The complexity is $O(d^6)$ (using your notation) as it is written because $A$ is $d^2times d^2$ and np.linalg.solve has cubic cost.
$endgroup$
– Federico Poloni
Oct 5 at 15:02

add a comment
|

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 $ mathrmvec$ operator this equation would be transformed to a standard linear equation as:

$$(K^T otimes K) mathrmvec(X) = mathrmvec(B)$$

Define: $K^T otimes K = A$

Finally:

$$A mathrmvec(X) = mathrmvec(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.

edited Oct 14 at 2:14

answered Oct 1 at 23:41

Alone Programmer

9782 silver badges16 bronze badges

1

$begingroup$
Thank you for this! I'll try to implement it.
$endgroup$
– InfProbSciX
Oct 1 at 23:47

$begingroup$
@InfProbSciX There was a typo and now everything is correct.
$endgroup$
– Alone Programmer
Oct 1 at 23:55

$begingroup$
I think cond(A) = cond(K)^2 so I’m a bit surprised this approach works.
$endgroup$
– Amit Hochman
Oct 2 at 14:30

$begingroup$
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.
$endgroup$
– Federico Poloni
Oct 5 at 9:36

$begingroup$
The complexity is $O(d^6)$ (using your notation) as it is written because $A$ is $d^2times d^2$ and np.linalg.solve has cubic cost.
$endgroup$
– Federico Poloni
Oct 5 at 15:02

add a comment
|

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 $ mathrmvec$ operator this equation would be transformed to a standard linear equation as:

$$(K^T otimes K) mathrmvec(X) = mathrmvec(B)$$

Define: $K^T otimes K = A$

Finally:

$$A mathrmvec(X) = mathrmvec(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.

edited Oct 14 at 2:14

answered Oct 1 at 23:41

Alone Programmer

9782 silver badges16 bronze badges

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 $ mathrmvec$ operator this equation would be transformed to a standard linear equation as:

$$(K^T otimes K) mathrmvec(X) = mathrmvec(B)$$

Define: $K^T otimes K = A$

Finally:

$$A mathrmvec(X) = mathrmvec(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.

edited Oct 14 at 2:14

answered Oct 1 at 23:41

Alone Programmer

9782 silver badges16 bronze badges

edited Oct 14 at 2:14

answered Oct 1 at 23:41

Alone Programmer

9782 silver badges16 bronze badges

answered Oct 1 at 23:41

Alone Programmer

9782 silver badges16 bronze badges

answered Oct 1 at 23:41

Alone Programmer

9782 silver badges16 bronze badges

1

$begingroup$
Thank you for this! I'll try to implement it.
$endgroup$
– InfProbSciX
Oct 1 at 23:47

$begingroup$
@InfProbSciX There was a typo and now everything is correct.
$endgroup$
– Alone Programmer
Oct 1 at 23:55

$begingroup$
I think cond(A) = cond(K)^2 so I’m a bit surprised this approach works.
$endgroup$
– Amit Hochman
Oct 2 at 14:30

$begingroup$
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.
$endgroup$
– Federico Poloni
Oct 5 at 9:36

$begingroup$
The complexity is $O(d^6)$ (using your notation) as it is written because $A$ is $d^2times d^2$ and np.linalg.solve has cubic cost.
$endgroup$
– Federico Poloni
Oct 5 at 15:02

add a comment
|

1

$begingroup$
Thank you for this! I'll try to implement it.
$endgroup$
– InfProbSciX
Oct 1 at 23:47

$begingroup$
@InfProbSciX There was a typo and now everything is correct.
$endgroup$
– Alone Programmer
Oct 1 at 23:55

$begingroup$
I think cond(A) = cond(K)^2 so I’m a bit surprised this approach works.
$endgroup$
– Amit Hochman
Oct 2 at 14:30

$begingroup$
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.
$endgroup$
– Federico Poloni
Oct 5 at 9:36

$begingroup$
The complexity is $O(d^6)$ (using your notation) as it is written because $A$ is $d^2times d^2$ and np.linalg.solve has cubic cost.
$endgroup$
– Federico Poloni
Oct 5 at 15:02

Thank you for this! I'll try to implement it.

– InfProbSciX
Oct 1 at 23:47

@InfProbSciX There was a typo and now everything is correct.

– Alone Programmer
Oct 1 at 23:55

I think cond(A) = cond(K)^2 so I’m a bit surprised this approach works.

– Amit Hochman
Oct 2 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 at 9:36

The complexity is $O(d^6)$ (using your notation) as it is written because $A$ is $d^2times d^2$ and np.linalg.solve has cubic cost.

– Federico Poloni
Oct 5 at 15:02

add a comment
|

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

$$X = K^-1hh^TK^-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 = fracX_direct-XX_direct$$

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.

edited Oct 7 at 23:23

answered Oct 5 at 20:23

Alone Programmer

9782 silver badges16 bronze badges

$begingroup$
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^-1CB^-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).
$endgroup$
– Federico Poloni
Oct 5 at 21:05

$begingroup$
@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...
$endgroup$
– Alone Programmer
Oct 5 at 21:12

$begingroup$
Sorry, the function is called np.linalg.solve(A, B). It computes $A^-1B$ 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)^-1C^T)^T$.
$endgroup$
– Federico Poloni
Oct 5 at 21:16

$begingroup$
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.
$endgroup$
– Federico Poloni
Oct 5 at 21:22

1

$begingroup$
@FedericoPoloni If we define the error as: $varepsilon = fracX_direct$. The error of your method is at least 50% lower than mine and it is really faster.
$endgroup$
– Alone Programmer
Oct 5 at 22:03

|
show 10 more comments

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

$$X = K^-1hh^TK^-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 = fracX_direct-XX_direct$$

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.

edited Oct 7 at 23:23

answered Oct 5 at 20:23

Alone Programmer

9782 silver badges16 bronze badges

$begingroup$
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^-1CB^-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).
$endgroup$
– Federico Poloni
Oct 5 at 21:05

$begingroup$
@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...
$endgroup$
– Alone Programmer
Oct 5 at 21:12

$begingroup$
Sorry, the function is called np.linalg.solve(A, B). It computes $A^-1B$ 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)^-1C^T)^T$.
$endgroup$
– Federico Poloni
Oct 5 at 21:16

$begingroup$
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.
$endgroup$
– Federico Poloni
Oct 5 at 21:22

1

$begingroup$
@FedericoPoloni If we define the error as: $varepsilon = fracX_direct$. The error of your method is at least 50% lower than mine and it is really faster.
$endgroup$
– Alone Programmer
Oct 5 at 22:03

|
show 10 more comments

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

$$X = K^-1hh^TK^-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 = fracX_direct-XX_direct$$

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.

edited Oct 7 at 23:23

answered Oct 5 at 20:23

Alone Programmer

9782 silver badges16 bronze badges

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

$$X = K^-1hh^TK^-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 = fracX_direct-XX_direct$$

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.

edited Oct 7 at 23:23

answered Oct 5 at 20:23

Alone Programmer

9782 silver badges16 bronze badges

edited Oct 7 at 23:23

answered Oct 5 at 20:23

Alone Programmer

9782 silver badges16 bronze badges

answered Oct 5 at 20:23

Alone Programmer

9782 silver badges16 bronze badges

answered Oct 5 at 20:23

Alone Programmer

9782 silver badges16 bronze badges

$begingroup$
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^-1CB^-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).
$endgroup$
– Federico Poloni
Oct 5 at 21:05

$begingroup$
@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...
$endgroup$
– Alone Programmer
Oct 5 at 21:12

$begingroup$
Sorry, the function is called np.linalg.solve(A, B). It computes $A^-1B$ 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)^-1C^T)^T$.
$endgroup$
– Federico Poloni
Oct 5 at 21:16

$begingroup$
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.
$endgroup$
– Federico Poloni
Oct 5 at 21:22

1

$begingroup$
@FedericoPoloni If we define the error as: $varepsilon = fracX_direct$. The error of your method is at least 50% lower than mine and it is really faster.
$endgroup$
– Alone Programmer
Oct 5 at 22:03

|
show 10 more comments

$begingroup$
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^-1CB^-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).
$endgroup$
– Federico Poloni
Oct 5 at 21:05

$begingroup$
@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...
$endgroup$
– Alone Programmer
Oct 5 at 21:12

$begingroup$
Sorry, the function is called np.linalg.solve(A, B). It computes $A^-1B$ 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)^-1C^T)^T$.
$endgroup$
– Federico Poloni
Oct 5 at 21:16

$begingroup$
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.
$endgroup$
– Federico Poloni
Oct 5 at 21:22

1

$begingroup$
@FedericoPoloni If we define the error as: $varepsilon = fracX_direct$. The error of your method is at least 50% lower than mine and it is really faster.
$endgroup$
– Alone Programmer
Oct 5 at 22:03

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^-1CB^-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 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 at 21:12

Sorry, the function is called np.linalg.solve(A, B). It computes $A^-1B$ 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)^-1C^T)^T$.

– Federico Poloni
Oct 5 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 at 21:22

@FedericoPoloni If we define the error as: $varepsilon = fracX_direct$. The error of your method is at least 50% lower than mine and it is really faster.

– Alone Programmer
Oct 5 at 22:03

|
show 10 more comments

I suggest you combine these two tricks.

edited Oct 8 at 11:29

answered Oct 5 at 9:49

Federico Poloni

4,69714 silver badges27 bronze badges

$begingroup$
Thank you! At the moment, I'm playing around with $alpha = 1, gamma = 10$.
$endgroup$
– InfProbSciX
Oct 5 at 9:53

$begingroup$
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.
$endgroup$
– Alone Programmer
Oct 5 at 14:33

$begingroup$
@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 $ntimes n$ matrices, and does not require vectorization, so it costs $O(n^3)$. I'll coment above about the cost of your solution.
$endgroup$
– Federico Poloni
Oct 5 at 15:05

add a comment
|

I suggest you combine these two tricks.

edited Oct 8 at 11:29

answered Oct 5 at 9:49

Federico Poloni

4,69714 silver badges27 bronze badges

$begingroup$
Thank you! At the moment, I'm playing around with $alpha = 1, gamma = 10$.
$endgroup$
– InfProbSciX
Oct 5 at 9:53

$begingroup$
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.
$endgroup$
– Alone Programmer
Oct 5 at 14:33

$begingroup$
@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 $ntimes n$ matrices, and does not require vectorization, so it costs $O(n^3)$. I'll coment above about the cost of your solution.
$endgroup$
– Federico Poloni
Oct 5 at 15:05

add a comment
|

I suggest you combine these two tricks.

edited Oct 8 at 11:29

answered Oct 5 at 9:49

Federico Poloni

4,69714 silver badges27 bronze badges

I suggest you combine these two tricks.

edited Oct 8 at 11:29

answered Oct 5 at 9:49

Federico Poloni

4,69714 silver badges27 bronze badges

edited Oct 8 at 11:29

answered Oct 5 at 9:49

Federico Poloni

4,69714 silver badges27 bronze badges

answered Oct 5 at 9:49

Federico Poloni

4,69714 silver badges27 bronze badges

answered Oct 5 at 9:49

Federico Poloni

4,69714 silver badges27 bronze badges

$begingroup$
Thank you! At the moment, I'm playing around with $alpha = 1, gamma = 10$.
$endgroup$
– InfProbSciX
Oct 5 at 9:53

$begingroup$
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.
$endgroup$
– Alone Programmer
Oct 5 at 14:33

$begingroup$
@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 $ntimes n$ matrices, and does not require vectorization, so it costs $O(n^3)$. I'll coment above about the cost of your solution.
$endgroup$
– Federico Poloni
Oct 5 at 15:05

add a comment
|

$begingroup$
Thank you! At the moment, I'm playing around with $alpha = 1, gamma = 10$.
$endgroup$
– InfProbSciX
Oct 5 at 9:53

$begingroup$
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.
$endgroup$
– Alone Programmer
Oct 5 at 14:33

$begingroup$
@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 $ntimes n$ matrices, and does not require vectorization, so it costs $O(n^3)$. I'll coment above about the cost of your solution.
$endgroup$
– Federico Poloni
Oct 5 at 15:05

Thank you! At the moment, I'm playing around with $alpha = 1, gamma = 10$.

– InfProbSciX
Oct 5 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 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 $ntimes 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 at 15:05

add a comment
|

draft saved

draft discarded

Thanks for contributing an answer to Computational Science Stack Exchange!

Please be sure to answer the question. Provide details and share your research!

But avoid …

Asking for help, clarification, or responding to other answers.

Making statements based on opinion; back them up with references or personal experience.

Use MathJax to format equations. MathJax reference.

To learn more, see our tips on writing great answers.

draft saved

draft discarded

Post as a guest

Name

Required, but never shown

Name

Required, but never shown

Name

Required, but never shown

This page is only for reference, If you need detailed information, please check here

搜尋此網誌

Bsrgvty

3 Answers
3

Your Answer

Post as a guest

3 Answers
3

3 Answers
3

Post as a guest

Popular posts from this blog

Genealogie vun de Merowenger Vum Merowech bis zum Chilperich I. | Navigatiounsmenü

3 Answers 3

Your Answer

Sign up or log in

Post as a guest

Post as a guest

3 Answers 3

3 Answers 3

Sign up or log in

Post as a guest

Post as a guest

Sign up or log in

Post as a guest

Sign up or log in

Post as a guest

Sign up or log in

Post as a guest

Popular posts from this blog

Genealogie vun de Merowenger Vum Merowech bis zum Chilperich I. | Navigatiounsmenü

3 Answers
3

3 Answers
3

3 Answers
3