Efficient computation of matrices involving large sums of KroneckerDelta's

Question

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}];

---

Edited after comments:

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

Answer 1

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 NDSolve`FEM`MapThreadDot):

f = u \[Function] NDSolve`FEM`MapThreadDot[
 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}