[Bug] sparse.linalg.LinearOperator.matmat cannot fallback properly to matvec

[homocomputeris]

Hi,
LinearOperator cannot properly use matvec to implement matmat.
Maybe somehow connected to issue #2434 because these calls seem to be a bit complicated:
matmat -> _matmat -> matvec -> _matvec -> matmat ...

Reproducing code example:

import numpy as np
import scipy as sp
import scipy.sparse.linalg as sla


dim = 8
shp = (dim, dim)

I = np.eye(*shp)

x0 = -sp.pi
xf = sp.pi
x = np.linspace(x0, xf, dim, endpoint=False)
# this is a vector representation of an operator with diagonal matrix
ref = (1-np.exp(-x))**2


def mv(y):
    # multiply this diag matrix by vector
    vv = ref
    return vv*y


def mm(z):
    # some crutches: explicit matmat function
   return np.array([mv(y) for y in z.T])

LO_with_mm = sla.LinearOperator(shp, matvec=mv, matmat=mm)
mm_test = LO_with_mm.matmat(I)
# multiplying by identity should result in the initial diagonal
print('works with explicit matmat:', mm_test - np.diagflat(ref))

# this won't work
LO_withot_mm = sla.LinearOperator(shp, matvec=mv)
mm_test = LO_withot_mm.matmat(I)

Error message:

traceback (most recent call last):

  File "<ipython-input-35-8eb4a3a9fdbf>", line 1, in <module>
    runfile('/home/user/Documents/Code/GniPy/scipy_bug.py', wdir='/home/user/Documents/Code/GniPy')

  File "/usr/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 705, in runfile
    execfile(filename, namespace)

  File "/usr/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 102, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/home/user/Documents/Code/GniPy/scipy_bug.py", line 43, in <module>
    mm_test = LO_withot_mm.matmat(I)

  File "/usr/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 326, in matmat
    Y = self._matmat(X)

  File "/usr/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 468, in _matmat
    return super(_CustomLinearOperator, self)._matmat(X)

  File "/usr/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 174, in _matmat
    return np.hstack([self.matvec(col.reshape(-1,1)) for col in X.T])

  File "/usr/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 174, in <listcomp>
    return np.hstack([self.matvec(col.reshape(-1,1)) for col in X.T])

  File "/usr/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 229, in matvec
    y = y.reshape(M,1)

ValueError: cannot reshape array of size 64 into shape (8,1)

Scipy/Numpy/Python version information:

1.0.0 1.14.0 sys.version_info(major=3, minor=6, micro=4, releaselevel='final', serial=0)

[lobpcg]

matmat falls back to matvec in

def _matmat(self, X):
# ...
    return np.hstack([self.matvec(col.reshape(-1, 1)) for col in X.T])

and that is where the given example breaks.

After a bit of investigation, it is not a bug, but a (probably undocumented?) feature. The trouble is that matvec(x) is designed to flexibly operate on a variety of different x: {matrix, ndarray} An array with shape (N,) or (N,1) (as well as a list? undocumented). Thus, matvec(x) is very forgiving and can be written in a sloppy way. However, when matvec needs to serve as a generator for matmat, suddenly the format of matvec must be strict in order to pass the line
np.hstack([self.matvec(col.reshape(-1, 1)) for col in X.T])
where X is a 2D array, because the argument x is now a 1D array of a specific shape, as far as I can see. E.g., in the example of this bug report, the matvec function call results in a 2D array output. This line can of course be edited to work on the provided example, by something like
return np.stack([self.matvec(col) for col in X.T])
but I do not see how to do it without braking the (undocumented?) backward compatibility.

Probably the best fix is to add a warning message that using matvec for matmat generation should be avoided for efficiency reasons if the line
np.hstack([self.matvec(col.reshape(-1, 1)) for col in X.T])
passes, and a detailed error message with a short example of a proper matmat function if it fails?

I guess, one can beat it to death by trying first the original
np.hstack([self.matvec(col.reshape(-1, 1)) for col in X.T])
if fails doing something like
return np.stack([self.matvec(col) for col in X.T])
if fails again trying something like
return np.stack([self.matvec(list(col)) for col in X.T])

[pv]

The LinearOperator documentation currently states:
"""
The user-defined matvec() function must properly handle the case
where v has shape (N,) as well as the (N,1) case. The shape of
the return type is handled internally by LinearOperator.
"""
The matvec defined in this issue handles the (N,1) case incorrectly, which causes the problem.

Also LO_withot_mm.dot(np.zeros([dim, 1])) fails, so what matmat does I think is not related.

[pv]

