Shouldn't scipy.linalg.sqrtm raise an error if matrix is singular?

[Acebulf]

I'm a bit confused by scipy.linalg.sqrtm not raising an error if matrix is singular, since it can't compute anything (returns NaN array). Other functions such as scipy.linalg.inv will raise scipy.linalg.LinAlgError if matrix is singular and operation can not be performed.

Why are errors ignored here?

*edit: I was under the assumption that it did not do ANY singular matrices, but it apparently does some of them. Even when it fails it won't raise an error, however.

[argriffing]

Just to be picky, it's not true that if the matrix is singular it can't compute anything. For example the matrix square root of [[0]] is [[0]]. But there are some matrices like [[0, 1], [0, 0]] for which the matrix square root doesn't exist.

[Acebulf]

It will return an array filled with NaN or infinite for some matrices which still have a square root. Which is why having an error returned would be nice when the function doesn't find a square root.

For example the matrix [[1,0,0,1],[0,0,0,0],[0,0,0,0],[1,0,0,1]] has a square root [[sqrt(0.5),0,0,sqrt(0.5)],[0,0,0,0],[0,0,0,0],[sqrt(0.5),0,0,sqrt(0.5)]] but scipy.linalg.sqrtm will yell about singularity and then return a matrix with every item being -9223372036854775808, yet without raising LinAlgError.

[pv]

Raising Linalgerrors is not useful behavior, as the concept "singular matrix" is ill-defined in floating point. Rather, the algorithm returns (set disp=False) the square root, and an error estimate. . Returning -9223372036854775808 instead of nans is due to supplying integer matrix as input. This is a bug. . That the algorithm fails to find square roots for certain matrices is also a bug.

[argriffing]

@Acebulf Thanks for reporting this, I can reproduce the problem.

>>> import numpy as np                                                          
>>> from scipy.linalg import sqrtm                                              
>>> a = np.array([[1, 0, 0, 1], [0, 0, 0, 0], [0, 0, 0, 0], [1, 0, 0, 1]], dtype=float)
>>> a
array([[ 1.,  0.,  0.,  1.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  1.]])
>>> b = np.sqrt(0.5) * a
>>> sqrtm(a)
Matrix is singular and may not have a square root.
array([[ nan,  nan,  nan,  nan],
       [ nan,  nan,  nan,  nan],
       [ nan,  nan,  nan,  nan],
       [ nan,  nan,  nan,  nan]])
>>> b.dot(b)
array([[ 1.,  0.,  0.,  1.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  1.]])
>>> from scipy.linalg import eigvalsh
>>> eigvalsh(a)
array([ 0.,  0.,  0.,  2.])

For what it's worth, your matrix is positive semi-definite so you could use https://github.com/scipy/scipy/pull/3556 to work around this problem.

[Acebulf]

There is already a failflag variable if no square root is found (independent of whether matrix is singular), but it does nothing.

In the meantime I submitted a patch which calls LinAlgError when a square root is not found. There is the ignoreerrors parameter which when True yields practically the same behavior as now (except for the failflag no longer being ignored if the matrix is singular.) This parameter defaults to False, but it might be wise to make it default to True in order to ensure scripts which already handled the NaNs still work.

[argriffing]

There is already a failflag variable if no square root is found (independent of whether matrix is singular), but it does nothing.

I just want to confirm that this flag does nothing. At first I thought it could possibly come into play if sqrtm failed in a particular way for a nonsingular input matrix, but now I see that this cannot happen (the square roots of two entries on the diagonal of the Schur form of the input matrix would have to sum to zero, but this can only happen if those entries are both zero, and this would mean the matrix has two eigenvalues that are zero so the matrix would be singular).

[argriffing]

Returning -9223372036854775808 instead of nans is due to supplying integer matrix as input. This is a bug.

I agree with this.

That the algorithm fails to find square roots for certain matrices is also a bug.

This is true, and it is a limitation of the current sqrtm algorithm not a mistake in its implementation. The reference http://eprints.ma.man.ac.uk/1926/01/covered/MIMS_ep2012_26.pdf assumes there are no eigenvalues on the closed negative real line. Maybe it should automatically fall back to a different implementation?