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[
    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

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

Creating random orientations:

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


The corresponding eigenvector is


{-0.134387, -0.74292, 0.655751}

  • $\begingroup$ 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). $\endgroup$ – user52181 Apr 25 '19 at 9:35
  • $\begingroup$ Q should always be diagonalizable because it is symmetric. $\endgroup$ – Henrik Schumacher Apr 25 '19 at 9:35
  • $\begingroup$ 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? $\endgroup$ – user52181 Apr 25 '19 at 9:41
  • $\begingroup$ Thanks again, very instructive as always! $\endgroup$ – user52181 Apr 25 '19 at 10:05
  • $\begingroup$ You're welcome! $\endgroup$ – Henrik Schumacher Apr 25 '19 at 10:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy