Coverage increased (+0.009%) to 74.042% when pulling6aee194 on mp4096:master into0b26e06 on python-control:master.

Travis fails exactly at thesame test asmaster, also see#70

Kudos for findingabs(det(F)). Determinant gives no information about the distance to singularity, therefore it cannot be used to numerically check for singularity. SVD is definitely much better. But, wouldn't it be even better to drop the singularity check and trustnumpy.linalg.solve to raiseLinAlgError if singularity is detected? Do you have a use case whennumpy.linalg.matrix_rank detects singularity, butnumpy.linalg.solve doesn't?

Thanks for the hint. Indeed, an additional SVD seems to be superfluous. However, it can be really helpful.

Consider this code:

importnumpyasnpimportmatplotlib.pyplotaspltmatrix_size=3forspreadin [5,10,20]:A1=np.random.rand(matrix_size,matrix_size)A2=np.diag(np.logspace(-spread,+spread,matrix_size))A3=np.random.rand(matrix_size,matrix_size)A=A1@A2@A3invA_A=np.linalg.solve(A,A)message="Matrix size: {:d}, matrix rank: {:d}"print(message.format(matrix_size,np.linalg.matrix_rank(A)))# Should be an identity matrixplt.figure()plt.imshow(np.log10(np.absolute(invA_A)),interpolation="none")plt.colorbar()

So instead of raising an exception,numpy.linalg.solve just returns some numerical garbage. The problem is,numpy.linalg.solve relies on LAPACK'sDGESV, which in turn callsDGETRF to perform an LUP decomposition.DGETRF raises an exception only if a diagonal element isexactly zero:

INFO is INTEGER= 0:  successful exit< 0:  if INFO = -i, the i-th argument had an illegal value> 0:  if INFO = i, U(i,i) is exactly zero. The factorization    has been completed, but the factor U is exactly    singular, and division by zero will occur if it is used    to solve a system of equations.

Hence, numerical rank loss cannot be detected innumpy.linalg.solve

Thanks for the detailed explanation. I assumed NumPy/SciPy would emulate Matlab with its warnings about large condition numbers, usingDGESVX. At least, there is an openNumPy issue concerning exactly this, although it is quite old. Until it is closed, I don't see a way to cheaply estimate the condition number in NumPy which could replace SVD.

Oh, I didn't knew about this issue, thanks! But it doesn't seem they're going to do anything about this behaviour.

Anyway, most technical systems have very few (in terms of linear algebra 😉) inputs/outputs. I think this additional SVD-based rank check will makes things slower only if there are at least 1000 inputs/outputs (we're invertingD here).

Thank you! I think one can close#84 and#85 now.