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.

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

1 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

  • $\begingroup$ Thanks. This seems the way to go. $\endgroup$
    – passerby51
    Jun 19, 2020 at 18:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.