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

Question

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

• Mathematica (MathWorld Notebook)
• Python

---

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

Answer 1

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.