6
$\begingroup$

I am familiar with Cholesky decomposition and LU factorization for solving systems of linear equations.

I have a problem where I have large sparse matrices (say, 1000x1000 or larger) where only one or two rows/columns (last couple rows/columns) break symmetry. I understand that Cholesky is nearly twice as efficient as LU, but I can't use Cholesky due to these outlier rows/columns.

Is there another method or algorithm that can handle this particular case - i.e., take advantage of the fact that nearly the entire matrix is self-adjoint/symmetric (real values only here)?

$\endgroup$
2
  • 1
    $\begingroup$ Performing an LU decomposition $O(2/3 n^{3})$ on a system with $O(n)$ entries is over expensive. As a consequence, Stationary method such as Jacobi iteration, Gauss-Seidel, and SOR method were developed but these methods were popular in 1950s and are considered old-fashion. Now you could resort to the highly advanced Krylov methods to solve sparse systems. On MATLAB you can use the command bicgstab() to solve this sparse system mathworks.com/help/matlab/ref/bicgstab.html $\endgroup$ Aug 22, 2021 at 10:35
  • 1
    $\begingroup$ Various strategies come to mind, including the one suggested in @StevenRoberts 's answer and Schur complements, but I think that the real challenge is coming up with something that is provably stable (or, at least, as unstable as solving the original system with LU factorization): if at any point you invert another matrix instead of $A$, that matrix may be ill-conditioned. $\endgroup$ Aug 22, 2021 at 19:41

1 Answer 1

6
$\begingroup$

Let's express the matrix $A \in \mathbb{R}^{n \times n}$ with which we want to solve linear systems as $$ A = S + U V $$ where $S$ is a symmetric matrix, $U \in \mathbb{R}^{n \times r}$, and $V \in \mathbb{R}^{r \times n}$. That is, $U V$ is a low rank update to account for the lack of symmetry. From your question, it appears $r$ is just 1 or 2.

The Woodbury matrix identity tells us that $$ A^{-1} x = (I_{n \times n} - \underbrace{S^{-1} U (I_{r \times r} + V S^{-1} U)^{-1}}_{W} V) S^{-1} x = (I_{n \times n} - W V) S^{-1} x. $$ Now, inverses appear for the symmetric matrix $S$. During the preprocessing step, you can compute the Cholesky factorization of $S$ and precompute $W = S^{-1} U (I_{r \times r} + V S^{-1} U)^{-1}$. Note that $W$ only requires $\mathcal{O}(r n^2)$ flops to compute, and that should be neglidgible.

After preprocessing, a linear system can be solved with $(I_{n \times n} - W V) S^{-1} x$. Working right to left, we can use forward/back substitution to compute $S^{-1} x$. Then, we need to perform a matrix-vector multiplication with $\left( I_{n \times n} - W V \right)$. At a cost of $\mathcal{O}(r n)$, that is neglidgible.

$\endgroup$

Your Answer

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

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