# Efficient computation of matrices involving large sums of KroneckerDelta's

I was wondering if one could benefit from Mathematica's rich linear algebra methods for diagonalizing 2nd rank tensors. Namely, in the context of systems (fluids) comprised of capsule-like particles, in order to quantify the amount of orientational order of the particles, the Nematic order parameter $$p$$ is chosen, which is $$\approx 1$$ when the particles are aligned towards a common vector (but positions uncorrelated) and $$\approx 0$$ when they are orientated randomly. The definition of $$p$$ is usually given as the largest eigenvalue of the 2nd rank tensor $$Q$$ given by:

$$Q_{ab} = \frac{1}{2n}\sum_{i=1}^n (3u_a^i u_b^i - \delta_{ab}) \tag{1}$$ where $$u_a^i$$ and $$u_b^i$$ are the a-th and b-th components of the orientation vector (normalized) of particle $$i$$ and $$\delta$$ is here the Kronecker delta.

An example to generate such a tensor with random oriented vectors (normalized):

randomVec = #*Normalize@RandomReal[{-1, 1}, 3] & (*[taken from here][1]*)
n = 2; (*particles count*)
u = Table[randomVec@1, {n}]
matQ = Table[
1/(2.*n)*
Sum[3*Part[u, i, a]*Part[u, i, b] - KroneckerDelta[a, b], {i, 1,
n}], {a, 3}, {b, 3}];


As pointed out by Henrik Schumacher in the comments, $$Q$$ is expected to be always diagonalizable as it is in fact a 3-by-3 symmetric matrix, therefore the eigendecomposition is not computationally challenging. Instead, the computation of the matrix itself may pose a performance problem, and can therefore be optimized (my approach given in the question is rather the most basic of going about it.). So the question is, given we know our target, the largest eigenvalue and its eigenvector, how could we optimize the computation and diagonalization of Q?

[1]: https://mathematica.stackexchange.com/a/13040/52181

• Could you add an example Q in Mathematica form for people to work with? Apr 25, 2019 at 8:41
• @ChrisK Sure, I've added one way of generating such matrix Q.
– user52181
Apr 25, 2019 at 9:24
• Q should always be diagonalizable because it is symmetric. Apr 25, 2019 at 9:35
• The performance problem here is not to diagonalize the matrix Q but to compute in the first place. Apr 25, 2019 at 9:41
• @HenrikSchumacher Thanks, I've re-iterated the question to incorporate your points.
– user52181
Apr 25, 2019 at 9:48

Creating random orientations:

SeedRandom[123];
n = 1000000;
u = RandomPoint[Sphere[{0, 0, 0}], n];


A helper function to compute many KroneckerProducts in a threaded way (utilizing the undocumented function NDSolveFEMMapThreadDot):

f = u \[Function] NDSolveFEMMapThreadDot[
ArrayReshape[u, Append[Dimensions[u], 1]],
ArrayReshape[u, Insert[Dimensions[u], 1, 2]]
];


That's how I interpret your definitions of Q and p:

Q = 0.5 (ConstantArray[3./n, n].f[u] - N@IdentityMatrix[3]);
{λ, U} = Eigensystem[Q];
order = Ordering[λ];
λ = λ[[order]];
U = U[[order]];
p = λ[[-1]]


0.00141718

The corresponding eigenvector is

U[[-1]]


{-0.134387, -0.74292, 0.655751}

• Can we still perform the eigendecomposition in an approximate manner even when $Q$ is not diagonalizable? (for particles in 3D, $Q$ is always a 3-by-3 matrix).
– user52181
Apr 25, 2019 at 9:35
• Q should always be diagonalizable because it is symmetric. Apr 25, 2019 at 9:35
• I completely missed that, nice observation! I checked your definition of $Q$ matches exactly mine given in the post now. Could you kindly add how you obtain the eigenvector corresponding to that eigenvalue?
– user52181
Apr 25, 2019 at 9:41
• Thanks again, very instructive as always!
– user52181
Apr 25, 2019 at 10:05
• You're welcome! Apr 25, 2019 at 10:22