Large matrix inversion

Question

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!

Answer 1

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.

Answer 2

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.