1

There's a problem I'm working on for school where I need to solve a very large system of equations (1million x 1million matrix) where the solution is a 1 million vector of all ones.

A = [  1    1/2    1/3    1/4   ...   1/n
     1/2    1/3    1/4    1/5   ...   1/n+1
     1/3    1/4    1/5    1/6   ...   1/n+2
      .      .      .      .     .     .
      .      .      .      .     .     .
     1/n    1/n+1  1/n+2  1/n+3 ...   1/2n-1 ]

b = [1,1,1,...,1]^T  (lenght = n)


n = 1x10^6


Ax = b  <-- solve this system of equations

I've been able to solve this using Python and the Numpy module with n = 1x10^5 in about 16 minutes. However, when I try to solve for n = 1x10^6 it takes a very (very!) long time, 18 hours and counting, ram at 100%.

Is there an algorithm to solve this system of equations more efficiently? I believe numpy.linalg..solve() uses Gaussian elimination.

Thanks in advance

ortunoa
  • 111
  • 1
    Search for the Cholesky factorization. – nejimban Jun 29 '22 at 20:36
  • 1
    It's a desperate attempt ! Hilbert matrix is known to be very "ill-conditionned" even for $n \times n$ with $n=10$. You say it has worked for $n=10^5$... But have you attempted to plug the results you find into your equations ? Is it satisfying ? I have some doubts... – Jean Marie Jun 29 '22 at 20:56
  • Yeah I plugged them back and got an error of 10^-7 when n was 1000 and 10^-6 when n was 10k. – ortunoa Jun 29 '22 at 21:01
  • Do you definitely want to do this numerically? The matrix is invertible (the reciprocal determinants are OEIS A005249) and you can probably find an explicitish expression for the inverse written down somewhere – M T Jun 30 '22 at 12:35
  • ... In fact this question seems to give the answer right away – M T Jun 30 '22 at 12:59
  • This might be a basic question, but if I have the inverse Hilbert matrix nxn, and the results vector, how can I compute my variables? Ax = b, x = A/b ? – ortunoa Jun 30 '22 at 14:02
  • Once you get to that point where you are having trouble, you have crossed a threshold where you must use double precision. Could this be a precision issue? – Disintegrating By Parts Aug 05 '22 at 01:39

0 Answers0