Efficient way of creating matrices from a very large array while avoiding singular values

Question

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.

Answer 1

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}},
   
   2. Sum[1./Abs[N[i]] Abs[n - Abs[i]], {i, 1., N[n - 1]}],
   
   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