ENH:linalg: expm overhaul and ndarray processing

[ilayn]

Closes #354
Closes #4897
Closes #12838

What does this implement/fix?

This is a complete overhaul of the function scipy.linalg.expm. Two important properties are the performance increase and being able to handle ndarrays with n>=2. Moreover, diagonal and triangular matrices are recognized and treated as such. The existing tests pass with very tiny adjustments.

If you can provide some benchmark feedback for the new and the old version, that'd be great. It doesn't have to be scientific, "yes it is fast" "no I don't see any difference" is also fine.

Produced via:

import perfplot
import numpy as np
from scipy.sparse.linalg import expm as EX
from scipy.linalg import expm


perfplot.show(
    setup=lambda n: -np.random.rand(n, n),
    kernels=[EX, expm],
    labels=["SciPy", "this PR"],
    n_range=[3*x for x in range(1, 30)],
    logx="auto",
    logy=True,
    equality_check=None,
    xlabel="input size (n)",
}

Additional information

There are some tests missing for the new functionality. I will try to add them asap.

Things left out

I ran out of steam at some point since the logic became too much for me to hold it in my head, but hermitian matrices through a Schur factorization can also benefit from some performance increases. At some point I hope I will come back to this.

[akhmerov]

Does this change not break backwards compatibility?

[ilayn]

Ideally people should be using scipy.sparse.linalg.expm and dense version kind of hijacks the sparse version.

[adeak]

Just some nitpicks. And I can't review the cython bits.

[adeak]

Do these changes signal

1. an API change (scipy.linalg.expm used to be sparse, now it's dense) or
2. a mistake in the original tests using dense expm?

[ilayn]

I should explain this somewhere indeed. Previously, we were using sparse implementation in case of a dense input. Here dense gets its own implementation such that sparse one can live in its own environment.

[adeak]

Yeah, that then sounds like something that can easily trip up sparse users. scipy.linalg.expm no longer being sparse, that is. And we'd prefer it for people to use the top-level version of a name where possible, so scipy.linalg.expm was probably the preferred version in the original.

[mdhaber]

Hi @ilayn! Here are the results on my machine:

perfplot seems handy.

[ilayn]

Whoa! That's really good and bad :) Good for small matrices bad for large ones because around 360 we switch to estimation. I'll check again. Thanks Matt good data point for me.

[adeak]

Well this is interesting...

Debian linux with openblas. That zigzag makes me think something weird's going on.

[ilayn]

Yes, OpenBLAS behaviors are not very easy to explain unfortunately. I've been hit with a different kind of an issue but still no idea what the problem might be. Once 1.8.0 is out I'll try to dissect more this interesting behavior. At least the small matrix behavior looks consistent across MKL and OpenBLAS so far which is promising.

[Kai-Striega]

here you go, I hope it's helpful

[ilayn]

Thanks @Kai-Striega . It looks like you have MKL right?

[Kai-Striega]

@ilayn I'm not sure how to check, looking at the output of show_config it looks like I'm not using MKL

>>> from scipy import show_config
>>> show_config()
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/home/kai/anaconda3/envs/scipy-dev/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/home/kai/anaconda3/envs/scipy-dev/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/home/kai/anaconda3/envs/scipy-dev/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/home/kai/anaconda3/envs/scipy-dev/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
    not found = AVX512F,AVX512CD,AVX512_KNL,AVX512_KNM,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL

[ilayn]

Ok apart from one consistently timing-out test this is all good to go from my side.

[tylerjereddy]

Could we get an approval from one of the reviewers? PR is pretty big, not trying to block it

[ilayn]

No worries @tylerjereddy if you feel like this is getting rushed, I'm OK with the milestone bump. Come to think of it maybe it is better to finish the full suite of matfuncs, solve and eig together in one release.

[rgommers]

No worries @tylerjereddy if you feel like this is getting rushed, I'm OK with the milestone bump. Come to think of it maybe it is better to finish the full suite of matfuncs, solve and eig together in one release.

Let's bump this one to 1.9.0 then, thanks @ilayn

[ilayn]

I tried to immitate the meson.build from _decomp_update structure but it fails if I use

  linalg_init_cython_gen.process('_matfuncs_expm.pyx'),

on one step and fails on another step if I use

  linalg_init_cython_gen.process('_matfuncs_expm_pyx'),

Can someone see my mistake here ?

Nevermind found it.

[ilayn]

Benchmarks were the last thing that needed tweaking. Hence this is final for review.

Current benchmark for sparse (existing expm) and dense (this PR)

              ===== ======== ==========
                n    format            
              ----- -------- ----------
                30   sparse   19.6±0ms 
                30   dense    520±0μs  
               100   sparse   76.0±0ms 
               100   dense    3.55±0ms 
               300   sparse   810±0ms  
               300   dense    50.7±0ms 
              ===== ======== ==========

[ilayn]

Hmm don't know why the sudden pytest mark errors appeared. Doesn't seem to be related.

[ilayn]

It's been quite a while that I got all green in CI so I'd like to celebrate it by merging tomorrow if everybody is OK with it :)

[rgommers]

I just had a quick browse through. I think that makes sense. It's too much code to review in detail for any other maintainer I think - but it's a rewrite of an existing function, so it's reasonable to rely on test coverage and green CI.

[ilayn]

Indeed I wish I had some more eyes on it but given the state of the overhaul (which is basically calling lots of LAPACK functions) I can understand that only thing to do for me is to brace for impact on 1.9 release for things that the tests don't cover.

[ilayn]

OK fingers crossed. Thanks all !