Improve performance of issym and ishermitian for sparse matrices

[KristofferC]

This implements what @toivoh mentioned here JuliaLang/LinearAlgebra.jl#217

My testing shows that it is about 3 times faster for hermitian matrices than the current one and for non hermitian it usually terminates right away so it gains a lot there. It also uses less memory.

The logic in this function is of course more complicated than the current one so it is important that it is tested properly. I prodded it with as many different cases I could come up with and it works for all those but maybe there should be more tests in base?

Edit:

Some small benchmarks:

sps = [0.0001, 0.001]
ns = [10^3, 10^4, 10^5]
n_iters = 10

for n in ns
    t = 0.0
    for sp in sps
        for i = 1:n_iters
            A = sprand(n, n, sp)
            A = A + A'
            t += @elapsed issym(A)
        end
        println("n: $n, sparsity: $sp, avg time: $(t/n_iters)")
    end
end

Results

# Master
n: 1000, sparsity: 0.0001, avg time: 1.69745e-5
n: 1000, sparsity: 0.001, avg time: 5.023520000000001e-5
n: 10000, sparsity: 0.0001, avg time: 0.0004951787999999999
n: 10000, sparsity: 0.001, avg time: 0.0035960212999999997
n: 100000, sparsity: 0.0001, avg time: 0.05512627459999999
n: 100000, sparsity: 0.001, avg time: 0.5852790912

# This branch
n: 1000, sparsity: 0.0001, avg time: 3.3673e-6
n: 1000, sparsity: 0.001, avg time: 1.22526e-5
n: 10000, sparsity: 0.0001, avg time: 8.872999999999999e-5
n: 10000, sparsity: 0.001, avg time: 0.0009253101999999999
n: 100000, sparsity: 0.0001, avg time: 0.020136009700000002
n: 100000, sparsity: 0.001, avg time: 0.20595624689999997

[tkelman]

maybe there should be more tests in base?

Always. Especially with sparse stuff, you should absolutely add a bunch of new test cases with this.

[KristofferC]

I updated with more test cases to issparse and ishermitian.

[KristofferC]

Build timed out on x64

[KristofferC]

Is there any consensus how functions should behave on sparse matrices where explicit zeros have been inserted into the structure?

[toivoh]

