logm fails for singular matrix

[pv]

scipy.linalg.logm fails to compute the logarithm of the following (non-diagonalizable) matrix:

https://dl.dropboxusercontent.com/u/5453551/scipy-logm-fail.mat

EDIT: it's also singular, not only non-diagonalizable

from scipy import io, linalg
d = io.loadmat('scipy-logm-fail.mat')
lnA, errest = linalg.logm(d['A'])

# not diagonalizable:
w, v = linalg.eig(d['A'])
s = linalg.svdvals(v)
print(s.min()/s.max())
# -> 3.0988115615e-17

Matlab can do it, so probably Higham has an algorithm that works here.
Fails for Scipy version 0.13.2 (also 0.12.0 fails, so it's not a regression).

The current version of logm can do simpler non-diag cases, e.g. [[1, 1], [0, 1]], however.

[argriffing]

I see that in the scipy development version this exception is raised in the error estimation. This is a part of the logm code that I have not updated using Higham's algorithms and which I kept only for backwards compatibility of the interface.

As a short-term solution maybe it would be worth adding a keyword argument flag that tells logm to not attempt to compute an error estimate? In the long term, the error estimate should be computed using a logm condition number. Such algorithms exist but have not yet been implemented. The implementation will need to use abstract linear operators.

I took a look at this matrix, and when I disable the part of the logm code that ignores exceptions and sets the matrix to nans for backwards compatibility, I see that I get the exception "cannot find logm of a singular matrix." Checking np.linalg.cond(A) I see that the condition number is 6.0491328449e+37 which is far from 1 meaning that the matrix is badly conditioned with respect to matrix inversion which means that it is kind of singular.

Edit: The short answer is that the problem is matrix singularity not non-diagonalizability. Rows A[13], A[14], A[41], and A[42] are all exactly zero in the sense that all of their entries have binary representation 0, so the matrix is singular so it has no meaningful matrix logarithm. I googled a link that might be relevant http://math.stackexchange.com/questions/549832/logarithm-of-singular-matrix.

Edit 2: Apparently a complex singular matrix like A=[[0, 0], [1+1j, 1+1j]] can have a non-principal matrix logarithm B such that expm(B) == A. This is not implemented in scipy.

[pv]

Thanks, I didn't think about the properties of the matrix in detail.

It might be nicer if the error estimation would be skipped when the matrix contains nans, however.
Raising an exception is not so bad, however, if it's unimplemented.

[argriffing]

Added a PR that returns something more useable than NaNs for logm of a singular matrix, while also giving a warning because logm of a singular matrix is not more meaningful than log of 0.

[pv]

PR merged.

Handling singular cases with warnings is probably the best option, since numerically singularity of a matrix is usually only fuzzily defined.

[argriffing]

It's kind of subtle. For example I think that logm([[0, 0], [1, 1]]) should return [[-inf, 0], [inf, 0]] (with a warning) if log(0) returns -inf (with a warning). But the logm implementation is not sophisticated enough to do this, so it just uses large numbers instead of inf.

[pv]

I'm wondering if the rank warning should also be emitted for matrices that are "almost" singular.

This resembles the discussion that was going on with numpy.linalg.inv, ie. it currently raises an exception for exactly singular cases and otherwise is silent. Whereas, it could be more useful for the user if a warning was raised when the matrix is already almost singular (because by then, nearly all precision is already lost).

The situation with logm singularity may be a bit different, though, since you might only care about the backward error (which is a not so common case with inversion).

[argriffing]

I agree that exact comparison zero is usually not the right thing for floating point, and I also agree that the situation for logm is a bit different than for inv because the connection between 'nearness to singularity' and 'largeness of condition number' are different between these two functions.

I've added this #3287 which raises a warning if a diagonal entry of a similar triangular matrix is smaller than the value it would have been replaced by had it been exactly zero. This is a little complicated but at least it is somewhat internally consistent.