# Julia can't solve sparse linear system, MATLAB can

**URL:** <https://discourse.julialang.org/t/julia-cant-solve-sparse-linear-system-matlab-can/81692>\
**Category:** General Usage\
**Tags:** matlab, linearalgebra, sparse\
**Created:** [May 26, 2022, 4:26am UTC](https://discourse.julialang.org/t/julia-cant-solve-sparse-linear-system-matlab-can/81692 "2022-05-26T04:26:51Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Arrigo\_Benedetti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arrigo_benedetti/32/25545_2.png) [@Arrigo\_Benedetti](https://discourse.julialang.org/u/Arrigo_Benedetti)\
**Post date:** [May 26, 2022, 4:26am UTC](https://discourse.julialang.org/t/julia-cant-solve-sparse-linear-system-matlab-can/81692/1 "2022-05-26T04:26:52Z")

</div>

I am trying to solve a sparse linear system, but Julia 1.7.2 crashes (I believe it’s an OutOfMemoryError()):

```julia
infil> size(A)
(2938746, 6906)

infil> nnz(A)
17591409

infil> nnz(A)/prod(size(A))
0.0008667862253365853

infil> typeof(A)
SparseArrays.SparseMatrixCSC{Float64, Int64}

infil> typeof(b)
Vector{Float64} (alias for Array{Float64, 1})

infil> @which A \ b
\(A::SparseArrays.AbstractSparseMatrixCSC, B::AbstractVecOrMat) in SparseArrays at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/SparseArrays/src/linalg.jl:1538

infil> delta = A \ b

Killed

```

MATLAB can solve the same linear system in less than 4 sec:

```julia
K>> tstart = tic;
x = A \ b;
toc(tstart)
Warning: Rank deficient, rank = 6267, tol = 4.718235e+01. 
> In test_Ab (line 115)
 
Elapsed time is 3.883479 seconds.

```

I thought that Julia would automatically use the QR decomposition of A, so the fact that A is not full rank would not be a problem… Any suggestion to debug this?

---

<div class="post-metadata">

**Author:** ![RobertGregg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robertgregg/32/22105_2.png) [@RobertGregg](https://discourse.julialang.org/u/RobertGregg)\
**Post date:** [May 26, 2022, 4:50am UTC](https://discourse.julialang.org/t/julia-cant-solve-sparse-linear-system-matlab-can/81692/2 "2022-05-26T04:50:15Z")

</div>

This is the function referenced by the `@which` macro:

```julia
function \(A::AbstractSparseMatrixCSC, B::AbstractVecOrMat)
    require_one_based_indexing(A, B)
    m, n = size(A)
    if m == n
        if istril(A)
            if istriu(A)
                return \(Diagonal(Vector(diag(A))), B)
            else
                return \(LowerTriangular(A), B)
            end
        elseif istriu(A)
            return \(UpperTriangular(A), B)
        end
        if ishermitian(A)
            return \(Hermitian(A), B)
        end
        return \(lu(A), B)
    else
        return \(qr(A), B)
    end
end

```

It looks like because A is not a square matrix that it does indeed default to to QR decomposition. [Matlab](https://www.mathworks.com/help/matlab/ref/mldivide.html) says it also defaults to QR given a rectangular sparse matrix. Maybe try calling `qr()` in both languages to see what you get?

---

<div class="post-metadata">

**Author:** ![aramirezreyes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aramirezreyes/32/42573_2.png) [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)\
**Post date:** [May 26, 2022, 3:47pm UTC](https://discourse.julialang.org/t/julia-cant-solve-sparse-linear-system-matlab-can/81692/3 "2022-05-26T15:47:47Z")

</div>

Seems more like a bug. Are you sure it is running out of memory? (are you seeing increased consumption?) having a way to replicate the problem the the post would be great.

---

<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 26, 2022, 5:04pm UTC](https://discourse.julialang.org/t/julia-cant-solve-sparse-linear-system-matlab-can/81692/4 "2022-05-26T17:04:29Z")

</div>

> [@Arrigo\_Benedetti](#):
>
> I thought that Julia would automatically use the QR decomposition of A, so the fact that A is not full rank would not be a problem… Any suggestion to debug this?

Even with the QR decomposition, if the matrix is rank-deficient (i.e. if the columns are almost linearly dependent) then you have to drop some terms with some tolerance. (Google “rank-deficient least-squares”, which is what you are effectively doing.) Both Julia’s and Matlab’s QR functions drop nearly-dependent columns automatically, but the question is with what tolerance.

Julia’s `A \ b` _is_ using `qr(A) \ b` as commented above, but Matlab might be using a larger drop tolerance than [Julia’s default](https://github.com/JuliaSparse/SuiteSparse.jl/blob/8e83ebc991c85a6d1c63bb4d6d6de7d46e464488/src/spqr.jl#L142-L143). You might try `qr(A, tol=4.72e1) \ b` in Julia, passing the same tolerance that is reported by Matlab.

Though [according to Tim Davis](https://dl.acm.org/doi/10.1145/2049662.2049670), Julia is using the same formula (and the same library) as Matlab (unless they’ve changed since 2011?):

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/c/5c74e1c28c799aee17a9e4ea8a4f3748f40963cd.png)

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [May 26, 2022, 8:01pm UTC](https://discourse.julialang.org/t/julia-cant-solve-sparse-linear-system-matlab-can/81692/5 "2022-05-26T20:01:24Z")

</div>

How does `Killed` happen? Is julia actually killed by the OS or through user intervention? Maybe there is some maximum default amount of RAM programs are allowed to use, which is set higher for Matlab than Julia?

---

<div class="post-metadata">

**Author:** ![Arrigo\_Benedetti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arrigo_benedetti/32/25545_2.png) [@Arrigo\_Benedetti](https://discourse.julialang.org/u/Arrigo_Benedetti)\
**Post date:** [May 26, 2022, 10:12pm UTC](https://discourse.julialang.org/t/julia-cant-solve-sparse-linear-system-matlab-can/81692/6 "2022-05-26T22:12:46Z")

</div>

MacOS kills the Julia kernel. I have noticed that Julia running under Linux reports more useful information in similar situations, so I will try to reproduce this problem on a Linux VM. I have also tried @stevengj’s suggestion to increase `tol` but that did not work. Eventually I want to create a repo with a link to the data and file an issue with the SuiteSparse.jl package.