Performant conversion of Symmetric(SparseMatrixCSC) to SparseMatrixCSC

Internals & Design
1

I would like to bring the following improvements to stdlib/SparseArrays.

julia> n = 10000
julia> Random.seed!(0); A = sprandn(n, n, 0.01); b = randn(n);  nnz(A)
1000598

julia> SA = Symmetric(A);
julia> @btime B = SparseMatrixCSC(SA);
  3.049 s (35 allocations: 40.03 MiB)

julia> using SparseWrappers
julia> @btime B = SparseMatrixCSC(SA);
  10.060 ms (7 allocations: 15.32 MiB)

I would probably feel better if somebody welcomed this work. The implementation exists in my package SparseWrappers but needs to be converted into a PR for julia.

3 Likes
2

stdlib/SparseArrays definitely shares the overall goal here. It would be great to improve the interoperability of the special matrix types and I agree that avoid nested special matrix types when possible is most likely part of the solution. We might also need other overall design changes. The current state is the result of a relatively organic process and not a grand design.

However, I disagree in some of the proposals in SparseWrappers and some of them go against conclusions in existing Julia issues. E.g we have already discussed your issue 2 and it’s very useful that Hermitian works as a Hermitian view that ignores the non-real part of the diagonal. It’s also unclear to me why Symmetric(Symmetric(A, :U), :L) would be Diagonal. Should it instead have been UpperTriangular(LowerTriangular(A))? I agree that the latter should be Diagonal.

3

Thanks for looking at the SparseWrsappers. Unfortunately the README.md is not maintained at its best.

I see, that this is comfortable sometimes. But I also think, it is not a good idea, because it makes some reasoning about combinations of wrappers unnecessarily complex.
For example (Symmetric(Hermitian(A))) == Symmetric(A) is true if isreal(diag(A)) but not in general. So with the proposed restriction the simplification is possible, otherwise not.

Other question: does it make sense to have Symmetric{Complex} or Adjoint{Real}. In all my mathematical experience I always saw either (Real, Symmetric, Transpose) or (Complex,Hermitian,Adjoint), but never mixtures like Transpose(Ajoint(A)) == conj.(A)…

Symmetric(Symmetric(A, :U), :L) would be Diagonal

That is a mistake in README.md only, sorry for that. It should be
Symmetric(Symmetric(A, :U), :L) == Symmetric(A, :U) and not throw an exception ArgumentError: Cannot construct Symmetric; uplo doesn't match as it is now.