Modify check_point_clouds to allow for sparse input

[ulupo]

Types of changes

• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to change)

Description
Following the addition of FlagserPersistence in #339 and the added support for list input in the simplicial homology transformers, it is a shame to keep excluding list of sparse matrix input, especially since the performance gains can be considerable for situations in which several edges are infinitely-valued (as @lewtun found out recently when using ripser). This PR modifies the check_point_clouds validation function to accept list of sparse matrix input, while still performing important checks. Furthermore, some docstrings have been edited for clarity/consistency.

Any other comments?
I tried to extend the changes to SparseRipsPersistence but got a mysterious error thrown when trying it on a single 2x2 sparse matrix with ones in the off-diagonal entries. @MonkeyBreaker @gtauzin, do you know why this might be happening? UPDATE: It is now clear that SparseRipsPersistence is not meant to receive sparse input yet.

Checklist

• I have read the guidelines for contributing.
• My code follows the code style of this project. I used flake8 to check my Python changes.
• My change requires a change to the documentation.
• I have updated the documentation accordingly.
• I have added tests to cover my changes.
• All new and existing tests passed. I used pytest to check this on Python tests.

[MonkeyBreaker]

Could you provide MWE code producing the error to be easier to profile ?

Thank you !

[ulupo]

@MonkeyBreaker sure!

when trying it on a single 2x2 sparse matrix with ones in the off-diagonal entries.

By this I meant:

from scipy.sparse import csr_matrix
from gtda.homology import SparseRipsPersistence

sparse_1 = csr_matrix(([1, 1], ([0, 1], [1, 0])))
SparseRipsPersistence(metric='precomputed')._gudhi_diagram(sparse_1)

which returns

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-2-edb7b0d0cfa2> in <module>
      3 
      4 sparse_1 = csr_matrix(([1, 1], ([0, 1], [1, 0])))
----> 5 SparseRipsPersistence(metric='precomputed')._gudhi_diagram(sparse_1)

~/Workspace/giotto-tda_ulupo/gtda/homology/simplicial.py in _gudhi_diagram(self, X)
    373         sparse_rips_complex = SparseRipsComplex(
    374             distance_matrix=Xdgms, max_edge_length=self.max_edge_length,
--> 375             sparse=self.epsilon)
    376         simplex_tree = sparse_rips_complex.create_simplex_tree(
    377             max_dimension=max(self._homology_dimensions) + 1)

~/Workspace/giotto-tda_ulupo/gtda/externals/python/rips_complex_interface.py in __init__(self, points, distance_matrix, max_edge_length, sparse)
     80             self.thisref.init_matrix_sparse(distance_matrix,
     81                                             max_edge_length,
---> 82                                             sparse)
     83         else:
     84             if points is None:

ValueError: allocator<T>::allocate(size_t n) 'n' exceeds maximum supported size

[ulupo]

@MonkeyBreaker just to be clear, I am not aware whether SparseRipsPersistence should be supporting sparse matrix input, even in principle!

[MonkeyBreaker]

@ulupo I just give it a look, self.thisref.init_matrix_sparse(distance_matrix, max_edge_length, sparse) expects this type of arguments:

const std::vector<std::vector<double>>& matrix, double threshold, double epsilon

But looking the values passed to the arguments:

distance matrix   (0, 1)        1.0
                  (1, 0)        1.0
max_edge_length inf
sparse 0.1

distance matrix isn't a vector of vector, hence the error that occurs.

[ulupo]

@MonkeyBreaker OK. So, quite simply, we just don't support sparse matrices in SparseRipsPersistence, correct?

[MonkeyBreaker]

If sparse format means this:

(0, 1)        1.0
(1, 0)        1.0

Then no. But what we could support is:

((0, 1), (1.0))
((1, 0), (1.0))

[MonkeyBreaker]

Also, the documentation of the Class SparseRipsComplex says:

"""
:param distance_matrix: A distance matrix (full square or lower
              triangular).
"""

I just don't understand the 1.0 at the end of the lines. Because ((0, 1), (1,0)) could also be a valid distance matrix.

[ulupo]

