Does scipy.linalg.inv check whether the matrix is diagonal?

user_4533188 0
Viewed 576

Getting the inverse of a diagonal matrix is very simple and does not require complex methods. Does scipy.linalg.inv check whether the matrix is diagonal before it applies more complex methods or do I need to check this myself?

scipy matrix python
2 Answers

As you can see the Github code of scipy.linalg.inv, function inv first calls 

getrf, getri, getri_lwork = get_lapack_funcs(('getrf', 'getri','getri_lwork'),

Then function getrf does it job to give the LU decomposition and so on. Now we have to investigate how getrf function gives the LU decomposition. Because if it checks if it's diagonal before to process the input matrix, then no need to check it yourself.

Function getrf is obtained by calling _get_funcs but I can't go further from there (_get_funcs is called with the following arguments _get_funcs(names, arrays, dtype, "LAPACK", _flapack, _clapack, "flapack", "clapack", _lapack_alias)).

I suggest that you run an experiment with a large diagonal matrix to compare the time given to spit the output with linalg and an inversion by hand.

Update (by question author):

import numpy as np
from scipy.linalg import inv
a = np.diag(np.random.random(19999))
b = a.copy()
np.fill_diagonal(a, 1/a.diagonal())
c = inv(b)

does not even require a time measuring tool: It it very obvious that inv is much slower... (that is surprisingly disappointing).

user_104847750

Please check: scipy.linalg.inv If you put scipy.linalg.inv in try except if it raises LinAlgError when matrix a is singular. The determinant for singular matrix it zero.

try:
    # your code that will (maybe) throw  scipy.linalg.inv(your matrix)

except np.linalg.LinAlgError as err:
    # It shows your matrix is singular
    # Its determinant of a matrix is equal to zero
    # The matrix does not have an inverse.
    # You can conclude if the matrix is diagonal or not

If the determinant of a matrix is equal to zero:

The matrix is less than full rank. The matrix is singular. The matrix does not have an inverse.

As manually:

def is_diagonal(matrix):
    #create a dummy matrix
    dummy_matrix = np.ones(matrix.shape, dtype=np.uint8)
    # Fill the diagonal of dummy matrix with 0.
    np.fill_diagonal(dummy_matrix, 0)

    return np.count_nonzero(np.multiply(dummy_matrix, matrix)) == 0

diagonal_matrix = np.array([[3, 0, 0],
                            [0, 7, 0],
                            [0, 0, 4]])
print is_diagonal(diagonal_matrix)
>>> True

random_matrix = np.array([[3, 8, 0],
                          [1, 7, 8],
                          [5, 0, 4]])
print is_diagonal(random_matrix)
>>> False

scipy.sparse.dia_matrix.diagonal returns the k-th diagonal of the matrix.

from scipy.sparse import csr_matrix
A = csr_matrix([[1, 2, 0], [0, 0, 3], [4, 0, 5]])
A.diagonal()
array([1, 0, 5])
A.diagonal(k=1)
array([2, 3])

Also, from scipy.linalg import block_diag creates diagonal matrix if input arrays are square, therefore if they are not square, it can not create diagonal matrix.

Please consider in Jupyter you can find out time complexity. %timeit yourfunctionname

user_22927670
Related