1

I want to inverse a 401x401 (upper triangular square) matrix in python but run into issues when the matrix size exceeds 153x153. For all sizes below this, everything works like a charm but then the product of the inverted matrix by the initial matrix is not the identity matrix anymore.

import numpy as np
import math
from scipy import special

N=1
check = np.identity(N)

def calculate_matrix(N):
    matrix = np.zeros( (N,N) ) 
    matrix[0]=[1/(2**(j))*special.binom(j,j/2) if j%2==0 else 0 for j in range(N)]
    for p in range(0,math.floor((N-2)/2)+1):
        i=2*p+1
        matrix[i]=[(-1)**(2*q-p)/(4**q)*special.binom(2*q+1,(q-p)) if (q-i>=0 and q%2==1) else 0 for q in range(N)]
    for p in range(1,math.floor((N-1)/2)+1):
        i=2*p
        matrix[i]=[(-1)**(2*q-p)/(2**(2*q-1))*special.binom(2*q,(q-p)) if (q-i>=0 and q%2==0) else 0 for q in range(N)]
    return matrix

while (np.allclose(check, np.identity(N))==True and N<=402):
    N=N+1
    matrix = calculate_matrix(N)
    inverse_matrix = np.linalg.inv(matrix)
    check=np.matmul(inverse_matrix, matrix)

print("Maximum size of the matrix that can be inverted:", N-1)

Any idea what goes wrong?

Addition: it appears (thx user@2357112) that the matrix is ill-defined, hence the issue. A workaround if it exists would be much appreciated as it appears that partial pivoting is not a solution here.

Jep
  • 59
  • 7
  • Is it a computer limit or program limit? – Solar Mike May 09 '23 at 08:34
  • The 401 size is imposed by the data. How would you suggest I test this ? – Jep May 09 '23 at 09:18
  • 1
    Your matrices are very [ill-conditioned](https://en.wikipedia.org/wiki/Condition_number#Matrices) - inverting them is intrinsically hard and highly affected by any source of error. Try calling `np.linalg.cond(matrix)` to see the condition number, or printing `np.max(np.absolute(inverse_matrix))` to see how big the elements of `inverse_matrix` get. – user2357112 May 09 '23 at 10:27
  • thanks @user2357112. The maximum condition number from `np.linalg.cond(matrix)` before the inversion is not valid (within the numerical precision of `check`) is `1.617507632158823e+23`. It goes up after that. And the largest element of `inverse_matrix` is `2.0076553962575547e+22`. – Jep May 09 '23 at 10:38
  • Why do you need the explicit inverse of that large matrix ? – Yves Daoust May 09 '23 at 12:37
  • @YvesDaoust: I am interested in the individual contributions from (non-linear) higher-order responses to a stimulus through Fourier Transform (carried out on N data points). I need to disentangle the various contributions (not necessarily decreasing in amplitude with higher orders) to the Fourier weights. This would also enable me to characterize the precision of the current approximate frameworks (and assess a precision over speed ratio) – Jep May 09 '23 at 14:26
  • I see, thanks. So you are essentially interested in the largest elements. – Yves Daoust May 09 '23 at 15:47
  • As a worst case scenario yes. Ideally, an exact expression would enable me 1/ to carry out an exact evaluation of both the quantities of interest to me and 2/ an evaluation of the error of the current framework. – Jep May 09 '23 at 15:53

0 Answers0