I just mean a sparse data structure in general. The ripser bindings support many scipy sparse formats (CRS, COO, ..., see https://docs.scipy.org/doc/scipy/reference/sparse.html), while for SparseRipsPersistence we seem to restrict the Python user to dense numpy arrays.

[ulupo]

I just don't understand the 1.0 at the end of the lines. Because ((0, 1), (1,0)) could also be a valid distance matrix.

I'm interpreting your question as being about the meaning of

(0, 1)        1.0
(1, 0)        1.0

What you see there is just a way of printing a CSR matrix in python. I means that entry [0, 1] (row 0, column 1) has value 1.0, and entry [1, 0] has value 1.0. All other entries ([0, 0] and [1, 1] in this case) are not set (this is what being sparse allows you to do) and hence are automatically assumed to be zero (or infinity or missing, according to the context).

[MonkeyBreaker]

Correct me if I'm wrong, but before calling ripser there is some preprocessing. If I modify the example like this:

from scipy.sparse import csr_matrix
from gtda.homology import SparseRipsPersistence
from sklearn.metrics.pairwise import pairwise_distances

sparse_1 = csr_matrix(([1, 1], ([0, 1], [1, 0])))
sparse_2 = pairwise_distances(sparse_1, metric='euclidean')
SparseRipsPersistence(metric='precomputed')._gudhi_diagram(sparse_2)

I get a new error related to the python code in simplicial.py:

Traceback (most recent call last):
  File "test.py", line 7, in <module>
    SparseRipsPersistence(metric='precomputed')._gudhi_diagram(sparse_2)
  File "/home/julian/workspace/TWS/giotto_learn_github/gtda/homology/simplicial.py", line 375, in _gudhi_diagram
    max_dimension=max(self._homology_dimensions) + 1)
AttributeError: 'SparseRipsPersistence' object has no attribute '_homology_dimensions'

[ulupo]

@MonkeyBreaker with your example, you are testing a different mathematical situation. You are interpreting sparse_1 as a list of 2D vectors, and then computing all pairwise distances between them, which will be a numpy array indeed. This was not the intended usage in my example. I intended to pass sparse_1 directly as the matrix of pairwise distances, and wanted this to be sparse. My initial example works with no problems when SparseRipsPersistence is exchanged with VietorisRipsPersistence or FlagserPersistence .

[ulupo]

I'm not sure how we handle sparse matrices in the ripser interface, we'd have to check. I don't know if we do something very trivial like converting them to dense numpy array format and feeding the resulting arrays to the ripser bindings, or we exploit the sparse structure in a more clever way to avoid expensive conversions.

[MonkeyBreaker]

In ripser.py, before calling the c++ code, the input is formatted to a specific format.

The same with FlagserPersistence, maybe we should update rips_complex_interface.py to allow more input formats and then convert them to the format expected by the implementation ?
Maybe an input from @gtauzin could help us in the right direction ?

[MonkeyBreaker]

@MonkeyBreaker with your example, you are testing a different mathematical situation. You are interpreting sparse_1 as a list of 2D vectors, and then computing all pairwise distances between them, which will be a numpy array indeed. This was not the intended usage in my example. I intended to pass sparse_1 directly as the matrix of pairwise distances, and wanted this to be sparse. My initial example works with no problems when SparseRipsPersistence is exchanged with VietorisRipsPersistence or FlagserPersistence .

Sorry for that, I just "made it work", I was looking at the reason why it didn't work.

[ulupo]

I believe this

giotto-tda/gtda/externals/python/ripser_interface.py

Line 20 in 188b675

def DRFDMSparse(I, J, V, N, maxHomDim, thresh=-1, coeff=2, do_cocycles=1):

is the function responsible for passing sparse input to C++ ripser. The input must be in COO format, so any other sparse format is converted to COO before calling this function.

[MonkeyBreaker]

You're right, from the C++ code:

//Initialize from COO format
	sparse_distance_matrix(int* I, int* J, float* V, int NEdges, int N, float threshold) :
							neighbors(N), vertex_births(N, 0) {
		int i, j;
		value_t val;
		for (int idx = 0; idx < NEdges; idx++) {
			i = I[idx];
			j = J[idx];
			val = V[idx];
			if (i < j && val <= threshold) {
				neighbors[i].push_back(std::make_pair(val, j));
				neighbors[j].push_back(std::make_pair(val, i));
			}
			else if (i == j) {
				vertex_births[i] = val;
			}
		}
	}

But in this case we can conclude that at the moment, SparseRipsPersistence cannot handle sparse data structure input.

We could before calling the c++ code in SparseRipsPersistence transform the data ?

[ulupo]

@MonkeyBreaker I need to think about this. At the moment, I think that if using sparse formats only gives a performance penalty because of conversions secretly performed by us, we should not fool the user into thinking that it might be a good idea to use these formats. The situation is different for VietorisRipsPersistence and FlagserPersistence where there can be large rewards for doing so.

[ulupo]

