4

I need to create matrix from a large array (>>10000 elements) for further use. I have tried Outer and ParallelTable commands, among which Outer seems to be faster when the array size is not so large;

Edit: Array ti below is given in terms Range for simplicity but in general it is intended to represent an arbitrary time series.

<< Developer`
tf = 100;
d = 1;
ti = N@Range[-tf, tf, d];
Nt = Length[ti];

m1 = Developer`ToPackedArray[Outer[Subtract, ti, ti]]; // RepeatedTiming


{0.007300590625, Null}


m2 = Developer`ToPackedArray[ParallelTable[Table[ti[[i]] - ti[[j]], {j, 1, Nt}], 
{i, 1, Nt}]]; // RepeatedTiming


{0.03681925625, Null}

But for larger array size I have got

tf = 1000;
d = 1;
ti = N@Range[-tf, tf, d];
Nt = Length[ti];

m1 = Developer`ToPackedArray[
    Outer[Subtract, ti, ti]]; // RepeatedTiming; // RepeatedTiming


{0.74297815, Null}


m2 = Developer`ToPackedArray[ParallelTable[Table[ti[[i]] - ti[[j]], {j, 1, Nt}], {i, 1, 
      Nt}]]; // RepeatedTiming


{0.1460876, Null}

which is run on a 64-core machine (ver. 12.3). I was wondering if there is a more efficient way of doing this. I also need to construct the matrix:

tf = 1000;
d = 1;
ti = N@Range[-tf, tf, d];
Nt = Length[ti];
m3 = Developer`ToPackedArray[
    ParallelTable[Table[If[i != j, 1/(ti[[i]] - ti[[j]]), 0], {i, 1, Nt}], {j, 1, 
      Nt}]]; // RepeatedTiming

{0.8566911, Null}

where diagonal values are to be excluded. Apparently, the if statement inside Table introduces a slow down. I have seen a similar problem posted here, but it is not obvious to me how to implement those solutions inside ParallelTable.

user91411
  • 400
  • 3
  • 11
  • 5
    I saw this in another thread but cannot find it. Note: ToPackedArray[Outer[Subtract, ti, ti]] is about 750 slower than ToPackedArray[Outer[Plus, ti, -ti]]; – Craig Carter Mar 05 '24 at 18:20
  • 2
    I am 97.5 % sure that this is an XY question and that you actually don't really need to form this matrix explictly. m1 is a rank-1-matrix and there are many ways to exploit this. m2 is a Toeplitz matrix, and there also many algorithms that avoid forming it in the first place (e.g., some use fast Fourier transform). So I think it is quite likely that you don't have to brute-force it. – Henrik Schumacher Mar 06 '24 at 04:03
  • Possibly the post that what Craig Carter refers to: https://mathematica.stackexchange.com/q/298086/38178 – Henrik Schumacher Mar 06 '24 at 04:08
  • @CraigCarter Thanks for the tip. Indeed it is quite faster and less memory consuming than ParallelTable. But I don't understand why on earth a minus would make a dramatic difference, it is an odd behavior. – user91411 Mar 06 '24 at 08:58
  • @HenrikSchumacher Thanks for pointing that out. It is indeed Toeplitz matrix. My problem is to evaluate the integral of the form $\int d x d y I (x-y)$. After having discretized the whole integrand my strategy is to compute the Grand sum of the matrix $M_{ij} = dx_i d x_j I_{i j}$. So I need to construct the matrices, I don't see atm how FFT could be relevant for the problem. – user91411 Mar 06 '24 at 09:05
  • I see. So you might (or might not) have to compute all I_{ij}. But you don't have to store them all at once in memory. Instead, you create on accumulation variable sum = 0. and then run Do[sum += f[i,j], {i, -tf, tf, s}, {j, -tf, tf, s}]. Problem is that this is a loop and loops are slow in Mathematica. But they are fast in Compile. That spares you from working on large memory chunk - which have to be loaded chunkwise into the processor, first. And this kind of bulk memory is much, much slower than the little scratch space that the processor can work on directly. – Henrik Schumacher Mar 06 '24 at 13:12
  • Note that the runtime of Do[sum += f[i,j], {i, -tf, tf, s}, {j, -tf, tf, s}] is still quadratic in ft. The point here is to involve only fast components of the machine. – Henrik Schumacher Mar 06 '24 at 13:16
  • @HenrikSchumacher Thank you for the suggestion. For the grand sum I intend to use the function Total[M,2] as $M$ given above, and keep the $M$ as a packed array (which can be turned into triangular form before the evaluation due to symmetry) . Total seems to work okay but haven't gauged its performance for larger dimension of $M$, (say over $10^4$) . Currently, I am scratching my head on getting m3 in the post in an efficient way. – user91411 Mar 06 '24 at 13:24

1 Answers1

4

From the comments I read that you actually only need to sum over all entries of a matrix like

A = Table[f[i-j], {i, 1,n}, {j, 1, n}]

That is, you actually mean to compute

Sum[f[i-j], {i, 1, n}, {j, 1, n}]

As I said already in the comments, A is a Toeplitz matrix, i.e., it is constant on each diagonal. So the sum over the i-th diagonal can be computed by evaluating f[i] and just multiplying with the length of the diagonal:

Sum[ f[i] Abs[n - Abs[i]], {i, -(n-1), (n-1) }]

This costs you $O(n)$ floating point operations instead of O(n^2). You can speed this up further by using Compile if f is compilable.

If f[i] is an even function, i.e. if f[-i] = f[i]``, then you can get rid of you can even reduce your overall work roughly by half:

f[0] n + 2 Sum[ f[i] Abs[n - Abs[i]], {i, 1, (n-1) }]

That can also help you to avoid the If check on the main diagonal. Of course, if f[0] = 0 you can use.

2 Sum[ f[i] Abs[n - Abs[i]], {i, 1, (n-1) }]

If f[-i] = f[i], then the net sum is 0, so you don't have to do anything. In general, every function f[i] can be decomposed into f[i] = g[i] + h[i], where g is even and h is odd. Sometimes you can exploit this by extracting g and just summing over g.

Finally, here is a runtime example:

f = x |-> If[x == 0, 0., 1./Abs[N[x]]];
n = 10000;

cSum = Compile[{{n, _Integer}},


    

  1. Sum[1./Abs[N[i]] Abs[n - Abs[i]], {i, 1., N[n - 1]}],
    2. 



CompilationTarget -> "C"
   ];


sum0 = Total[Table[f[i - j], {i, 1, n}, {j, 1, n}]]; // AbsoluteTiming // First


sum1 = Sum[f[i - j], {i, 1, n}, {j, 1, n}]; // AbsoluteTiming // First


sum2 = 2. Sum[f[i] Abs[n - Abs[i]], {i, 1, (n - 1)}]; //

RepeatedTiming // First

sum3 = cSum[n]; // RepeatedTiming // First

126.279

125.448

0.000680404

0.000126996

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309
  • Very interesting! and thanks for posting this, but I should apologize for the confusion; I should have pointed out that although specific part of the integral at hand can be written as Toeplitz my concern here is not specifically related to the Toeplitz matrix. The other integrands I need to construct in the matrix form do not necessarily enjoy such symmetry. That's why I intend to use Total[M,2] for general applicability. That being said, I will take note of this feature for future use. – user91411 Mar 06 '24 at 13:55
  • D'oh. Well, that is why it is important to work use representative examples in posts like these... ;) – Henrik Schumacher Mar 06 '24 at 13:59