Addressing #84 and #85

[mp4096]

• Compute poles of a state-space system using numpy.linalg.eigvals as suggested in Calculate poles more smarter for state-space systems #84
• Replaced matrix inversion in the feedback computation by numpy.linalg.solve as suggested in Avoid using inv when not necessary #85
• Also replaced the singularity check with a more stable one. numpy.linalg.matrix_rank uses a reasonable numerical threshold, so that we don't need to check abs(det(F)) by hand. Determinant is actually very sensitive to matrix scaling, so better not use it. Both the determinant and the SVD required for matrix_rank are cubic, so there should be no slowdown.

[coveralls]

Coverage increased (+0.009%) to 74.042% when pulling 6aee194 on mp4096:master into 0b26e06 on python-control:master.

[mp4096]

Travis fails exactly at the same test as master, also see #70

[pmli]

Kudos for finding abs(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 trust numpy.linalg.solve to raise LinAlgError if singularity is detected? Do you have a use case when numpy.linalg.matrix_rank detects singularity, but numpy.linalg.solve doesn't?

[mp4096]

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

Consider this code:

import numpy as np
import matplotlib.pyplot as plt

matrix_size = 3

for spread in [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@A3   

    invA_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 matrix
    plt.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's DGESV, which in turn calls DGETRF to perform an LUP decomposition. DGETRF raises an exception only if a diagonal element is exactly 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 in numpy.linalg.solve

[pmli]

Thanks for the detailed explanation. I assumed NumPy/SciPy would emulate Matlab with its warnings about large condition numbers, using DGESVX. At least, there is an open NumPy 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.

[mp4096]

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 inverting D here).

[mp4096]

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