How do I find the rotation axis that maximises the angular momentum for a set of discrete points and velocities? I have a set of independent particles (we can assume with equal mass), distributed pseudo-randomly in 3D space, each with its own individual velocity.
What is the process by which I could determine the orientation of an axis that would maximise the angular momentum about that axis?
I want a method that it isn't just trial and error, looping over a grid of possible positions and rotation angles (I know how to do this, and it takes a long time!).
NB: the rotation axis should pass through the centre of mass.
 A: 
Angular momentum is an extensive quantity; i.e. the total angular momentum of any composite system is the sum of the angular momenta of its constituent parts

from here
So, if the COM of the set of particles were found and its velocity, we could find a position relative to that to maximise the angular momentum.
Since the magnitude of the angular momentum is $rmv$ where $m$ is the total mass and $v$ the velocity of the COM, it would be a maximum for a point furthest away from the COM on any line through the COM and at right angles to the velocity of the COM.
At this point, the orientation of the axis would be at right angles to both the velocity vector and the line joining the point to the COM.
A: From this derivation the total angular momentum of any number of particles is
$\vec{L} = \vec{R} \times M \vec{V} + \sum_{i} \vec{r_{i}}\times m_{i}\vec{v_{i}}$
where
$\vec{R}$ is the vector from the origin (or position of rotation axis) to the center of mass
$M\vec{V}$ is the center of mass momentum (calculate once)
$\vec{r_{i}}$ is the position of particle $i$ relative to the center of mass.
$\vec{v_{i}}$ is the velocity of particle $i$ relative to the center of mass (just the raw velocities are fine, if you don't want to rotate your coordinate system).
(The above is essentially what John Hunter's answer describes in words)
You can calculate the center of mass momentum ($M\vec{V}$) and spin angular momentum $\vec{S} = \sum_{i} \vec{r_{i}}\times m_{i}\vec{v_{i}}$ once and cheaply recalculate $\vec{L}$ for any desired $\vec{R}$.
The best axis will be parallel with the spin angular momentum so that it contributes maximally. This provides the optimal orientation given any distance. Then, choose a position so that $\vec{R}\times\vec{V}$ is also parallel to $\vec{S}$. You can use $\vec{R} \propto \vec{V}\times\vec{S}$.
This is where the question as currently asked is not uniquely specified, since you can arbitrarily extend the distance to arbitrarily increase $L = |\vec{L}|$. However, it's possible that you want to choose $\vec{R}=0$, by placing the axis at the position of the center of mass of the particles.
To calculate center of mass position, you use $\frac{\sum m_{i}\vec{x}_{i}}{\sum m_{i}}$ with data $\vec{x}_{i}$. This is needed to calculate $r_{i}$.
A: There may be a way, it's still a computer program way, but may be quicker.
The method would be to let the direction of the axis vary in small steps 0-1 for $a$, $b$, $c$ in a nested loop and compute the total angular momentum for an axis.  To do that we would need to find the vector of the red line for each particle $Q_i$
If the COM is (0,0,0) and the equation of the axis is $A = \left[\begin{array}{ccc}0\\0\\0\end{array}\right]+\left[\begin{array}{ccc}a\\b\\c\end{array}\right]t$
where $t$ is a parameter, then the value of $t$ for the red line to be perpendicular to the axis is found as follows:
$$\left[\begin{array}{ccc}a\\b\\c\end{array}\right].\left[\begin{array}{ccc}x_i-at\\y_i-bt\\z_i-ct\end{array}\right]=0$$
leading to the value of $t_i$ for each particle
$$t_i = \frac{ax_i+by_i+cz_i}{a^2+b^2+c^2}$$
that would find the right $t_i$ then the red line has vector $$r_i = \left[\begin{array}{ccc}x_i-at_i\\y_i-bt_i\\z_i-ct_i\end{array}\right]$$
and the angular momentum for particle $Q_i$ is then $mv_i \times r_i$ (cross product), then sum for all particles using vector addition and find the magnitude of the result.

Store the value and repeat in the nested loop, and keep values exceeding the best previous.
Let's say the best $\left[\begin{array}{ccc}a\\b\\c\end{array}\right]$ was $\left[\begin{array}{ccc}0.8\\0.4\\0.2\end{array}\right]$, then repeat with smaller step size between $\left[\begin{array}{ccc}0.7-0.9\\0.3 - 0.5\\0.1-0.3\end{array}\right]$
then reduce the step size, repeat etc...maybe it could be made automatic.  The program would then home in on the direction of the best axis.
A: If no rigid body rotation is involved here then there isn't such a thing. The system has one  angular momentum vector defined as
$$ \vec{L} = \sum_i \vec{r}_i \times (m_i \vec{v}_i) $$
where $\vec{r}_i$ is the position of each particle relative to the center of mass.
The angular momentum vector has its own direction and magnitude.
Now you ask,

What is the process by which I could determine the orientation of an axis that would maximise the angular momentum about that axis?

But there is no process to find angular momentum about some other axis, as this implies a rotation about some axis, and the premise of the question is without any specific rotation defined (each particle having their own independent velocity).
If you could specify rotational kinematics for the particles about arbitrary axis, then the process would as follows.

*

*Find the three angular momentum vectors resulting from three unit orthogonal rotation axis, $\vec{L}_x$, $\vec{L}_y$, $\vec{L}_z$.


*Combine the three vectors as the columns of a mass moment of inertia matrix
$$ \mathrm{I} = \left\{ \begin{matrix} \vec{L}_x & \vec{L}_y & \vec{L}_x \end{matrix} \right\} $$


*Find the eigen values and eigenvectors of the matrix, and the vectors form the principal directions of rotation. The direction associated with the largest eigenvalue is the direction that maximizes angular momentum.
In summary, you cannot apply rigid body properties to things that aren't rigid bodies. At least not in a simplistic way of just adding things up. Think of a chain floating in space, and what is the significance of the center of mass point? A force applied through the COM of a chain is both going to translate and rotate the chain, so the meaning of COM of a non-rigid body is different from what you think, and the same goes for concepts like mass moment of inertia tensors which are well defined for rigid bodies, but kind of meaningless for non-rigid bodies.
