Solving an m x m symmetric linear system involving a matrix multiplication versus an (n+m) x (n+m) system

Suppose that $$R$$ and $$D$$ are an $$n \times m$$ and $$m \times m$$ matrices. Assume that $$m \ll n$$ and that $$D$$ is positive definite. We would like to solve the system $$(R^T R + D) x = R^T b$$. This involves the matrix multiplication $$R^T R$$ which is quite costly, $$O(m^2 n)$$ in the worst-case.

Instead, letting $$z = R x - b$$, we can try to solve $$\begin{bmatrix} D & R^T \\ R & - I_n \end{bmatrix} \begin{bmatrix} x \\ z \end{bmatrix} = \begin{bmatrix} 0 \\ b \end{bmatrix}$$ at the worst-case cost of $$(n+m)^3$$ which seems to be worse that the original approach. But the latter system has a big block of identity. Is it possible to solve the larger system more efficiently than the smaller system, say by an efficient Cholesky decomposition or an iterative procedure. Let us assume that $$R$$ is generally dense, although we might be able to approximate it by a sparse matrix, by thresholding.

I should add that, any solution of the original problem that avoids direct matrix multiplication is also interesting.

• How large are $m$ and $n$, typically? Jun 9 '20 at 11:51
• You've entered the realm of symmetric indefinite linear systems. MinRes is your friend. Jun 9 '20 at 15:30
• Would we expect minres on the second system to perform better than CG on the first? Jun 9 '20 at 19:13
• @LedHead, let us say $m \sim 1000$ and $n \in [100K, 1M]$. Jun 9 '20 at 21:08
• You might also consider Cholesky factorization for the reduced system and sparse symmetric factorization for the larger system. How badly conditioned are your systems? Jun 10 '20 at 12:56

In my opinion, you need a really good excuse to take a symmetric positive definite system and then turn it into an indefinite system.

The first thing I would try is the conjugate gradient method with your original $$R^{T}R+D=R^{T}b$$ system, except perform 3 matrix-vector products each iteration. Treat $$R^{T}R+D$$ as an operator instead of explicitly forming the matrix.

If that's too slow, you can perform the Cholesky decomposition on $$D$$.

$$R^{T}R+D=R^{T}R+LL^{T} = \begin{bmatrix}R^{T}&L\end{bmatrix}\begin{bmatrix}R&\\L^{T}&\end{bmatrix}=C^{T}C$$

From here, use conjugate gradient with two matrix-vector operations.

Finally you could compute the QR factorization of $$C$$, which essentially yields a cholesky factorization of $$R^{T}R+D$$.

$$C^{T}C=r^{T}Q^{T}Qr=r^{T}r$$

Its hard to know what the best approach is without testing your actual matrices.

Edit: corrected QR decomposition formula

• Thanks. This seems the way to go. Jun 19 '20 at 18:54