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

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 = DeveloperToPackedArray[Outer[Subtract, ti, ti]]; // RepeatedTiming

{0.007300590625, Null}

m2 = DeveloperToPackedArray[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 = DeveloperToPackedArray[
Outer[Subtract, ti, ti]]; // RepeatedTiming; // RepeatedTiming

{0.74297815, Null}

m2 = DeveloperToPackedArray[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 = DeveloperToPackedArray[
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.

• 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]]; Commented Mar 5 at 18:20
• 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. Commented Mar 6 at 4:03
• Possibly the post that what Craig Carter refers to: mathematica.stackexchange.com/q/298086/38178 Commented Mar 6 at 4: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. Commented Mar 6 at 8: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. Commented Mar 6 at 9:05

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

• 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. Commented Mar 6 at 13:55
• D'oh. Well, that is why it is important to work use representative examples in posts like these... ;) Commented Mar 6 at 13:59