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 at 11:51
  • 3
    $\begingroup$ You've entered the realm of symmetric indefinite linear systems. MinRes is your friend. $\endgroup$ – Wolfgang Bangerth Jun 9 at 15:30
  • 2
    $\begingroup$ Would we expect minres on the second system to perform better than CG on the first? $\endgroup$ – rchilton1980 Jun 9 at 19:13
  • $\begingroup$ @LedHead, let us say $m \sim 1000$ and $n \in [100K, 1M]$. $\endgroup$ – passerby51 Jun 9 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$ – Brian Borchers Jun 10 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

| cite | improve this answer | |
  • $\begingroup$ Thanks. This seems the way to go. $\endgroup$ – passerby51 Jun 19 at 18:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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