# Large matrix inversion

I'm still learning Mathematica, and I'm trying to write a program that needs to invert a large (nxn, n=160 or possibly larger than this) matrix with numerical entries, multiply with to another matrix, diagonalize the resulting matrix and then compute the n ipr (inverse participation ratio) values from its n eigenvectors. Because at the end of the program I also need to find some statistical quantities (mean,median, min, max) on these ipr values, I am using a Do loop to repeat the aforementioned matrix operations 100 times, so that I obtain a 100-dim list (each element is itself a list of n ipr values), and average over each column, which will give me a list of n average ipr values.

Here's the program:

n = 160; t = 100;
H = Table[0.0, {2 n}, {2 n}];
ipr[v_] := Sum[Abs[v[[i]]]^4, {i, Length[v]}];
iprv = Table[0.0, {i, t}];

Do[
Print[s];
m = Table[RandomReal[{0, 0.5}], n];
k = RandomReal[{0, \[Pi]}];
M = Table[If[i == j, 2, 1/(Abs[i - j])^3], {i, 1, n}, {j, 1, n}];
F = Inverse[M];
Print["y"];
Psi = Abs[F.m];
h1 = Table[If[i == j, (1 - Cos[k]) + Psi[[i]], Sqrt[Psi[[i]] Psi[[j]]]/(2*(Abs[i - j])^3)],
{i, 1, n}, {j, 1, n}];
h2 = Table[If[i == j, Psi[[i]], Sqrt[Psi[[i]] Psi[[j]]]/(2*(Abs[i - j])^3)], {i, 1, n}, {j,
1, n}];
H = ArrayFlatten[{{h1, h2}, {-h2, -h1}}];
{h, u} = Eigensystem[H]; v = Pick[u, NonNegative[h]];
iprv [[s]] = Table[ipr[v[[i]]], {i, Length[v]}],
{s, t}
]

mipr = Mean[iprv];

max = Max[mipr]
min = Min[mipr]
me = Mean[mipr]
med = Median[mipr]
stat = {me, med, min, max};
Print[stat];


The problem is that this program runs very slow on my computer (it takes approx 5 min to do each loop), and I believe the problem is the large matrix inversion. Do you have any experience on inverting large matrices on Mathematica and any suggestion on how to speed it up? Is there anything else in the Do loop that can be optimized (even the Do loop itself), that might be slowing down the program?

Thank you so much for the help!

• You're taking the inverse of a matrix of rationals, and this is very slow. Since you take the dot product of the inverse with a vector of machine numbers, you should convert your matrix to machine numbers first. So, something like: M = Table[If[i == j, 2., 1./(Abs[i - j])^3], {i, 1, n}, {j, 1, n}]. Also, it seems you're taking an inverse of the same matrix during each iteration, that seems unnecessary. – Carl Woll Mar 18 at 23:32
• Please include the necessary definitions in your program (n=, etc.) – A.G. Mar 18 at 23:34
• Sorry I cut the first line of code by mistake. n = 160; t = 100; – ofey Mar 18 at 23:57
• Carl Woll thank you for the answer, you're right about matrix of rationals and repeating the matrix inversion, I'll try to change that – ofey Mar 19 at 0:01
• Actually you are solving a system of linear equations M.Psi=m using matrix inversion. This is very inefficient. For this purpose you should not use Inverse[M], but rather Psi=LinearSolve[M,m]. – yarchik Mar 19 at 8:20

To get the matrix inversion to work on machine numbers just replace

F = Inverse[0.0 + M];


That brings the running time to under 150 seconds on my MBPro -- replace the Do line by

Timing@Do[


to have Mathematica tell you about running time.

You can indeed shave a few seconds off by moving the matrix definition and inverse out of the Do loop, but not much.

Here are a couple of further suggestions to make the code faster and a bit shorter.

n = 160;
t = 100;
id = IdentityMatrix[n, SparseArray, WorkingPrecision -> MachinePrecision];
A = 0.5 ((DistanceMatrix[Range[1., n]] + id)^-3 - id);
M = 2. (A + id);
F = LinearSolve[M];
iprv = ConstantArray[0., {n, t}];

ParallelDo[
Block[{m, k, Psi, B, h1, h2, H, h, u, v},
m = RandomReal[{0, 0.5}, {n}];
k = RandomReal[{0, \[Pi]}];
Psi = Abs[F[m]];
B = Sqrt[KroneckerProduct[Psi, Psi]] A;
h1 = DiagonalMatrix[SparseArray[(1 - Cos[k]) + Psi]] + B;
h2 = DiagonalMatrix[SparseArray[Psi]] + B;
H = ArrayFlatten[{{h1, h2}, {-h2, -h1}}];
{h, u} = Eigensystem[H];
v = Pick[u, UnitStep[h], 1];
iprv[[s]] = Total[Abs[v]^4, {2}]
]
, {s, 1, t}, Method -> "CoarsestGrained"];


This takes about 1.5 seconds on my 4 Core machine.

That's what I did:

• What Carl Woll and A.G. pointed out already and which makes the biggest difference is converting M to a matrix of machine doubles.

• Instead of using Table, I used various ways to create the matrices. In contrast to Table, all these guarantee that the resulting matrices are stored as packed arrays.

• In general, it is a good idea to use LinearSolve instead of Inverse, both for performance and accuracy. LinearSolve computes an LU-factorization (stored in F within a LinearSolveFunction object), that can be reused multiple times. So factorizing M once in the beginning definitely saves some time.

• Also the matrix A is reused a couple of times. This avoids having to recompute the distance matrix in each iteration of the loop.

• Most of these methods (e.g., componentwise addition and maultiplication, KroneckerProduct, applying Sqrt or Power componentwise) are also vectorized.

• In particular, this applies to the replacement Total[Abs[v]^4, {2}] for ipr.

• If I am not mistaken, the computations in each loop iteration are fully independent. Thus, we can use ParallelDo to parallelize it. Note that there is a considerable overhead when a Parallel construct is called for the first time in a Mathematica session, because some packages are loaded on the fly.

I also observed that the eigenvalues seem to come always as pairs (on positive, one negative). This is due to the special structure of the matrix H. By exploiting this structure, it is likely that one can reduce your problem to solving a smaller $$n \times n$$ eigensystem. But I did not attempt to do it.

• These are all very useful comments, thank you! I will update my code with your suggestions. – ofey Mar 20 at 19:03
• You're welcome! – Henrik Schumacher Mar 20 at 19:11