Chapel: Support LinearAlgebra.isSymmetric for sparse matrices

[buddha314]

/*
 chpl -I/usr/local/Cellar/openblas/0.2.20/include -L/usr/local/Cellar/openblas/0.2.20/lib -lblas issym.chpl
 */
use 
    LinearAlgebra;

var nv: int = 8,
    D: domain(2) = {1..nv, 1..nv},
    SD: sparse subdomain(D),
    W: [SD] real;

SD += (1,2); W[1,2] = 1;
SD += (1,3); W[1,3] = 1;
SD += (1,4); W[1,4] = 1;
SD += (2,2); W[2,2] = 1;
SD += (2,4); W[2,4] = 1;
SD += (3,4); W[3,4] = 1;
SD += (4,5); W[4,5] = 1;
SD += (5,6); W[5,6] = 1;
SD += (6,7); W[6,7] = 1;
SD += (6,8); W[6,8] = 1;
SD += (7,8); W[7,8] = 1;

writeln(isSymmetric(W));

And then...

./issym 
true

Wha happen...?

[cassella]

isSymmetric() has an optimization to avoid considering both A[i,j] != A[j,i] and A[j,i] != A[i,j]. That's causing trouble for your sparse domain -- all the off-diagonal elements have indices (i,j) with i < j so it never gets to compare to (j,i):

/* Return `true` if matrix is symmetric */
proc isSymmetric(A: [?D]) : bool {
  if D.rank != 2 then
    compilerError("Rank size is not 2");
  for (i, j) in D {
    if i > j {
      if A[i, j] != A[j, i] then return false;
    }
  }
  return true;
}

(For the same reason, it reports that this dense array is symmetric:

var A : [1..2, 1..10] real;
A[2,7] = 3.9;

But maybe one is expected to know better than to ask if non-square matrices are symmetric.)

The transpose of your array should be reported as isNotSymmetric.

[bradcray]

@cassella: One implication of your observation is that the optimization doing only one triangle is incorrect, or at least too naive (it might be correct if it was combined with a check that the number of indices above and below the diagonal was equal?)

Is another implication that isSymmetric() should check if the matrix is square? Or is it reasonable to have an isSymmetric() call only check those indices in the square prefix of the complete matrix? (and in that case, should it check that the remaining indices are 0?)

Sign me,
Forgot Linear Algebra Ways

[buddha314]

Is there an efficient way to check super-linearity of the (i,j) pairs? Or perhaps a spiral method? Is there an efficient way to do (i,j) =? (j,i), (j+1,i+i) =? (i+1, j+1) That way you'd be alternating super and sub-diagonal postulates.

[cassella]

@bradcray Unless I'm missing something, I don't think you can depend on counting indices. Consider replacing all of the indices in the original example with just these two, (where the second's array value is set to 0, not 1):

SD += (1,2); W[1,2] = 1;
SD += (3,1); W[3,1] = 0;

I don't think I can contribute to the other parts of the discussion,

[bradcray]

@casella: I wasn't suggesting merely counting indices, I was suggesting combining the current approach (comparing only one half of the triangle) with counting indices in each triangle. I believe if all lower triangular indices have a matching upper triangular index _and_ the number of indices in each triangle is the same, then the matrix could be considered symmetric. Sorry if I wasn't clear about that.

That said, I'm also skeptical that, at least with a vanilla COO or CSR/CSC format, this would be any cheaper than just doing the full O(n**2) traversal.

[casella]

I guess the right person is @cassella with 2 's'

[bradcray]

Oh wait, should your example have included a third line:

SD += (2,1);  W[2,1] = 1;

such that the (2,1) and (1,2) entries were symmetric, but the number of indices differed, but the difference was a 0 element, and therefore not a difference in terms of matrix symmetry? In that case, I understand what you were trying to say better now and retract my previous answer that counting indices is sufficient. It would be sufficient for a quick exit in some cases I think (if all indices in a triangle were symmetric in the sense of "both triangles having that index and the arrays having the same value stored in it) but not others (the other triangle could have more indices but they could have 0 array elements and it could still be symmetric. So I think you were more right than I thought. In any case, I think we all agree that the current implementation is incorrect / buggy. And even with my renewed understanding, I'm not implying that I think counting the indices is a good / efficient approach.

@casella / @cassella: Sorry for the mix-up. Not that it's an excuse, but learning Italian ruined me for double consonants.

[ben-albrecht]

In #8088, I have addressed the non-square bug, and have updated the documentation to be more explicit that none of the functions outside of LinearAlgebra.Sparse are expected to support sparse formats, unless explicitly noted in their documentation.

Beyond that, I think this issue should remain open until isSymmetric supports CSR arrays.

[bradcray]

Oh, I didn't realize that the current isSymmetric was only intended to work with dense arrays... In the interest of closing this issue and/or preventing others from stepping into this hole, I'd suggest either implementing the naive n**2 implementation for CSR (replacing the i < j with i != j, I think the same code would work, and faster versions could come later) or else have the dense-only version include a where-clause or type check to make sure it's not called with a sparse domain?