ENH: sparse.linalg: Implement matrix_power()

[ljwolf]

Reference issue

No specific reference issue

What does this implement/fix?

This adds a scipy.sparse.linalg.matrix_power() analogue to numpy.linalg.matrix_power(). It uses the old _spmatrix.__pow__() recursive implementation, since MatrixPowerOperator() was quite a bit slower.

[ljwolf]

Whitespace issues fixed. 37 tests fail expecting _spmatrix.__pow__ to raise an exception for a sparse matrix raised to a negative power.

[ljwolf]

We decided to keep the default behaviour in both MatrixPowerOperator() and _spmatrix.__pow__, which raises an error when the user raises a matrix to a negative power.

If we allowed the user to pass a negative power, we would have to invert the input array A. That would usually create a dense matrix unless there is special structure in the input matrix. It's easy enough for an interested user to invert the matrix themselves and use abs(pow), so we'll keep the current behaviour consistent by raising an error with a negative power.

[rossbar]

Just to note - it looks like #513 introduced a recursive implementation that saves on the number of matmuls for larger exponents. It would be worthwhile to ensure there is a benchmark in place for sparse matrix power to ensure that there are no performance regressions with the updated implementation. @perimosocordiae is working on a benchmark in a separate PR!

[ljwolf]

Happy to pull in the old implementation, too. I figured that the MatrixPowerOperator() implementation would be faster than the recursive call, especially if structure is known.

[perimosocordiae]

I just opened gh-18553 with a benchmark covering this case. Once that merges, it'll be really easy to evaluate the impact of this PR's changes using dev.py bench --compare main.

[perimosocordiae]

My benchmark PR is now merged.

[ljwolf]

Yeah, this is not gonna fly. The LinearOperator() version is an order of magnitude slower than the recursive solution on my computer:

on main
              === ============== ==============

              --           N / density
              --- -----------------------------
               x   1000 / 1e-06   1000 / 0.001
              === ============== ==============
               0     49.7±3μs       52.0±6μs
               1    66.7±20μs       58.9±7μs
               2     209±20μs       268±10μs
               3     338±20μs       451±10μs
               8     516±50μs       662±50μs
               9     697±50μs       856±20μs
              === ============== ==============
in PR: 
              === ============== ==============
              --           N / density
              --- -----------------------------
               x   1000 / 1e-06   1000 / 0.001
              === ============== ==============
               0     424±10μs       428±30μs
               1     826±30μs       845±40μs
               2     932±20μs     1.05±0.06ms
               3   1.13±0.06ms     1.52±0.3ms
               8   1.76±0.05ms     2.13±0.1ms
               9    2.11±0.2ms     2.35±0.2ms
              === ============== ==============

I'll swap back to the _spmatrix.__pow__ implementation, rather than scipy.sparse.linalg.MatrixPowerOperator().

[ljwolf]

OK, this now uses the older recursive __pow__() implementation for scipy.sparse.linalg.matrix_power(), and transitions _spmatrix.__pow__ to use scipy.sparse.linalg.matrix_power(), so the benchmark is about the same as on main:

=== ============== ==============
--           N / density
--- -----------------------------
 x   1000 / 1e-06   1000 / 0.001
=== ============== ==============
 0     56.2±5μs       52.1±4μs
 1     60.2±8μs      70.8±20μs
 2     197±20μs       254±30μs
 3     351±20μs       442±30μs
 8     533±50μs      769±200μs
 9     627±80μs      905±100μs
=== ============== ==============

[perimosocordiae]

Commenting here to say that I haven't forgotten about this PR! I still need to do one more read through the changes, but in the meantime, could you update the PR description to reflect the current state of the change?

[ljwolf]

Great, the description of the PR and docstring is now fixed.

[jjerphan]

Thank you, @ljwolf!

Here are a few comments.

[jjerphan]

What is the benefit of proceeding as follow, here?

[ljwolf]

This is part of the existing matrix_power implementation, which (as I understand it) was used to reduce the number of required multiplications. This performance improvement was significant then, and seems so again.

[jjerphan]

Ah yes, that's just a recursive implementation.

We could further reduce the cost of this function using inline computations with more or less np.log2(power) iterations not to have the checks just above be rerun several times.

What do you think?

[ljwolf]

OK, I took all of @jjerphan's suggestions, and clarified the comments! let me know where to take this further!

[jjerphan]

Here is a second pass.

I just have another comment regarding the implementation of matrix_power which I think can be improved.

[jjerphan]

Ah yes, that's just a recursive implementation.

We could further reduce the cost of this function using inline computations with more or less np.log2(power) iterations not to have the checks just above be rerun several times.

What do you think?

[WarrenWeckesser]

A couple docstring issues:

• As @jjerphan already noted, structure must be described in the "Parameters" section. If there is nothing implemented yet, then remove the parameter.
• Add an "Examples" section (see DOC: Add "Examples" to docstrings #7168).
• Add a "Notes" section with the appropriate versionadded markup.

For comparison, see the expm docstring in this file. A "Notes" section with just the versionadded markup is fine.

[ljwolf]

Thanks! Have done. "Structure" was a holdover from using the MatrixPowerOperator implementation, and it has been removed. I also refer now to the idea of your stackoverflow comment in the notes... I think it's pretty important to disclose this to the user.

[WarrenWeckesser]

The recursive/logarithmic implementation is probably a good default, and I don't suggest we change it, but folks involved in this PR might be interested to know that it is not always the best approach. Its performance can be substantially worse than simple iteration, depending on how the number of nonzero values increases as powers are computed. See my answer to a question on stackoverflow for an example where the recursive calculation of A**16 is much slower than computing A*A*A*A*A*A*A*A*A*A*A*A*A*A*A*A.

[ljwolf]

thanks for the comment @WarrenWeckesser! I think I've addressed the concerns mentioned. The point on your stackoverflow post is useful to disclose (imho) so I've added the gist of the point to the notes.

[ljwolf]

I should say, I'm happy to move forward with an inline-based computation, but my goal was only to bring forward the existing implementation. This recursive strategy was adopted in #513, and we probably should develop a more stringent benchmark to characterise the relative performance of a recursive vs. inline strategy over different sizes and sparsity levels. For now, we should probably proceed with the current _sp_matrix.__pow__ implementation, and optimise its implementation once it's brought forward to the sparse array.

[jjerphan]

LGTM.

I am also in favor of studying both possible implementations in a dedicated PR.

I just have one comment regarding the docstring.

[perimosocordiae]

As discussed, there are a few ways that this implementation could be made more efficient, but as it stands this is already a very nice improvement to the library: without a standalone matrix_power function, users of sparse arrays would be out of luck if they were migrating old sparse matrix code using the ** operator.

Thanks for the PR, @ljwolf !