2
$\begingroup$

I need to compute the eigenvalues of matrices of size 1000x1000. I create two symmetric matrices $A$ and $B$ where each entry $(i,j)$ is given by a function ($f$ for $A$ and $g$ for $B$) of $i,j$ and a parameter $n$ over which I will later sum. The functions $f$ and $g$ are essentially Gauss Hypergeometric functions:

f[i_,j_,n_] := Hypergeometric2F1[0.5, 0.5, N[i-j+1], 1.-(L+1.)/n];
g[i_,j_,n_] := Hypergeometric2F1[-0.5, 1.5, N[i-j+1], 1.-(L+1.)/n]

Creating these two 1000x1000 matrices takes already a fair amount of time. Then I repeat the process for some different values of $n$, computing the eigenvalues of the matrix product $C=A.B$. At the end, I use these eigenvalues to compute some quantity $R$.

Here is the code I have written:

Do[
  L = 10*k;
  R = 0;
  ti = SessionTime[];

  Do[   
   (* Creating symmetric matrices A and B *)
   Atmp = ParallelTable[f[i,j,n], {i, 1, 1000}, {j, 1, i}];
   Btmp = ParallelTable[g[i,j,n], {i, 1, 1000}, {j, 1, i}];
   A = MapThread[Join, {Atmp, Rest /@ Flatten[Atmp, {{2}, {1}}]}];
   B = MapThread[Join, {Btmp, Rest /@ Flatten[Btmp, {{2}, {1}}]}];

   C = A.B;

   (* Eigenvalues of C *)
   Ev = Re[Eigenvalues[C]];

   (* Some computation with Ev *)
   Rint = 0;
   Do[
    If[Ev[[i]] > 1/2, Rint += Log[Ev[[i]]-1/2]];
    , {i, 1, 1000}];

   R += Rint;

   If[n == L, Print[L , "\t", R, "\t", SessionTime[] - ti, "s"]];

   , {n, 1, L}];

  , {k, 1, 20}];

For $L=10$ (ie $k=1$ in the first Do), it takes around 9500s, for $L=40$ it takes 61 000s, and so on.. For larger values of $L$, it will take several days! I cannot seem to find how to speed things up (ParallelDo does not improve the speed). An idea may be to compute the eigenvalues of my $L$ matrices at the same time in some sort of vectorized way, but I am not sure if it is possible or if it even makes sense. Do you have any idea?

$\endgroup$
3
  • 2
    $\begingroup$ Would you please also post the defintions of f and g? $\endgroup$ Apr 29, 2018 at 13:53
  • $\begingroup$ A quick thought. Are you using exact arithmetic when you don't need it? $\endgroup$
    – mikado
    Apr 29, 2018 at 13:53
  • $\begingroup$ No, I use MachinePrecision numbers everywhere. I will post simplified $f$ and $g$. $\endgroup$
    – Kaio
    Apr 29, 2018 at 13:59

1 Answer 1

1
$\begingroup$

Because of floating point underflow, it might not be a good idea to use machine precision numbers. This is why I define f and g like this.

f[i_, j_, n_, L_] := 
  Hypergeometric2F1[1/2, 1/2, i - j + 1, 1 - (L + 1)/n];
g[i_, j_, n_, L_] := 
  Hypergeometric2F1[-1/2, 3/2, i - j + 1, 1 - (L + 1)/n];

The matrices A and B are symmetric Toeplitz matrices and hence only their first row has to be computed. This cuts down the number of evaluations of f and g considerably.

We can obtain the matrices like this

prec = 20;
size = 1000;
n = 1;
L = 10;

a = Developer`ToPackedArray[N@f[Range[size], N[1, prec], n, L]];
b = Developer`ToPackedArray[N@g[Range[size], N[1, prec], n, L]];
A = ToeplitzMatrix[a, a];
B = ToeplitzMatrix[b, b];

Afterwards, you compute Tr[A.B] in a very expensive way by computing A.B and its eigenvalues. However, the trace is nothing but the Frobenius inner product of A with B and it can be computed from a and b without even assembling A and B:

mult = 2. Range[N[size], 1., -1];
mult[[1]] = N[size];

(a mult).b == Tr[A.B]

True

This is how the whole code looks like after refactorization (I also dropped the printing of intermediate results).

f[i_, j_, n_, L_] := Hypergeometric2F1[1/2, 1/2, i - j + 1, 1 - (L + 1)/n];
g[i_, j_, n_, L_] := Hypergeometric2F1[-1/2, 3/2, i - j + 1, 1 - (L + 1)/n];

size = 100;
mult = 2. Range[N[size], 1., -1];
mult[[1]] = N[size];

prec = 20;

data = ParallelTable[
Table[
 Block[{a, b, time, val, tic, toc},
  tic = SessionTime[];
  a = Developer`ToPackedArray[N@f[Range[size], N[1, prec], n, L]];
  b = Developer`ToPackedArray[N@g[Range[size], N[1, prec], n, L]];
  val = (a mult).b;
  toc = SessionTime[];
  {L, n} -> {"Value" -> val, "Time" -> toc - tic, "Kernel" -> $KernelID}
  ]
 , {n, 1, L}]
, {L, 10, 100, 10}]; // AbsoluteTiming

For size = 100, this runs through in about 12.5 seconds on my Quad Core CPU.

Note that this algorithm has linear complexity in size while yours had cubic complexity due to the use of Eigenvalues and the matrix-matrix product (if you don't count these in then it was quadratic due to the way the matrices A and B were generated).

The slowest part here remains the evaluation Hypergeometric2F1. Maybe one can utilize some functional equations to reduce the number of evaluation needed. I don't know.

$\endgroup$
6
  • $\begingroup$ I should have said that I actually require the eigenvalues of $A.B$ because I compute some quantity with it, not only the trace as it seems. I will modify my post to clarify this. Henrik, is there a way to get the eigenvalues kand keeping a linear complexity as you did? $\endgroup$
    – Kaio
    Apr 29, 2018 at 17:23
  • $\begingroup$ Maybe there is. I don't know. You can do the weirdest things with Toeplitz matrices and I am definiately not an expert in those. $\endgroup$ Apr 29, 2018 at 17:28
  • $\begingroup$ I will look into it. Can you also explain your first comment "Because of floating point underflow, it might not be a good idea to use machine precision numbers."? $\endgroup$
    – Kaio
    Apr 29, 2018 at 17:29
  • $\begingroup$ Ah, that was because at least g threw underflow errors when I was trying to call it with size = 1000. That happens when a floating point operation returns a value whose modulus is below $MinMachineNumber (that's approximately 2.22507*10^-308). $\endgroup$ Apr 29, 2018 at 17:34
  • $\begingroup$ Henrik, this line: {L, n} -> <|"Value" -> val, "Time" -> toc - tic, "Kernel" -> $KernelID|> produces an error with MMA 9.0 ( Syntax::sntxb: Expression cannot begin with "<|Value->val".) $\endgroup$
    – Kaio
    Apr 29, 2018 at 18:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.