Sparse matrix plus its adjoint is converted to a dense array

[PetrKryslUCSD]

I ran into the following reason of a major slowdown with

Julia Version 0.7.0-DEV.3207
Commit 4d166fa76d* (2017-12-28 21:38 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-6650U CPU @ 2.20GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, skylake)
Environment:

Adding together a (very) sparse matrix with the adjoint of this sparse matrix, S+S', resulted in a dense array, with attendant memory and processing hit. Adding together a sparse matrix with its TRANSPOSE does not get converted to a dense array and remains sparse.

Is this the intended behavior? I don't see why it should be so…

[Sacha0]

Can you provide a MRE? Thanks! :) (Edit: Of both the S + S' and contrasting transpose behavior that you mention.)

[fredrikekre]

julia> A = sprand(2,2,1.);

julia> A + A'
2×2 Array{Float64,2}:
 1.16695   0.967356 
 0.967356  0.0754022

[fredrikekre]

It seems to be a more general problem with broadcasting with the new Adjoint and Transpose wrappers:

julia> A = sprand(2,2,0.0);

julia> A*0
2×2 SparseMatrixCSC{Float64,Int64} with 0 stored entries

julia> Transpose(A)*0
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0

julia> Adjoint(A)*0
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0

[Sacha0]

Yes, at the moment I expect S + Adjoint(S) and S + Transpose(S) to hit AbstractArray fallbacks for + that yield Array, whereas S + adjoint(S) and S + transpose(S) should perform adjoint(S)/transpose(S) eagerly and then hit a specialized + method yielding SparseMatrixCSC. (I expect this kind of behavior for most operations involving Adjoint/Transpose at the moment, including e.g. the matrix * scalar behavior illustrated above.)

I conjecture that the OP meant S + Adjoint(S) (via S + S') and S + transpose(S) behave as just described (rather than S + Adjoint(S) and S + Transpose(S) behaving differently); the purpose of my request for an MRE is to check that conjecture, because if S + Transpose(S) does not behave like S + Adjoint(S) then something curious is going on :).

[PetrKryslUCSD]

I would expect both operations of transpose and conjugate transpose (adjoint) performed on a SPARSE matrices to result in a SPARSE matrices. That was the point of my post.

[Sacha0]

Yes, ultimately that behavior should obtain. Independent of future expected behavior, my question is which of the two present transpose operations you refer to in the OP: S + transpose(S), or S + Transpose(S)? The former I imagine?

[KristofferC]

Presumably whatever S + S' lowers to. I think there are very few people who have any idea what the difference between Transpose and transpose is yet (inluding me).

[PetrKryslUCSD]

In my original code the S+S' statement produced a dense array. I corrected this by S+transpose(S). I haven't tried Transpose, but I believe Fredrik illustrated that above.

[Sacha0]

Thank you for confirming :).

[Sacha0]

Presumably whatever S + S' lowers to.

(S + S' effectively lowers to S + Adjoint(S) presently as mentioned above :). )