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

Question

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.

Answer 1

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