# 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:** 14
**Page:** 1
### Author:  [@Andrea\_Vigliotti](https://discourse.julialang.org/u/Andrea_Vigliotti)
#### Post date: [May 24, 2022, 2:22pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/1 "2022-05-24T14:22:23Z")
is there any news on this? It appears it is still not possible solving problems with sparse rhs. This is a problem when dealing with very large, and very sparse, matrices, where it is not possible allocate memory for a full rhs, and solving for each column of the rhs separately it would be too slow.
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [May 24, 2022, 3:26pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/2 "2022-05-24T15:26:19Z")
You mean sparse-matrix rhs?
---
### Author:  [@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")
(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.
---
### Author:  [@Andrea\_Vigliotti](https://discourse.julialang.org/u/Andrea_Vigliotti)
#### Post date: [May 24, 2022, 3:39pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/4 "2022-05-24T15:39:08Z")
yes, the situation is the type
\mathbf{A} \mathbf{X} = \mathbf{Y}
with \mathbf{A} and \mathbf{Y} sparse matrices, where \mathbf{Y} has fewer columns than \mathbf{A} but still too many for storing the full matrix. But, as @stevengj, just pointed out in general \mathbf{X} will not be sparse even if \mathbf{A} and \mathbf{Y} are.
I have a large FE problem and I would like to solve it in chunks, but plain sectioning of the domain will not solve memory problems, on the contrary it seems it will make it worse.
---
### Author:  [@Andrea\_Vigliotti](https://discourse.julialang.org/u/Andrea_Vigliotti)
#### Post date: [May 24, 2022, 3:41pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/5 "2022-05-24T15:41:03Z")
thanks I will look into it, I realized just after posting the message that I was not going to be able to store the result anyway even if I could solve the linear system.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 24, 2022, 5:50pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/6 "2022-05-24T17:50:04Z")
> [@stevengj](#):
>
> (They applied it to scattering problems, but it’s quite a general linear-algebra trick.)
Apparently @miles.lubin used a similar trick in a [2014 paper](https://epubs.siam.org/doi/10.1137/130908737) (equations 3.11–3.13), where they also commented that it’s surprisingly little known in the context of direct solves with sparse inputs/outputs.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 24, 2022, 5:54pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/7 "2022-05-24T17:54:28Z")
> [@stevengj](#):
>
> He writes 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.
Wade Hsu tells me that SuiteSparse probably cannot do this efficiently. His group currently uses MUMPS, which (because it is multifrontal) is optimized for this sort of thing—in fact, both MUMPS and PARDISO can return the desired Schur complement Y without even _storing_ the other factors U and F. See the [`mumps_schur_complement`](https://github.com/JuliaSmoothOptimizers/MUMPS.jl/blob/23a422732fa870dfff812643745aa61ba5ec65ec/src/convenience.jl#L181-L199) function in MUMPS.jl.
---
### Author:  [@fph](https://discourse.julialang.org/u/fph)
#### Post date: [November 14, 2022, 11:22pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/8 "2022-11-14T23:22:14Z")
Isn’t this equivalent to non-backward-stable Gauss-Jordan?
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [November 14, 2022, 11:27pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/9 "2022-11-14T23:27:10Z")
No.
---
### Author:  [@fph](https://discourse.julialang.org/u/fph)
#### Post date: [November 15, 2022, 8:32am UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/10 "2022-11-15T08:32:12Z")
Well that’s a short answer. 🙂 I disagree, and I have sample code to show it:
```julia
using LinearAlgebra
using Random
Random.seed!(0)
n = 5
A = rand(n,n)
b = rand(n)
# Gauss-Jordan
U = copy(A)
w = copy(b)
for k = 1:n
ind = vcat(1:k-1, k+1:n)
mult = U[ind,k] / U[k,k]
U[ind,:] = U[ind,:] - mult * U[k,:]'
w[ind] = w[ind] - mult*w[k]
# diagonal scaling of the pivot row (this is usually done at the end, but it is equivalent to do it now)
U[k,:] = U[k,:] / U[k,k]
w[k,:] = w[k,:] / U[k,k]
display(U)
display(mult)
end
display(norm(A*x - b) / norm(b))
# OP's method with P=I (one can omit rows in the second half of the matrix to change P)
K = [A b; -I zeros(n)]
U = copy(K)
for k = 1:n
mult = U[k+1:end,k] / U[k,k]
U[k+1:end,:] = U[k+1:end,:] - mult * U[k, :]'
display(U[k+1:n+k, 1:end-1])
display(mult[1:n-1])
end
x2 = U[n+1:end, end]
display(norm(A*x2 - b) / norm(b))
```
One can see that the two methods compute the exact same multipliers and entries in the elimination matrix, just with the rows in a different order.
The stability of Gauss–Jordan elimination is studied in detail, for instance, in Higham’s _Accuracy and stability of numerical algorithms_, Section 14.4. The algorithm is forward stable, but not backward stable.
The reason for this instability, informally, is that, even after pivoting, you have no control over the magnitude of the multipliers that create the lower half of K during the elimination.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [November 15, 2022, 1:46pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/11 "2022-11-15T13:46:17Z")
You _can_ use (unstable) unpivoted elimination to compute the Schur complement, but you don’t _have_ to. You are free to use as many row permutations as you want _within the A block_.
(_A_ is the matrix you would be solving with in the original form of the problem, so it’s exactly equivalent to solves using _A_ alone. And of course, sparse-direct solvers do use pivoting.)
The point here, by the way, is not the Schur complements are a new idea — this is indeed first-semester linear-algebra stuff! — but rather that people don’t seem to have appreciated that this is a powerful way to _exploit sparse outputs and right-hand sides_ in multifrontal sparse direct solvers.
---
### Author:  [@fph](https://discourse.julialang.org/u/fph)
#### Post date: [November 15, 2022, 4:22pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/12 "2022-11-15T16:22:37Z")
Sorry, I made a terrible job at explaining myself. That code was not meant to display a numerical problem; it was meant to display that the quantities you compute in that algorithm are the same as in GJ. I left off pivoting since it wasn’t the main point.
My point is that this method is _not_ equivalent to a system solve using LU factorization; you make the same LU steps, but after those you do something different than backward substitution to complete the algorithm, and this ‘something different’ is forward but not backward stable. Of course you may decide that this is sufficient stability for your large-scale applications (especially since you don’t compute all entries of the solution, so you can’t even form the residual Ax-b to check if it is small), but it is an interesting remark.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [November 15, 2022, 7:42pm UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/13 "2022-11-15T19:42:57Z")
You don’t (need to) do the “ upward” elimination steps you would do with GJ. Instead, you can just do (partial) LU with a few more rows.
---
### Author:  [@amontoison](https://discourse.julialang.org/u/amontoison)
#### Post date: [November 29, 2022, 7:32am UTC](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/14 "2022-11-29T07:32:34Z")
[BasicLU.jl](https://github.com/JuliaSmoothOptimizers/BasicLU.jl) supports sparse right-hand side for information.
```julia
m = 1000
F = basiclu_object(m)
A = sprand(m, m, 5e-3) + I
err = factorize(F, A)
b = sparsevec([1], [1.0], m)
x = solve(F, b, 'N') # sparse solution
```
Référence: [L. Schork and J. Gondzio, Permuting Spiked Matrices to Triangular Form and its Application to the Forrest-Tomlin Update](https://www.maths.ed.ac.uk/~gondzio/reports/LU-update.html)