RFC: sparse: establish an official policy for duplicate elements

[perimosocordiae]

Background
As an admission to performance and ease of construction, most of our sparse matrix formats admit duplicate entries. For example:

>>> foo = coo_matrix(([4,5], ([0,0],[1,1])), shape=(2,2))
>>> foo.nnz
2
>>>foo.A
array([[0, 9],
       [0, 0]])

Unfortunately, this causes a lot of problems for sparse matrix operations (#4409, #5394, #5806), and causes confusion regarding the true meaning of nnz, even ignoring the issue of explicit zeros (#3343).

Most sparse matrix formats have a method sum_duplicates() which operates in-place to canonicalize the internal storage, but it's unclear whether other methods are allowed to call this without first making a copy (see #5741 (comment)).

Proposal
I think that allowing duplicate entries in internal sparse matrix representations was an API mistake, but now that we allow it, it's difficult to disallow without breaking lots of existing user code. Therefore, I propose that we make it a policy that:

1. Duplicate entries are not preserved. That is, it's okay to canonicalize in-place.
2. Whenever a method other than sum_duplicates() triggers in-place canonicalization, a SparseEfficiencyWarning is thrown, to alert the user that something potentially unexpected is going on.
3. The presence/lack of duplicate entries is remembered with a boolean flag, which we will document and encourage users to toggle if they manually modify a sparse matrix's internal members.

If this gets traction here, I'll send it along to the Scipy-dev mailing list as well. Thoughts?

[rgommers]

+10 for documenting this kind of thing clearly.

(1) and (2) of the proposal seem fine to me, I'm not sure about (3). I was trying to find a place in our docs where we explain the structure of the matrix formats and if/how to manually modify data. But there's very little in the tutorial or refguide - the best description is probably this one: http://www.scipy-lectures.org/advanced/scipy_sparse/index.html

[pv]

This is blocking: (at least) gh-5394

Agreed on that it would have been better to not permit duplicates. I understand they are useful at matrix assembly phase, but I think this would be better handled by a specialized matrix format rather than forcing it everywhere. It is unclear to me if the present support for duplicates has good uses; I guess most either use lil/dok (which are also inefficient due to Python overhead), or accumulate coo-format arrays manually.

I would go allowing canonicalization at any time (1).

[pv]

The other RFC to decide is gh-3343

[ewmoore]

Adding a flag that user is suppose to toggle is just asking for subtle broken programs. +1 to canonicalization at any time (1).

[perimosocordiae]

Yeah, it seems like (3) could be better handled by adding some documentation, rather than more complexity in the code.

It looks like (1) is agreed to be a good choice. We should figure out a good place to document it, then inform the mailing list.

I think we can implement (2) pretty cheaply by adding a warn=False keyword argument to sum_duplicates(), and force internal callers to pass warn=True. It's not a particularly elegant solution, but I like it better than using inspect to peek at the caller's stack frame.

[hpaulj]

I first used sparse matrices in MATLAB for finite element problems (in the 1990s). Coming from that background I regard the summation of duplicates as a useful, if not essential, feature. I just used it in a Stackoverflow answer, http://stackoverflow.com/a/37040831/901925

The only documentation that I'm aware of is on the coo_matrix page

By default when converting to CSR or CSC format, duplicate (i,j) entries will be summed together. This facilitates efficient construction of finite element matrices and the like. (see example)

I haven't studied the compile code much, so don't know exactly where and how the summation is implemented. My observation is that the .row and .col attributes of a coo format might actually be the arrays provided as input in the (data, (i, j)) pattern (if they are arrays of the correct dtype). So at least in some cases duplicates are present simply because inputs have not been processed or checked.

I haven't observed whether duplicates are allowed with the other formats. bmat uses coo matrices as intermediates, but doesn't allow overlapping blocks. It is obvious how a lil could have duplicates. dok is a dictionary, so it can't have duplicates, can it?

I haven't taken time yet to look at the cited issues where this seems to be a problem. But I just wanted to express support for the existing practice. Documentation can always be improved without affecting legacy practices.

[jni]

Hi all,

I'd like to contribute a big -1 to this discussion. #5394 already sent me on a bug hunt and reduced efficiency in my own code (since I have to explicitly make a copy of a coo_matrix before converting to CSR). As mentioned in #6603, in that PR, canonization occurs even when copy=True is passed to tocsr().

One of the motivating examples is #4409, which in my reading points to an application (csgraph) that absolutely cannot work when sparse is canonizing arrays left, right, and center, namely, the existence of graph edges with weight 0 (as opposed to non-edges).

Another, #5806, handled the non-canonical case cleanly.

I can't count right now the number of places I use the internal structure of COO and CSR. It's certainly used within scikit-image. In my own code I've implemented a CSR matrix that can get expanded with new rows (here). I'd expect arbitrary canonization to wreak havoc with it. Along with any number of user subclasses.

My suggestion is not to mess with users' objects unless explicitly allowed by the user. For example, you could add canonize=False to all the places where canonization would be useful. If False, make a copy of the original array if necessary. This mirrors the common pattern of having an output=input kwarg for in-place operations.

In short, this isn't a small case of "we might break some users' code but it's ok since it's undocumented". This has massive downstream ramifications.

[jni]

To be clear, my argument is -1 for the policy being in-place canonization by default, not for establishing a policy. =P

[perimosocordiae]

@jni: Note that "canonicalization" in this context refers only to duplicate entries, not explicit zeros. The discussion on explicit zeros (gh-3343) seems to have converged on preservation, specifically for the csgraph example you mentioned. If there are operations which don't preserve explicit zeros (other than, say, eliminate_zeros()) those should be reported as bugs.

That said, there is also a case for using duplicate entries in a sparse matrix to represent a multigraph. I personally think that this goes a little too far, and that a specialized data structure would be more appropriate. I'd like to hear what others think on the topic.

My ideal is to allow duplicate entries only during construction via __init__: once the spmatrix subclass is initialized, all of its internal members should be duplicate-free and unchanged by normal operations. This would allow us to get rid of the numerous self.has_canonical_format checks and calls to sum_duplicates() which currently litter the codebase, eliminating all the problems with in-place canonicalization. It does mean that the arguments passed to the __init__ method may end up getting copied in some cases, but judicious use of kwargs should allow the user to say "I know what I'm doing, just use these arrays directly" with the understanding that duplicate entries will cause incorrect results for some operations.

I'm interested to find out which use-cases would break under this "canonicalize-at-init" proposal. I can't think of any at the moment, but I'm sure some exist!

[pv]

Canonicalization on init may have the same issue, because then with the (default) copy=False, `__init__` would modify the provided arrays. With copy=True, the data arrays could contain different data, which also may break assumed invariants. . The question here is deciding where the line of "spacebar heating" goes. Strict backward compatibility doesn't really allow canonicalization at any stage AFAICS so the internal format would be stuck with duplicates etc., with the costs what that entails.