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