Solve many linear equations of similar structure

Question

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 &gt;= 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

Is there any relation between the matrices $G$ and vectors $v$ of the disctinct systems? — Ben Grossmann, May 30 '20 at 04:04
Not much. G can be a random hermitian matrix, and v represents a linear combination of two elements of G. n iterations cycles through all n combinations. — rwalsh3750, May 30 '20 at 05:56

score 4 · Accepted Answer · answered May 30 '20 at 04:15

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$.

On initial pass, this seems to checkout. One last question for reference is what is the name for this type of factoring? I want to look up this and similar materials. Other than that, I will revise my post, later today, with a basic python demonstration of the method and verify your post that this is a good solution. Thanks! — rwalsh3750, May 30 '20 at 16:01
@Walden95 I would refer to this trick/technique as block (or blockwise) Gaussian elimination or block $LU$ factorization. Another fun application of this idea is the proofs by "blockwise elimination and $LDU$ factorization of the Woodbury matrix identity. — eepperly16, May 30 '20 at 21:32
Nice! This helps a lot. Also, I've revised my initial post and further generalized your solution to account for $v$ being a NxM matrix. — rwalsh3750, May 31 '20 at 01:17

Solve many linear equations of similar structure

for eigenvalue decomposition

lambda_G, Q_G = np.linalg.eigh(A_G)

lambda_v, Q_v = np.linalg.eigh(A_v)

Solution:

Setup problem and randomize

General solution to compare against

Generalized solution created from eepperly16's solution Part 1

Generalized solution crated from eepperly16's solution Part 2

Final Edit:

1 Answers1