Sparse matrix and sparse subarray subtraction returns full arrays

Specific Domains Numerics
1

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:

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

@time f(2000);

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?

2

I can’t reproduce:

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?

3

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

4

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

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

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?

5

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.

6

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!

7

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.)