Given
- G: real and symmetric square matrix
- v: real column vector
I need to solve n linear systems of the form
\begin{align} A = \begin{pmatrix} G & v \\\ v^T & 0 \end{pmatrix}\end{align} \begin{align} Ax = b\end{align}
Where
- n is large
- G: real and symmetric square matrix, constant for all n systems
- v: real column vector, changes for each system (Combination vector where at most 2 values are nonzero)
- b: is zero column vector with exception of the last element
I want to know if there is a fast method to solve these many systems via exploiting this structure and suspect that there is a way to do this via eigenvalue decomposition of sums of hermitian matrices. However, I am unsure of how to combine the results.
I currently solve n systems via a hermitian solver which doesn't scale well.
For convenience, I provide the following equivalent python code
import numpy as np
import scipy.linalg as sp_linalg
np.set_printoptions(threshold=np.inf, linewidth=100000, precision=3, suppress=True)
N = 10 # Size of A-1
G = np.random.random(size=(N, N))
G += G.T
G *= 2
v = np.zeros((N, 1))
v[np.random.choice(N, 2)] = 1.0
A = np.block([[G, v], [v.T, 0.0]])
A_G = np.block([[G, np.zeros((N, 1))], [np.zeros((1, N+1))]])
A_v = np.block([[np.zeros((N, N)), v], [v.T, 0.0]])
b = np.concatenate((np.zeros((N, 1)), np.random.random((1,1))))
###
x = sp_linalg.solve(A, b, assume_a='sym') # General solution to compare against
###
# for eigenvalue decomposition
# lambda_G, Q_G = np.linalg.eigh(A_G)
# lambda_v, Q_v = np.linalg.eigh(A_v)
Thanks!
Solution:
I've taken the solution mentioned by eepperly16 and further generalized the problem. Now
- G: NxN random symetric matrix constant for all n systems
- v: NxM matrix of random variables
The big idea is since v is now a matrix, an inverse of $-v^\top G^{-1} v$ rather than doing a simple divide. These changes include...
- $x_2 = -y_2 / (v^\top G^{-1}v)$ Becomes $x_2 = (v^\top G^{-1}v)^{-1} -y_2$
- $x_1 = y_1 - x_2G^{-1}v$ Becomes $x_1 = y_1 - G^{-1}vx_2$
Since the result of this is always symmetric, that can be exploited with similar factorization. Note, however, that now the time complexity of the second stage expands proportionately to $O(M^2)$.
And finally the code with benchmark
import numpy as np
import scipy.linalg as sp_linalg
import timeit
np.random.seed(40)
np.set_printoptions(threshold=8, linewidth=1000, precision=3, suppress=True)
N = 100 # Size of square matrix G
M = 10 # Number of columns in v
# Setup problem and randomize
def setup_and_randomize():
# Create random symmetric matrix G on range (-1.0, 1.0)
G = 2.0 * np.random.random(size=(N, N)) - 1.0
G += G.T
G *= 0.5
# Create random rectangular matrix v on range (-1.0, 1.0)
v = 2.0 * np.random.random(size=(N, M)) - 1.0
A = np.block([[G, v], [v.T, np.zeros((M, M))]])
b_1 = np.zeros((N, 1))
b_2 = np.ones((M, 1))
b = np.concatenate((b_1, b_2), axis=0)
return A, G, v, b, b_1, b_2
# General solution to compare against
def naive_method(A, b):
return sp_linalg.solve(A, b, assume_a='sym')
# Generalized solution created from eepperly16's solution Part 1
def answer_method_precompute(G, b_1, b_2):
P, L, U = sp_linalg.lu(G, overwrite_a=True, check_finite=False)
L_inv = sp_linalg.solve_triangular(L, np.eye(N), lower=True, trans='N', overwrite_b=True)
U_inv = sp_linalg.solve_triangular(U, np.eye(N), lower=False, trans='N', overwrite_b=True)
G_inv = U_inv @ L_inv @ P.T
y_1 = G_inv @ b_1
y_2 = b_2 - v.T @ y_1
return y_1, y_2, G_inv
# Generalized solution crated from eepperly16's solution Part 2
def answer_method_main(v, y_1, y_2, G_inv):
G_inv_dot_v = G_inv @ v
# IF M >= 1 -----------------------------------------------------
B = v.T @ G_inv_dot_v
P, L, U = sp_linalg.lu(B, overwrite_a=True, check_finite=False)
L_inv = sp_linalg.solve_triangular(L, np.eye(M), lower=True, trans='N', overwrite_b=True)
U_inv = sp_linalg.solve_triangular(U, np.eye(M), lower=False, trans='N', overwrite_b=True)
B_inv = U_inv @ L_inv @ P.T
x_2 = B_inv @ -y_2
x_1 = y_1 - G_inv_dot_v @ x_2
# IF M == 1 -----------------------------------------------------
# x_2 = -y_2 / (v.T @ G_inv_dot_v)
# x_1 = y_1 - (x_2 * G_inv_dot_v)
return np.concatenate((x_1, x_2), axis=0)
if __name__ == "__main__":
# Verify Same Solution ------------------------------------------
A, G, v, b, b_1, b_2 = setup_and_randomize()
x_naive = naive_method(A, b)
y_1, y_2, G_inv = answer_method_precompute(G, b_1, b_2)
x_answer = answer_method_main(v, y_1, y_2, G_inv)
print('Naive Solution:\t', x_naive.T)
print('Final Solution:\t', x_answer.T)
# Benchmark Performance ----------------------------------------------
n_tests = 1000
A, G, v, b, b_1, b_2 = setup_and_randomize()
print('\nTimeit on naive_method', timeit.timeit('naive_method(A, b)', globals=globals(), number=n_tests))
print('Timeit on answer_precompute', timeit.timeit('answer_method_precompute(G, b_1, b_2)', globals=globals(), number=n_tests))
print('Timeit on answer_main', timeit.timeit('answer_method_main(v, y_1, y_2, G_inv)', globals=globals(), number=n_tests))
Which yields the following on my machine for 1000 iterations of N=100, M=10
Naive Solution: [[ 0.33 -1.518 0.434 ... -0.394 -0.569 0.824]]
Final Solution: [[ 0.33 -1.518 0.434 ... -0.394 -0.569 0.824]]
Timeit on naive_method 0.39002
Timeit on answer_precompute 0.46521499999999993
Timeit on answer_main 0.14545809999999992
Final Edit:
I understand that with scipy, there are better ways to compute the inverse that better tie into one of many BLAS style libraries. Below are 2 ways to compute the inverse of G that work better than the initial solution. Also, enabling more flags on the naive solver also makes that timing calculation fairer.
G_inv = sp_linalg.lu_solve(
sp_linalg.lu_factor(G, overwrite_a=True, check_finite=False),
np.eye(N), overwrite_b=True, check_finite=False)
L, D, perm = sp_linalg.ldl(G, overwrite_a=True, hermitian=True, check_finite=False)
L_inv = sp_linalg.solve_triangular(L[perm, :], np.eye(N), lower=True, trans='N', overwrite_b=True, check_finite=False)[:, perm]
G_inv = (L_inv.T / D.diagonal()) @ L_inv
Notice that $A$ can be factored as
$$ A = \begin{bmatrix} G & v \\ v^\top & 0 \end{bmatrix} = \begin{bmatrix} G &0 \\ v^\top & 1 \end{bmatrix}\begin{bmatrix} I & G^{-1}v \\ 0 & -v^\top G^{-1} v\end{bmatrix}. $$
Using this we can devise a scheme to solve $A$ for lots of different $G$'s. First, factorize $G$ using an $LU$ factorization (or a Cholesky factorization or $LDL^\top$ factorization or whatever). This requires time proportional to the cube of the size of $G$ ($O(n^3)$ operations), but once you have such a factorization you can compute $G^{-1}u$ in time proportional to the square of the size of $G$ ($O(n^2)$ operations). Now suppose you want to solve $Ax = b$. Write $x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}$, where $x_2$ is the last entry of $x$. Write
$$ y = \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} = \begin{bmatrix} I & G^{-1}v \\ 0 & -v^\top G^{-1} v\end{bmatrix}x. $$
Then we have that
$$ Ax = \begin{bmatrix} G &0 \\ v^\top & 1 \end{bmatrix}\begin{bmatrix} y_1 \\ y_2 \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \end{bmatrix}. $$
Then we have that $Gy_1 = b_1$. Use your precomputed $LU$ factorization to solve $Gy_1 = b_1$ for $y_1$. Then we have that $v^\top y_1 + y_2 = b_2$ so $y_2 = b_2 - v^\top y_1$.
Next we compute $x$ from $y$. Write
$$ \begin{bmatrix} I & G^{-1}v \\ 0 & -v^\top G^{-1} v\end{bmatrix}\begin{bmatrix}x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \end{bmatrix}. $$
Use your precomputed $LU$ factorization to compute $G^{-1}v$. Then we have that $(-v^\top G^{-1} v)x_2 = y_2$ so $x_2 = -y_2 / (v^\top G^{-1}v)$. We also have that $x_1 + x_2G^{-1}v = y_1$ so $x_1 = y_1 - x_2G^{-1}v$. We've now solved $Ax = b$ by using only two linear solves with $G$, which are much faster when we've precomputed the factorization of $G$.