# Sparse matrix and sparse subarray subtraction returns full arrays

**URL:** <https://discourse.julialang.org/t/sparse-matrix-and-sparse-subarray-subtraction-returns-full-arrays/88239>\
**Category:** Numerics\
**Tags:** sparse\
**Created:** [October 4, 2022, 3:46pm UTC](https://discourse.julialang.org/t/sparse-matrix-and-sparse-subarray-subtraction-returns-full-arrays/88239 "2022-10-04T15:46:19Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![vivian-rogers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vivian-rogers/32/26032_2.png) [@vivian-rogers](https://discourse.julialang.org/u/vivian-rogers)\
**Post date:** [October 4, 2022, 3:46pm UTC](https://discourse.julialang.org/t/sparse-matrix-and-sparse-subarray-subtraction-returns-full-arrays/88239/1 "2022-10-04T15:46:19Z")

</div>

Hi! So, I am writing a routine to recursively invert a big, sparse matrix which is known to be nice and reasonably invertible with this recursive method. Specifically, this pertains to the green’s function of a sparse hamiltonian for nanoelectronics modelling.

Unfortunately, when I take some @view of a sparse matrix, the product of a sparse matrix and the subarray return full arrays, greatly reducing the efficiency of this method.

The function is defined below, and here is a test case to profile it:

```julia
Using SparseArrays; Using LinearAlgebra; f(n) = recInv(sprand(n,n,1/n))

@time f(2000);

```

```julia
function recInv(M::Union{SparseMatrixCSC, SubArray}, sizecutoff::Int=128)
    nx = size(M)[2]
    ny = size(M)[1]
    if(nx <= sizecutoff || ny <= sizecutoff)
        return sparse(inv(Array(M)))
    else
        hxu = Int(round(nx/2))+1; hyu = Int(round(ny/2))+1
        hxl = hxu -1; hyl = hyu -1
        A = @view M[1:hyl,1:hxl]; U = @view M[1:hyl,hxu:nx];
        V = @view M[hyu:ny,1:hxl]; C = @view M[hyu:ny,hxu:nx];
        Cinv = recInv(C,sizecutoff)
        Γ = Cinv*V; Δ = sparse(U*Cinv)
        expThatBreaks = A .- Δ*V
        Σ = recInv(sparse(expThatBreaks),sizecutoff)
        Minv = [Σ -Σ*Δ;
                -Γ*Σ Cinv .+ Γ*Σ*Δ]
        return sparse(Minv)
    end
end

```

If one takes a look at expThatBreaks, Δ\*V is returning full Arrays though I see no reason why it should be. Do I need to write my own routines for sparse array - subarray multiplication?

---

<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:** [October 4, 2022, 4:44pm UTC](https://discourse.julialang.org/t/sparse-matrix-and-sparse-subarray-subtraction-returns-full-arrays/88239/2 "2022-10-04T16:44:08Z")

</div>

> [@vivian-rogers](#):
>
> when I take some @view of a sparse matrix, the product of a sparse matrix and the subarray return full arrays

I can’t reproduce:

```julia
julia> using SparseArrays

julia> A = sprand(10,10, 0.1);

julia> @views A * A[:, 1:3]
10×3 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
 0.0486689 ⋅ 0.537533
  ⋅ ⋅ ⋅ 
  ⋅ ⋅ ⋅ 
  ⋅ 0.367537 ⋅ 
  ⋅ ⋅ ⋅ 
  ⋅ ⋅ ⋅ 
  ⋅ ⋅ ⋅ 
  ⋅ ⋅ ⋅ 
  ⋅ ⋅ ⋅ 
  ⋅ ⋅ ⋅ 

```

Could you make a similar simple example that illustrates your problem?

---

<div class="post-metadata">

**Author:** ![vivian-rogers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vivian-rogers/32/26032_2.png) [@vivian-rogers](https://discourse.julialang.org/u/vivian-rogers)\
**Post date:** [October 4, 2022, 4:47pm UTC](https://discourse.julialang.org/t/sparse-matrix-and-sparse-subarray-subtraction-returns-full-arrays/88239/3 "2022-10-04T16:47:38Z")

</div>

Huh. Okay, let me mess around with this for a bit, definitely might be on my end.

---

<div class="post-metadata">

**Author:** ![vivian-rogers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vivian-rogers/32/26032_2.png) [@vivian-rogers](https://discourse.julialang.org/u/vivian-rogers)\
**Post date:** [October 10, 2022, 7:39pm UTC](https://discourse.julialang.org/t/sparse-matrix-and-sparse-subarray-subtraction-returns-full-arrays/88239/4 "2022-10-10T19:39:42Z")

</div>

Okay. I apologize, I realize my premise was false. My issue arises in the following MWE with sparse matrix / sparse subarray subtraction, not multiplication:

```julia
exmat(n) = sprand(n,n,1/n)+5I(n);
A = exmat(50); B = exmat(100);
Bview = @view B[1:50,1:50]; 
C = A .- Bview
display(typeof(A))
display(typeof(Bview))
display(typeof(C))

```

Which returns

```julia
SparseMatrixCSC{Float64, Int64}
SubArray{Float64, 2, SparseMatrixCSC{Float64, Int64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}
Matrix{Float64} (alias for Array{Float64, 2})

```

Can anything to be done to prevent the subtraction from returning a non-sparse matrix?

---

<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:** [October 10, 2022, 10:04pm UTC](https://discourse.julialang.org/t/sparse-matrix-and-sparse-subarray-subtraction-returns-full-arrays/88239/5 "2022-10-10T22:04:53Z")

</div>

> [@vivian-rogers](#):
>
> Can anything to be done to prevent the subtraction from returning a non-sparse matrix?

Not using views seems to be the current workaround: `A .- B[1:50,1:50]` works fine. Of course, this allocates a bit extra memory.

Probably there is a missing broadcast method in the SparseArrays package — can you file an issue?

PS. Note that `sprand(n,n,1/n)+5I` should be faster than `sprand(n,n,1/n)+5I(n)`. In Julia, `I` can be used as a “magic matrix” that automatically chooses its size according to what it is added to, and there are specialized methods for adding it to most matrix types. `I(n)` allocates a `Diagonal` matrix, which isn’t too bad (only the diagonal entries are stored), but it doesn’t take advantage of the fact that all of the diagonal entries are the same.

---

<div class="post-metadata">

**Author:** ![vivian-rogers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vivian-rogers/32/26032_2.png) [@vivian-rogers](https://discourse.julialang.org/u/vivian-rogers)\
**Post date:** [October 10, 2022, 11:37pm UTC](https://discourse.julialang.org/t/sparse-matrix-and-sparse-subarray-subtraction-returns-full-arrays/88239/6 "2022-10-10T23:37:53Z")

</div>

Ah well enough.That is somewhat unfortunate to hear, since I think my code is scaling with O(n^2) on account of the memory allocations.

Sure, I will file an issue then.

Thank you for the insight!

---

<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:** [October 11, 2022, 12:24am UTC](https://discourse.julialang.org/t/sparse-matrix-and-sparse-subarray-subtraction-returns-full-arrays/88239/7 "2022-10-11T00:24:49Z")

</div>

> [@vivian-rogers](#):
>
> Ah well enough.That is somewhat unfortunate to hear, since I think my code is scaling with O(n^2) on account of the memory allocations.

The copying workaround won’t hurt the asymptotic scaling, since `A + B` has to allocate anyway for the result — it only changes the constant factor if you convert the `B` view into a copy.

(If you really want to cut down on allocations for this type of code, you need to think more deeply about the CSC data structure used to store sparse matrices.)