How to get inverse matrix for large size matrix in Mathematica?

Solutions to Hundred-dollar, Hundred-digit Challenge problems can be found easily:

I'm trying to solve the 7th problem:

Let A be the 20000×20000 matrix whose entries are zero everywhere except for the primes 2, 3, 5, 7, ..., 224737 along the main diagonal and the number 1 in all the positions $${\displaystyle a_{ij}}$$ with $${\displaystyle |i-j|=1,2,4,8,\dots ,16384}$$. What is the (1, 1) entry of $${\displaystyle A^{-1}}$$ ?

Below is my attempt:

Clear["Global*"];
matrixSize = 20;
(*Generate the prime numbers along the diagonal*)
primes = Prime[Range[1, matrixSize]];

(*Create a sparse diagonal matrix with the primes on the diagonal*)
A = SparseArray[Band[{1, 1}] -> primes, {matrixSize, matrixSize}];
A

(*Create the off-diagonal elements with 1s*)
kSet = Table[2^m, {m, 0, Floor[Log[2, matrixSize - 1]]}];
offDiagonals =
Flatten[Table[{Band[{1 + k, 1}] -> 1, Band[{1, 1 + k}] -> 1}, {k,
kSet}]];

(*Add the off-diagonal elements to the matrix*)
A = A + SparseArray[offDiagonals, {matrixSize, matrixSize}];
(A = Normal@A) // MatrixForm
Dimensions@A
(*Compute the inverse of the matrix A*)
Ainv = Inverse[A, Method -> "DivisionFreeRowReduction"];

(*Extract the (1,1) entry of the inverse matrix*)
Ainv[[1, 1]]


How to get InverseMatrix for large size Matrix (e.g., 20000) in Mathematica?

Any more effective method? Like Singular Value Decomposition (SVD) (LAPACK: GESVD)?

• Strassen's Algorithm (and related fast matrix multiplication methods)
• Blocked LU Decomposition with Partial Pivoting
• Cholesky Decomposition (for positive definite matrices)
• Iterative Methods (for large, sparse matrices)
• Approximate Methods (for very large matrices)
• Parallel and GPU-based Implementations
• Sparse Inverse Methods

As you can observe in the tiny instance of the problem, the matrix Ainv is dense. So you would need $$O(\mathrm{matrixSize}^2)$$ memory for storing the result. That becomes quickly infeasible. Moreover, computing all entries of the inverse matrix may cost even more than $$O(\mathrm{matrixSize}^3)$$, in particular when you use exact arithmetic.

A much faster way is to use a sparse direct solver with machine precision numbers:

S = LinearSolve[N[A]];


To extract the entry $$A^{-1}_{ij} = e_i^T A^{-1} e_j$$ you can do

UnitVector[matrixSize, i] . S[ UnitVector[matrixSize, j] ]


or

S[ UnitVector[matrixSize, j] ][[i]]


There is a good chance that this requires substantially less than $$O(\mathrm{matrixSize}^2)$$ time. (How much it really costs depends pretty much on the sparsity structure of the matrix.)

This might not be the most accurate method. Compared to the "solution" 0.7250783462 posted on Wikipedia, the relative error is 9.43361*10^-11. But the solutions is given only with 10 significant digits, so one cannot say how good my solution really is.

• With LinearSolve[N[A, 1000], UnitVector[matrixSize, 1]][[1]] (for example) you can get 1000 digits of accuracy in not much more time. Commented Aug 15 at 7:14
• Good point. Hm. I tried 100 digits of precision with n = 20000. I let it run in the background for some hours. Nothing came out of it but hot air. When I run it with 1000 digits for n = 800, then it finishes after a bit less than 2 minutes, taking more than 5000 times longer than in douple precision. It is truely amazing that it does work, but I find "in not much more time" a bit far fetched... ;) Commented Aug 15 at 12:23
• Your code does S = LinearSolve[N[A]], which is very slow when called with extra precision like S = LinearSolve[N[A, 1000]]. The code I gave avoids this step by calling LinearSolve with two arguments, and finishes very much faster. Commented Aug 15 at 18:48
• Ah, good point! Now I am curious what LinearSolve` is doing in this case. Hm. For example, it could do a factorization in machine precision and use that as preconditioner for GMRES... Commented Aug 15 at 20:18