If something would be changed, matvec could in principle always reshape the input to (N,) shape (unless it's a np.matrix) to relax the shape requirement.

But this would be enhancement request and not a bug.

[lobpcg]

To preserve the strong backward compatibility, try first the original for (N,1):
np.hstack([self.matvec(col.reshape(-1, 1)) for col in X.T])
if fails turn into (1,N) via something like (untested)
return np.stack([self.matvec(col) for col in X.T])
if fails again reshape the input to like (N,) by something like (untested)
return np.stack([self.matvec(list(col)) for col in X.T])

Please advise, and I can make a quick PR to close this issue. (These are my Python exercises :-)

[pv]

The current behavior is correct AFAIK, as matvec is supposed to return (N,1) output for (N,1) input.
The default _matmat implementation appears correct to me.

In the example in this issue above, note that the matvec implementation there returns (8,8) size result for (8,1) input, so it's not a matrix-vector product.

[lobpcg]

yes. my proposal is for an enhancement, to allow sloppy matvec like this one

[pv]

As I see it, the enhancement would be to change y = self._matvec(x) to y = self._matvec(x.ravel() if not isinstance(x, np.matrix)) in def matvec. This would force ndarray inputs to (N,) shape, which the matvec implementations are supposed to be able to handle, which would remove the need to do also the other special case that caused problems here.

[lobpcg]

that might break strong backward compatibility

[lobpcg]

also would need to make rmatvec similar?

[lobpcg]

there is a similar bug/issue for rmatvec
ValueError: shapes (200,200) and (1,200) not aligned: 200 (dim 1) != 1 (dim 0)
#5546

[pv]

Right, it would break matvec implementations that handle (N,1) inputs but fail to handle (N,) ones. As you say, this would also be true for similar changes in matmat. Trying to work around it via try/except/etc fallbacks is generally not good design, so I think the choices are either status quo or potential breakage in corner case.

The iterative solvers etc. in scipy.sparse.linalg mostly feed in (N,) shape vectors, which makes that shape more probable to be OK in implementations.

---

Re: gh-5546, the bug in that case is I think in the implementation of matfuncs.py:MatrixPowerOperator._rmatvec which is incorrect when self._A is np.matrix --- it's incorrect both for (N,) and (N,1) inputs, I think all the code there just has not been written with np.matrix in mind so this might not be the only issue. (The right fix would be np.asarray on input)

[lobpcg]

I personally would not make a PR potentially breaking tons of codes that are not strictly complaint, but currently work. But you can of course feel free to try it.

If you do not like try/fail, my other proposal is:

Probably the best fix is to add a warning message that using matvec for matmat generation should be avoided for efficiency reasons if the line
np.hstack([self.matvec(col.reshape(-1, 1)) for col in X.T])
passes, and a detailed error message with a short example of a proper matmat function if it fails?

[pv]

I think there's not so much code that would break, but that is generally hard to find out since most code is not public. Efficiency warning here is probably quite noisy, and unclear if helpful for this issue (incorrectly implemented matvec). A safe improvement would be to make the error message clearer.

The general opinion around here is that np.matrix is for legacy code only (numpy nowadays emits a PendingDeprecationWarning if you try to create a np.matrix).

[rileyjmurray]

I read through this issue before creating a vaguely related GitHub issue of my own, and just wanted to make the following comment.

@pv regarding changing the matvec semantics from (n, 1) -> (n, 1) to (n,1) -> (n,):

I think there's not so much code that would break, but that is generally hard to find out since most code is not public.

The bigger issue is that the matvec codepath is accessible by python's @ operator. So changing matvec to follow (n,1) --> (n,) would violate the semantics of @ as introduced in PEP 465. It can confuse downstream developers when packages use non-compliant @ for matrix multiplication. If SciPy wasn't compliant with PEP 465, I think that would be a big issue.

I see that the current consensus for this issue is updating an error message, so I realize that maybe the discussion about (n, 1) --> (n, 1) vs (n,1) --> (n,) might no longer be on the table. In any case I feel strongly enough about SciPy conforming to PEP 465 that I wanted to express my concern.

[pv]

@rileyjmurray: note that the question here is about the internal _matvec API in LinearOperator, not the public __matmul__. The shape in the former has no effect on the latter.

[rileyjmurray]

Okay I see my mistake. I was worried because A @ x can induce the call sequence:

A.__matmul__ -> A.__mul__ -> A.__dot__ -> A.matvec -> A._matvec,

and so I thought the possible change to _matvec would be visible to the user. I see now that matvec performs the relevant reshaping ((n,1) vs (n,), depending on input), so the shape returned by _matvec isn't relevant when it comes to __matmul__.

Sorry for the bother!