I think what makes most sense usually is to preserve them (others will have to correct me if I'm wrong). Our in what respect?

[KristofferC]

In the respect of correctness. If you allow explicit zeros, a symmetric sparse matrix does not necessarily have a symmetric sparse structure.

For example, with my function in the PR:

A = speye(10,10)
A2 = speye(10,10)
A[8,9] = 1.0
A.nzval[9] = 0
A == A2 # True
ishermitian(A) # True
ishermitian(A2) # False(!)

In order for this to work, you need insert some checks for zero.

[tkelman]

For sparse matrices, "is numerically symmetric" is a distinct property from "is structurally symmetric." It's possible to construct examples of all 4 combinations. Numerical symmetry should treat explicit zeros as equal to implicit zeros in the other triangle, structural symmetry should treat explicit zeros as a structurally stored location.

[ScottPJones]

@KristofferC The cases I've tried will not represent them, and on the other hand, if you set a location to None, it preserves it... 👎
Another issue, the documentation states the following, which doesn't seem to be true:

In some applications, it is convenient to store explicit zero values in a SparseMatrixCSC. These are accepted by functions in Base (but there is no guarantee that they will be preserved in mutating operations). Such explicitly stored zeros are treated as structural nonzeros by many routines. The nnz() function returns the number of elements explicitly stored in the sparse data structure, including structural nonzeros. In order to count the exact number of actual values that are nonzero, use countnz(), which inspects every stored element of a sparse matrix.

julia> z = sparse([1,2], [1,2], ["hello", 1.5])
2x2 sparse matrix with 2 Any entries:
    [1, 1]  =  "hello"
    [2, 2]  =  1.5

julia> z[1,2] = None
Union()

julia> z
2x2 sparse matrix with 3 Any entries:
    [1, 1]  =  "hello"
    [1, 2]  =  Union()
    [2, 2]  =  1.5

julia> nnz(z)
3

julia> z[1,2] = 0
0

julia> nnz(z)
2

julia> z
2x2 sparse matrix with 2 Any entries:
    [1, 1]  =  "hello"
    [2, 2]  =  1.5

This is rather unfortunate, because it means these are not suitable for database work...
(if "missing" data was always null, and null was treated as 0 for mathematical operations, you'd get the same mathematical results, and the ability to store 0s, as per the documentation)

Am I missing something here?

[tkelman]

setindex with a rhs of 0 does eliminate the value from the nonzero structure right now (cf. "there is no guarantee that they will be preserved in mutating operations"), which is not always what you want, I agree.

[ScottPJones]

OK, when I read that about mutating operations, I didn't think it meant a plain set... especially after reading:

it is convenient to store explicit zero values

[KristofferC]

Please @ScottPJones, if you start your post with "Another issue" it should probably be in "Another issue".

@tkelman So as it is, the function in this PR is bugged because it fails on no structurally symmetric, symmetric matrices. I will fix and update.

[ScottPJones]

@KristofferC Well, it seemed perfectly on point about your question... the docs say one thing, but the code seems to do something else...

Is there any consensus how functions should behave on sparse matrices where explicit zeros have been inserted into the structure?

[KristofferC]

@ScottPJones Yes, I should have said that it was more concretely about the function in this PR. A general discussion about dealing with explicit zeros in sparse matrices should probably be held in a separate issue.

Anyway, I fixed the problem now so that

A = speye(10,10)
A2 = speye(10,10)
A[8,9] = 1.0
A.nzval[9] = 0
A == A2 # True
ishermitian(A) # True
ishermitian(A2) # True

Also rebased.

[KristofferC]

Test fail seem to be unrelated?

ERROR (unhandled task failure): readcb: connection reset by peer (ECONNRESET)
 in yieldto at ./task.jl:21
 in wait at ./task.jl:309
 in wait at ./task.jl:225
 in wait_full at ./multi.jl:572
 in remotecall_fetch at multi.jl:672
 in remotecall_fetch at multi.jl:677
 in anonymous at task.jl:1392
    From worker 4:       * examples             in  77.77 seconds
    From worker 5:       * subarray             in 347.51 seconds
exception on 7: ERROR: LoadError: ProcessExitedException()
 in worker_from_id at ./multi.jl:326
 in call_on_owner at ./multi.jl:719
 in wait at ./multi.jl:724
 in SharedArray at sharedarray.jl:73
 in include at ./boot.jl:252
 in runtests at /tmp/julia/share/julia/test/testdefs.jl:197
 in anonymous at multi.jl:838
 in run_work_thunk at multi.jl:589
 in anonymous at task.jl:838
while loading /tmp/julia/share/julia/test/parallel.jl, in expression starting on line 20
ERROR: LoadError: ProcessExitedException()
 in anonymous at task.jl:1394
while loading /tmp/julia/share/julia/test/runtests.jl, in expression starting on line 5

[ViralBShah]

@ScottPJones I request sticking to the discipline of discussing relevant matters in an issue, or else it becomes inconvenient for everyone else. Like I had suggested earlier, the issue of zeros in sparse matrices for non-numerical data is best discussed on julia-users, or through a PR, or perhaps even a separate issue.

[ScottPJones]

@ViralBShah I was simply responding to this, and bringing up the point that the code and documentation didn't seem to agree... That did seem relevant to the question.

Is there any consensus how functions should behave on sparse matrices where explicit zeros have been inserted into the structure?

[ViralBShah]

Fair enough. The documentation should be updated too.

BTW, issym by default should be about numerical symmetry, which means stored zeros should be ignored. We should certainly have optional arguments for structural symmetry.

[KristofferC]

What is the actual definition of structural symmetry?

From http://docs.oracle.com/cd/E19205-01/819-5268/plug_sparse.html:

"A structurally symmetric sparse matrix has nonzero values with the property that if a(i, j) not equal 0, then a(j, i) not equal 0 for all i and j."

This seems to exclude cases like

[1 1]
[Ø 0]

where Ø is a stored zero.

But from @tkelman's comment it seems like this should be regarded as structural symmetric.

[ViralBShah]

IMO, structural symmtery ignores the numerical values.

[KristofferC]

Ok, I agree. There are two ways to implement the isstructsym, one way is just to manually remove the numerical checks, like here: https://gist.github.com/KristofferC/5b02b179208f72bbf591

Another way is to make ishermsym take a bool and do the numerical checks based on the bool. This leads to less code duplication at the cost of a bit larger ishermsym.

[toivoh]

But you would also like to have the option to check structural and numeric symmetry at the same time I guess.

[mauro3]

Whether to store or not store zeros has been discussed before, but without definitive conclusion. The most discussion happened here #9928, where most wanted to keep zeros.
See also JuliaLang/LinearAlgebra.jl#60 and #9906

[KristofferC]

If we are going to have the store zeros or not discussion here, I would say I am for keeping the sparse structure when you can.

Stuff like this feels silly:

julia> A = sprand(10^5,10^5, 0.001);

julia> A[20, 20] = 1.0
1.0

julia> 1 / ( (@elapsed A[20, 20] = 0.3) / (@elapsed A[20, 20] = 0.0))
32.11705573550765

Regarding the PR, I don't think we should add too much. The symmetric structure checking could easily be added in a separate PR.

[tkelman]

I've been thinking that setindex with a zero rhs should introduce an explicit zero when it's assigning to a location that is already stored in the matrix, but that'll have to be argued about in a different PR.

I'm not sure I'd take sunperf's documentation as the definitive truth here, I don't know what their policy is on explicit zeros there. MKL might be slightly more informative?

Structural symmetry shouldn't look at the contents of nzval at all.

[KristofferC]

Updated with some actual benchmarks.

[andreasnoack]

@KristofferC Is this ready? LGTM.

[KristofferC]

Yeah, should be ready for merge.

[andreasnoack]

Great. Thanks for the work.