# Sparse solve with sparse rhs

**URL:** https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595
**Category:** Numerics
**Tags:** linearalgebra, sparse
**Created:** [May 24, 2022, 2:22pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595 "2022-05-24T14:22:23Z")
**Posts on this page:** 1
**Showing post:** 3

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 24, 2022, 3:32pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/3 "2022-05-24T15:32:17Z")

</div>

(Moved to a new thread — please don’t resurrect ancient threads.)

> [@Andrea\_Vigliotti](#):
>
> This is a problem when dealing with very large, and very sparse, matrices, where it is not possible allocate memory for a full rhs

In general, even if you have a sparse rhs and a sparse matrix, the solution will be dense, so if you don’t have enough memory to allocate the rhs then you will be in trouble with the solution, no?

However, if you have a sparse rhs (or set of right-hand-sides) and you only want a _sparse subset_ of the solution vectors, then there is a nice technique recently identified by a colleague of mine ([Lin, Wang, & Hsu (2022)](https://arxiv.org/abs/2205.07887)) that may give enormous speedups. (They applied it to scattering problems, but it’s quite a general linear-algebra trick.)

Suppose you are solving AX=B for some matrix B of right-hand sides, but you are only interested in a sparse “projection” Y = PX of the solutions involving some small subset of the outputs (rows P \<\< cols P). Then, you can form the following “augmented” sparse matrix K and only do a **partial** LU factorization to introduce zeros in the block corresponding to -P:

K = \begin{pmatrix} A & B \\ -P & 0\end{pmatrix} \rightsquigarrow \begin{pmatrix} U & F \\ 0 & Y\end{pmatrix}

at which point the block Y (the [Schur complement](https://en.wikipedia.org/wiki/Schur_complement)) is exactly the desired solution Y = PA^{-1}B. They find that this can often be vastly faster than a traditional sparse solve, mainly because it doesn’t need to do a full LU factorization of A (the matrix U can be discarded, and does not even need to be upper triangular). It also fully exploits sparsity of B.

They write that MUMPS and PARDISO already support such partial factorizations, so you could presumably do this via the Julia [MUMPS.jl](https://github.com/JuliaSmoothOptimizers/MUMPS.jl) or [Pardiso.jl](https://github.com/JuliaSparse/Pardiso.jl) packages (possibly with a little digging to access the lower-level API). I haven’t checked yet whether SuiteSparse supports this kind of partial factorization.

---

_[View the full topic](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595)._