Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Inconsistencies with transpose and Diagonal * Sparse #420

Closed
georgemarrows opened this issue Apr 18, 2017 · 6 comments
Closed

Inconsistencies with transpose and Diagonal * Sparse #420

georgemarrows opened this issue Apr 18, 2017 · 6 comments
Labels
Hacktoberfest Good for Hacktoberfest participants sparse Sparse arrays

Comments

@georgemarrows
Copy link

This bug report is split out from @fredrikekre's https://github.com/JuliaLang/julia/issues/21286#issuecomment-292077402:

Here are some other faulty things with Diagonal/SparseMatrixCSC I found looking at JuliaLang/julia#20367 but then forgot about:

D = Diagonal(ones(4))
S = sprand(4, 4, 0.4)

D\S  # returns a "full" SparseMatrixCSC

S'*D # errors (works if eltype(D) = Complex)
D*S' # errors (works if eltype(D) = Complex)

S*D' # returns an Array
D'*S # interestingly enough this works

I (George) have checked and the issues are still present at

julia> versioninfo()
Julia Version 0.6.0-pre.beta.203
Commit b0c708424b* (2017-04-18 17:33 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin16.4.0)
  CPU: Intel(R) Core(TM) i7-3740QM CPU @ 2.70GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, ivybridge)

This obviously is related to #136 and JuliaLang/julia#21431

@georgemarrows
Copy link
Author

I'll have a look at this in the next day or two.

@mfalt
Copy link
Contributor

mfalt commented Apr 19, 2017

I'm not sure if it is relevant to this issue, similar problems arise for the Symmetric type of sparse matrices, ex:

julia> s = sprandn(3,3,0.5);
julia> ss = Symmetric(s'*s) #Everything is still sparse
3×3 Symmetric{Float64,SparseMatrixCSC{Float64,Int64}}:
 1.78954    0.0318829  0.0
 0.0318829  0.166128   0.0
 0.0        0.0        0.0
julia> ss*2  # Symmetric{Float64,SparseMatrixCSC{Float64,Int64}} Works fine
julia> 2*ss  # Array{Float64,2} Not sparse and lost Symmetric tag 
julia> ss+ss # Array{Float64,2} Same problem
julia> ss*ss # Array{Float64,2} Same problem

and as noted above (but with additional missing Symmetry tag)

ss'*Diagonal(randn(3)) # Array{Float64,2} Not sparse no Symmetry tag
Diagonal(randn(3))'*ss # Same

And I don't understand the reasoning behind

julia> ss[1,2] = 2
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix

which causes ArgumentError: Cannot set a non-diagonal index in a symmetric matrix on calls like

ss*Diagonal(randn(3))
Diagonal(randn(3))*ss

I guess this is related to JuliaLang/julia#8240, but for sparse matrices.
(This was on 0.6.0-pre.beta.182)

@georgemarrows
Copy link
Author

The cause of this part of the issue

S'*D # errors (works if eltype(D) = Complex)
D*S' # errors (works if eltype(D) = Complex)

is that the A[c]_mul_B[c] instances here are provided by diagonal.jl (lines 190 and 205). Those methods expect to be able to use similar to create a new matrix and then immediately ctranspose! into it. That's not a contract that sparse matrices honour because of the checks in the specialised ctranspose! in sparsematrix.jl.

When eltype(D) = Complex, we end up in generic ctranspose!, which doesn't have the checks and all works.

So who is at fault?

  • Diagonal, because there's no "can always ctranspose! after similar" contract?
  • SparseMatrixCSC, because it breaks that contract?

@Sacha0 - as far as I can tell, the ctranspose! checks came in as part of the JuliaLang/julia#13001 rewrite you did to remove LGPL. Can/should ctranspose! be changed to "fatten up" its target matrix if necessary so it can be transposed into? Or would that hide performance costs too much?

@Sacha0
Copy link
Member

Sacha0 commented Apr 20, 2017

Hm, tricky that. Not automagically resizing the destination sparse matrix's storage seems like it breaks a worthwhile abstraction, so perhaps automagically resizing is best? Opinions from the usual sparse suspects? Best!

@ettersi
Copy link

ettersi commented May 8, 2017

For all those who need a fix right now, D*ctranspose(S) works.

@fredrikekre fredrikekre added Hacktoberfest Good for Hacktoberfest participants linear algebra and removed linear algebra labels Sep 28, 2017
@andreasnoack
Copy link
Member

This is fixed in 0.7

@KristofferC KristofferC transferred this issue from JuliaLang/julia Nov 26, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Hacktoberfest Good for Hacktoberfest participants sparse Sparse arrays
Projects
None yet
Development

No branches or pull requests

7 participants