Reason for this: Sparse matrices can easily end up being non-square, but that is OK in all current applications for which distance_matrices might be set to True, i.e. the persistent homology transformers in gtda.homology.simplicial. Missing non-diagonal entries are consistently interpreted as np.inf in those cases.

[gtauzin]

@ulupo As I mentioned earlier, I believe the check_point_clouds function should not be the one checking adjacency matrices of graphs. A long time ago a check_graphs function was started but never finished. I suggest the logic for checking adjacency matrices is offloaded there. Also, all preprocessing transformers in gtda.graphs could make use of it.

[ulupo]

@gtauzin in principle, I agree with your suggestion. In practice, it is a slightly complicated one to implement: check_point_clouds is already being used to check adjacency matrices of graphs, because it is already being used by VietorisRipsPersistence and SparseRipsPersistence for instance, when metric='precomputed'. Presumably, the function would have to be rewritten substantially (and split into two) even in the current state. The current changes are simply minimal with respect to the existing (debatable) decisions.

What do you think about making this the subject of a separate PR?

[ulupo]

Another smaller note: if #422 is merged, we won't need a check function for the graph transformers. This is because in fact only one existing transformer in gtda.graphs receives graphs as inputs, namely GraphGeodesicDistance. But, as per #422, a thorough scipy check_graph-type function is run automatically on each graph before computing geodesic distances, so it seems to me that any batch-level check_graph function would just end up doing duplicate work in that specific case.

For the future, there remains of course a general case for having check_graph.

[gtauzin]

@ulupo I think it is quite good! I have also a general comment. In check_array, if the input is a list but the shape of its element allow it to be formatted as a numpy array, it returns a numpy array. I think this would be good to have as well to enhance performance.

[CLAassistant]

Thank you for your submission! We really appreciate it. Like many open source projects, we ask that you all sign our Contributor License Agreement before we can accept your contribution.
1 out of 2 committers have signed the CLA.

✅ ulupo
❌ Umberto

---

Umberto seems not to be a GitHub user. You need a GitHub account to be able to sign the CLA. If you have already a GitHub account, please add the email address used for this commit to your account.
_{You have signed the CLA already but the status is still pending? Let us recheck it.}

[ulupo]

@gtauzin I have now addressed two of the points you raised (conversion to an overall array when possible, removal of dimensionality warnings). Concerning copies, this is the current situation:

When copy=False and the input is a "regular" list – of, say, float arrays – the return value of check_point_clouds is a new list "typically" pointing to the same objects as the input list. So e.g. doing

X_list = [np.ones((3,4)), np.zeros((3,3))]
X_arr = check_point_clouds(X_list, copy=False)
X_arr[0][0, 0] = 2

would also modify the top-left entry of the first array in X_list. Notice that this is not true (it cannot be, as far as I can tell) when the arrays in X_list have the same shapes because then the addition made in 1b110ea basically triggers a copy regardless of whether copy was set to True or False.

So in a sense we are only "half" copying in these scenarios (and the copying we do is a "light" one, i.e. a new list of pointers, with little new memory allocation).

The reason I think it's not so straightforward to do better is as follows: One could think that the new behaviour should be something like this:

1. When copy is set to True, keep the current behaviour;
2. When copy is set to False, the loop could be as follows:

        has_check_failed = False
        messages = []
        for x in X:
            try:
                _check_array_mod(x, **kwargs_, ensure_2d=True)
                if distance_matrices:
                    if not x.shape[0] == x.shape[1]:
                        raise ValueError(
                            f"All arrays must be square: {x.shape[0]} rows "
                            f"and {x.shape[1]} columns found in this array.")
            except ValueError as e:
                has_check_failed = True
                messages.append(f"Entry {i}:\n{e}")
        if has_check_failed:
            raise ValueError(
                "The following errors were raised by the inputs:\n\n" +
                "\n\n".join(messages))
        else:
            Xnew = X

i.e. we just run _check_array_mod (which is virtually the same as check_array) on each entry, and if no errors are thrown we simply make Xnew an alias of the original data.

The problem with 2 I think is that not always does passing copy=False in sklearn's check_array mean that no copies are triggered. Sometimes they can be triggered by data type conversions. In these cases, we would ignore the copies which are presumably triggered by check_array for a good reason. What do you think?

[ulupo]

@gtauzin a particular case of my point above is the following: suppose (God forbid!) that a user passes a nested 3-level list instead of a list of 2D arrays. Then, in the current design, thanks to check_array's behaviour we would return a new list of 2D arrays which is important for downstream tasks. However, if the alternative in the above code block is used instead, we would just return the original nested list.