Is something like `expm_multiply` available in PyLops?

[SimoneGasperini]

This is just to ask whether PyLops implements a function to compute something like $e^A \cdot v$ (i.e. the matrix exponential of a linear operator $A$ applied to the vector $v$). In scipy.sparse.linalg there is a function expm_multiply (see here) implementing the algorithm but taking a quick look to the documentation I was not able to find anything like that in PyLops.

Any kind of help would be appreciated!😃

[aman-coder03]

PyLops does not currently implement a function equivalent to scipy.sparse.linalg.expm_multiply.
PyLops operators are compatible with SciPy’s LinearOperator interface, though, so you can simply pass a PyLops operator directly to scipy.sparse.linalg.expm_multiply.

example:-
from scipy.sparse.linalg import expm_multiply
import pylops
A = pylops.SomeLinearOperator(...)
v = ...
y = expm_multiply(A, v)

as long as your PyLops operator implements matvec (and optionally rmatvec), expm_multiply will work out of the box.
so, the functionality isn’t provided inside PyLops itself, but you can use the SciPy routine seamlessly with PyLops operators.

[SimoneGasperini]

Thank you very much @aman-coder03 for your suggestion. However I think that trying to call the SciPy expm_multiply function passing a pylops.LinearOperator actually fails raising a TypeError. Take a look for example at the following:

import numpy as np
from pylops import LinearOperator
from scipy.sparse.linalg import expm_multiply

class Identity(LinearOperator):

    def __init__(self, shape):
        super().__init__(dtype=float, shape=shape)

    def _matvec(self, vec):
        return vec

dim = 2
linop = Identity(shape=(dim,dim))
vec = np.ones(dim)

expm_multiply(linop, vec, traceA=dim)

TypeError: unsupported operand type(s) for -: 'Identity' and 'float'

[mrava87]

@SimoneGasperini sorry for my late reply.

Whilst we strive to be 100% compatible with SciPy operators, it seems like something is not working here as expected... I suspect because adding a print before scipy/sparse/linalg/_expm_multiply.py:270

print(type(A), mu, type(ident))
A = A - mu * ident

gives

<class 'pylops.basicoperators.identity.Identity'> 0.9890340370149784 <class 'numpy.ndarray'>

and we can't sum an operator to a numpy array. Apparently if A is a scipy LinearOperator you get

<class 'scipy.sparse.linalg._interface._CustomLinearOperator'> 0.9424000000000001 <class 'scipy.sparse.linalg._interface.IdentityOperator'>

and so you can sum them. I think the culprit is in how _ident_like in scipy/sparse/linalg/_expm_multiply.py works... and they even say it is a hack (# A compatibility function which should eventually disappear.)... if they had added a method to the actual LinearOperator object we could overwrite it, but they didn't, so I am not sure we should try to fix it on our side...

As a workaround, although it is not identical to passing a PyLops operator directly, but with a simply wrapping of PyLops' matvec/rmatvec into a scipy LinearOperator you can make this method work.

Hope this helps 😄

import numpy as np
from pylops import Identity

from scipy.sparse.linalg import expm_multiply
from scipy.sparse.linalg import LinearOperator

dim = 2
linop = Identity(100)
linop = LinearOperator((100, 100), matvec=linop._matvec, rmatvec=linop._rmatvec)
vec = np.ones((100, dim))

expm_multiply(linop, vec)