Rotations for geology

PLEASE DO NOT REDISTRIBUTE WITHOUT PERMISSION
2014 May 7 by Joshua R. Davis, jdavis@carleton.edu

Conversion among representations of rotations

You can compute the exponentials and logarithms for rotations of 3D space using special methods (see Fletcher et al. (2003), “Gaussian distributions on Lie groups and their application to statistical shape analysis”). The rotationLog function here is about 2 or 3 times as fast as diagonalization. The rotationExp function here is about 1 or 2 times as fast as Mathematica’s built-in MatrixExp function.
	rotationFromAngleAxis returns a rotation matrix for the specified angle (in radians) and axis (unit Cartesian vector), according to the right hand rule. rotationFromVectors takes two unit Cartesian vectors, and returns the smallest rotation matrix that takes the first to the second. rotationFromVectorPairs takes independent vectors u0, v0 and independent vectors u1, v1, and returns the rotation that takes the initial vectors to the final vectors (assuming that the angles match, so that the rotation exists).

Clear[rotationExp, rotationLog, rotationFromAngleAxis, angleAxisFromRotation, volumetricFromAngleAxis, angleAxisFromVolumetric, volumetricFromRotation, rotationFromVolumetric, quaternionFromAngleAxis, angleAxisFromQuaternion, rotationFromQuaternion, quaternionFromRotation, rotationFromVectors, rotationFromVectorPairs];

(* Converts between an antisymmetric W and its rotation matrix exp W. *)
rotationExp[w_] :=
	Module[{theta},
		theta = Sqrt[Tr[Transpose[w].w] / 2];
		If[theta == 0,
			IdentityMatrix[3],
			IdentityMatrix[3] + (Sin[theta] / theta) w + ((1 - Cos[theta]) / theta^2) w.w]];
rotationLog[rot_] :=
	Module[{a, u},
		{a, u} = angleAxisFromRotation[rot];
		a {{0, -u[[3]], u[[2]]}, {u[[3]], 0, -u[[1]]}, {-u[[2]], u[[1]], 0}}];

(* Converts between rotation matrices and angle-axis representations, which are *)
(* of form {angle, axis}, where angle is in radians and axis is unit 3D Cartesian. *)
rotationFromAngleAxis[{alpha_, axis_}] :=
	Module[{u},
		u = {{0, -axis[[3]], axis[[2]]},
			{axis[[3]], 0, -axis[[1]]},
			{-axis[[2]], axis[[1]], 0}};
		IdentityMatrix[3] + Sin[alpha] u + (1 - Cos[alpha]) u.u];
angleAxisFromRotation[r_] :=
	Module[{cos, u1, u2, u3},
		cos = 0.5 * (Tr[r] - 1.0);
		If[cos == 1.0,
			{0, {0, 0, -1}},
			u1 = Sqrt[Max[0, 1 + (r[[1, 1]] - 1) / (1 - cos)]];
			If[r[[3, 2]] < r[[2, 3]],
				u1 = -u1];
			u2 = Sqrt[Max[0, 1 + (r[[2, 2]] - 1) / (1 - cos)]];
			If[r[[1, 3]] < r[[3, 1]],
				u2 = -u2];
			u3 = Sqrt[Max[0, 1 + (r[[3, 3]] - 1) / (1 - cos)]];
			If[r[[2, 1]] < r[[1, 2]],
				u3 = -u3];
			{arcCos[cos], unit[{u1, u2, u3}]}]];

(* Converts between the angle-axis and volumetric representation. The latter is *)
(* a vector equaling rho times the axis of the rotation, where rho increases *)
(* from 0 to (3 / (4 Pi))^(1 / 3) as the angle increases from 0 to pi. *)
(* Volumetric coordinates give a volume-respecting depiction of SO(3) as the *)
(* ball of radius (3 / (4 Pi))^(1 / 3) in R^3. *)
volumetricFromAngleAxis[{theta_, u_}] :=
	((theta -Sin[theta]) * 3 / (4 Pi^2))^(1 / 3) * u;
angleAxisFromVolumetric[vol_] :=
	Module[{rho, a, alpha},
		rho = Norm[vol];
		If[rho == 0,
			{0, {1, 0, 0}},
			Off[Solve::ratnz];
			alpha = Solve[{a - Sin[a] == rho^3 (4 Pi^2 / 3), 0 ≤ a ≤ Pi + 10^-6}, a, Reals];
			On[Solve::ratnz];
			{alpha[[1, 1, 2]], unit[vol]}]];
volumetricFromRotation[rot_] :=
	volumetricFromAngleAxis[angleAxisFromRotation[rot]];
rotationFromVolumetric[vol_] :=
	rotationFromAngleAxis[angleAxisFromVolumetric[vol]];

(* Converts between rotation matrices and angle-axis pairs and *)
(* quaternions q = q[[1]] + q[[2]] i + q[[3]] j + q[[4]] k. *)
quaternionFromAngleAxis[{alpha_, axis_}] :=
	{Cos[alpha / 2], Sin[alpha / 2] axis[[1]], Sin[alpha / 2] axis[[2]], Sin[alpha / 2] axis[[3]]};
angleAxisFromQuaternion[q_] :=
	Module[{alpha, sine, v, axis},
		alpha = 2 arcCos[q[[1]]];
		sine = Sin[alpha / 2];
		v = q[[2;;4]];
		If[sine==0.0,
			If[Abs[v[[1]]] > Abs[v[[2]]] && Abs[v[[1]]] > Abs[v[[3]]],
				If[v[[1]] > 0,
					{alpha, {1, 0, 0}},
					{alpha, {-1, 0, 0}}],
				If[Abs[v[[2]]] > Abs[v[[1]]] && Abs[v[[2]]] > Abs[v[[3]]],
					If[v[[2]] > 0,
						{alpha, {0, 1, 0}},
						{alpha, {0, -1, 0}}],
					If[v[[3]] > 0,
						{alpha, {0, 0, 1}},
						{alpha, {0, 0, -1}}]]],
			{alpha, v / sine}]];
rotationFromQuaternion[q_] :=
	{{q[[1]]^2 + q[[2]]^2 - q[[3]]^2 - q[[4]]^2,
		2(q[[2]] q[[3]] - q[[1]] q[[4]]),
		2(q[[1]] q[[3]] + q[[2]] q[[4]])},
	{2(q[[1]] q[[4]] + q[[2]] q[[3]]),
		q[[1]]^2 + q[[3]]^2 - q[[2]]^2 - q[[4]]^2,
		2(q[[3]] q[[4]] - q[[1]] q[[2]])},
	{2(q[[2]] q[[4]] - q[[1]] q[[3]]),
		2(q[[1]] q[[2]] + q[[3]] q[[4]]),
		q[[1]]^2 + q[[4]]^2 - q[[2]]^2 - q[[3]]^2}};
quaternionFromRotation[rot_] :=
	quaternionFromAngleAxis[angleAxisFromRotation[rot]];

(* Returns the smallest rotation matrix that takes u0 to u1. *)
(* Inputs are assumed unit 3D Cartesian. *)
rotationFromVectors[u0_, u1_] :=
	Module[{axis, norm},
		axis = Cross[u0, u1];
		norm = Norm[axis];
		If[Chop[norm] != 0, 
			(* This is the typical case. *)
			rotationFromAngleAxis[{ArcCos[Dot[u0, u1]], axis / norm}],
			If[Dot[u0, u1] > 0,
				(* The vectors are parallel. *)
				IdentityMatrix[3],
				(* The vectors are antiparallel. *)
				rotationFromAngleAxis[{Pi, orthogonalVector[u0]}]]]];

(* Returns the rotation taking u0 to u1 and v0 to v1. *)
(* Inputs are 3D Cartesian, with equal angles between them. *)
rotationFromVectorPairs[u0_, u1_, v0_, v1_] :=
	Module[{rot0, rot1},
		rot0 = Orthogonalize[{u0, v0, Cross[u0, v0]}];
		rot1 = Orthogonalize[{u1, v1, Cross[u1, v1]}];
		Transpose[rot1].rot0];

Geodesic methods

rotationDistance computes the distance between two given rotations, in the usual metric on SO(3). In other words, it returns the angle of the smallest rotation that takes one rotation to the other.
	rotationDiameter computes the diameter of a set of rotations (in the usual metric on SO(3)). That is, it finds the two rotations that are farthest apart, and reports the distance between them. This function is not fast.
	rotationMean computes the mean of a set of rotations. The first argument is a nonempty list of rotation matrices. The second argument is an error tolerance epsilon; 10^-10 is reasonable. The third argument is a bound on how many iterations are allowed; 100 is reasonable, although I’ve never seen more than 14 used. The function returns a list of three outputs. The first output is the mean rotation matrix. The second is the error achieved. The third is the number of iterations used. The latter two outputs can be used to diagnose problems. Problems may arise if the diameter is larger than pi. The mean here is the Frechet mean (also called the geodesic L2 mean). The algorithm is from Manton (2004), reproduced in Hartley et al. (2013). Supposedly Krakowski et al. (2007) has a faster algorithm.
	rotationVariance computes the variance of a set of rotations. The first argument is the set of rotations. The second argument is the mean rotation, probably computed by rotationMean. For the definition of variance see Krakowski et al. (2007).
	This is the general K-means clustering algorithm of Asgharbeygi and Maleki (2008, Algorithm 1). A clustering gamma is a list of lists of integers 1, ..., n, where n is the number of rotations. Each integer appears exactly once in gamma. Each list in gamma is non-empty.

Clear[rotationDistance, rotationDiameter, rotationMean, rotationVariance, rotationGeodesicSegmentSteps, rotationGeodesicSteps, rotationClustersWGD, rotationClustersGamma, rotationClusters];

(* Returns the distance in SO(3) between the given rotations. *)
(* That's the angle of the smallest rotation taking one to the other. *)
rotationDistance[r_, q_] :=
	arcCos[(Tr[Transpose[r].q] - 1) / 2];

(* Given a set of rotations, returns the distance between the two most *)
(* distant ones. This function is not fast. *)
rotationDiameter[rots_] :=
	Module[{diameter, i, j},
		diameter = 0;
		For[i = 1, i ≤ Length[rots], i += 1,
			For[j = i + 1, j ≤ Length[rots], j += 1,
				diameter = Max[diameter, rotationDistance[rots[[i]], rots[[j]]]]]];
		diameter];

(* Computes the Fre'chet (geodesic L^2) mean of the given set of *)
(* rotations. epsilon is an error tolerance; 10^-10 seems to work *)
(* well. n is a bound on the number of iterations used; 100 *)
(* works, although I've never seen more than 14 used. *)
(* Returns a triple: mean, error achieved, iterations used. *)
(* Works best if rotations are contained within an open ball of *)
(* radius pi / 2, which is a stronger condition than the *)
(* rotations' having diameter < pi. *)
rotationMean[rots_, epsilon_, n_, seed_] :=
	Module[{mean, errorSquared, k, r},
		mean = seed;
		k = 0;
		r = Mean[Table[rotationLog[Transpose[mean].rots[[i]]], {i, 1, Length[rots]}]];
		errorSquared = frobenius[r, r];
		While[errorSquared ≥ epsilon^2 && k < n,
			mean = mean.rotationExp[r];
			r = Mean[Table[rotationLog[Transpose[mean].rots[[i]]], {i, 1, Length[rots]}]];
			errorSquared = frobenius[r, r];
			k += 1];
		{mean, Sqrt[errorSquared], k}];
rotationMean[rots_, epsilon_, n_] :=
	rotationMean[rots, epsilon, n, rots[[1]]];
rotationMean[rots_] :=
	rotationMean[rots, 10^-10, 100];

(* Returns the variance of a set of rotations about any given center. *)
(* Usually the center you want is the mean given by rotationMean. *)
(* The variance is half of the average of the distances-squared; *)
(* see Krakowski et al. (2007). *)
rotationVariance[rots_, center_] :=
	Mean[Table[rotationDistance[center, rots[[i]]]^2, {i, 1, Length[rots]}]] / 2;

(* Approximates the geodesic from rotation q to rotation r as *)
(* a sequence of rotations, using the indicated step sizes. *)
(* The sequence always begins with q and ends with r. *)
rotationGeodesicSegmentSteps[q_, r_, step_] :=
	Module[{aa, angle},
		aa = angleAxisFromRotation[r.Transpose[q]];
		Append[Table[
			rotationFromAngleAxis[{angle, aa[[2]]}].q,
			{angle, 0, aa[[1]], step}], r]];

(* Approximates the geodesic through q in the direction of the *)
(* rotation axis u (a unit vector) as a sequence of rotations, *)
(* using the indicated step size, going the full 2 pi distance *)
(* around SO(3). *)
rotationGeodesicSteps[q_, u_, step_] :=
	Append[Table[
		rotationFromAngleAxis[{angle, u}].q,
		{angle, 0, 2 Pi, step}], q];

(* Implements the W_GD general distance loss function, based *)
(* on Eq. (3). Combining Theorem 3.3 with Eq. (5) leads to *)
(* the same W_GD, but multiplied by a factor of 2. *)
rotationClustersWGD[dd_, gamma_] :=
	Sum[dd[[gamma[[l, i]], gamma[[l, j]]]],
		{l, 1, Length[gamma]}, {i, 1, Length[gamma[[l]]]}, {j, 1, Length[gamma[[l]]]}];
(* Given an old gamma, returns the new gamma. *)
rotationClustersGamma[dd_, gamma_] :=
	Module[{seconds, l, q, r, i, fulls, argmins, result},
		(* Compute the second term once for each cluster. *)
		seconds = Table[Sum[
				dd[[gamma[[l, q]], gamma[[l, r]]]],
				{q, 1, Length[gamma[[l]]]}, {r, 1, Length[gamma[[l]]]}] / Length[gamma[[l]]]^2,
			{l, 1, Length[gamma]}];
		(* Compute a full table indexed over i and then l. *)
		fulls = Table[Sum[
				dd[[i, gamma[[l, r]]]],
				{r, 1, Length[gamma[[l]]]}] * 2 / Length[gamma[[l]]] - seconds[[l]],
			{i, 1, Length[dd]}, {l, 1, Length[gamma]}];
		(* Sort each i into the cluster with the least fulls-value. *)
		argmins= Map[Function[list, argMin[list][[1]]], fulls];
		result = Table[{}, {Length[gamma]}];
		For[i = 1, i ≤ Length[dd], i += 1,
			result[[argmins[[i]]]] = Append[result[[argmins[[i]]]], i]];
		result];
(* Given a set of rotations and a starting clustering (with no cluster empty). *)
(* Also given termination conditions: epsilon for change in fit, bound on *)
(* number of iterations. Returns a triple containing the supposedly optimal *)
(* clustering, the error achieved, and the number of iterations used. *)
rotationClusters[rots_, gamma0_, epsilon_, bound_] :=
	Module[{dd, t, gamma, wgd0, wgd1, error},
		(* Build a table of distances squared between the rotations *)
		dd = Table[rotationDistance[rots[[i]], rots[[j]]]^2, {i, 1, Length[rots]}, {j, 1, Length[rots]}];
		(* Initialize the loop. *)
		gamma = gamma0;
		wgd1 = rotationClustersWGD[dd, gamma];
		error = epsilon + 1;
		t = 0;
		While[t < bound && error ≥ epsilon,
			(* Prepare for this iteration. *)
			t += 1;
			wgd0 = wgd1;
			(* Update the clustering. *)
			gamma = rotationClustersGamma[dd, gamma];
			(* Test the loop termination condition. *)
			wgd1 = rotationClustersWGD[dd, gamma];
			error = Abs[wgd1 - wgd0]];
		{gamma, error, t}];

Projected arithmetic mean and related computations

Clear[rotationProjected, rotationProjectedMean, rotationSignedSVD, rotationFisherNormalizer, rotationFisherMeanConcentration, rotationFisher3By2PValue, rotationPrenticeMeanRegion, rotationUniformRayleigh, rotationUniformRayleighPValue];

(* Projects the given matrix M onto SO(3); returns the nearest rotation to M. *)
(* May produce an answer with det = -1 if M is not close to SO(3). *)
rotationProjected[m_] :=
	Module[{vals, vecs, q, dSqrt, mHat},
		(* Q D^(1/2) Q^T = (M^T M)^(1/2). *)
		{vals, vecs} = Eigensystem[Transpose[m].m];
		dSqrt = DiagonalMatrix[Sqrt[vals]];
		q = Transpose[vecs];
		(* Mhat = R Q D^(-1/2) Q^T. *)
		mHat = m.q.Inverse[dSqrt].Transpose[q];
		mHat];

(* Computes the projected arithmetic mean of the given rotations. *)
rotationProjectedMean[rots_] :=
	rotationProjected[Mean[rots]];

(* Given 3x3 matrix M, returns a particular singular value *)
(* decomposition {Q, D, R} such that M = Q D R^T, *)
(* det Q == det R == 1, and |D11| ≥ D22 ≥ D33 ≥ 0. *)
(* See Sei et al. (2013). *)
rotationSignedSVD[m_] :=
	Module[{q, d, r, order, p},
		(* Get a singular value decomposition M = Q D R. *)
		{q, d, r} = SingularValueDecomposition[m];
		r = Transpose[r];
		(* Ensure that the singular values are in decreasing order. *)
		order = Reverse[Ordering[{d[[1, 1]], d[[2, 2]], d[[3, 3]]}]];
		p= Map[Function[k, If[k == 1, {1, 0, 0}, If[k == 2, {0, 1, 0}, {0, 0, 1}]]], order];
		q = q.Transpose[p];
		d = p.d.Transpose[p];
		r = p.r;
		(* Ensure that det Q = det R = 1. *)
		If[Det[q] < 0,
			q[[1, 1]] = -q[[1, 1]];
			q[[2, 1]] = -q[[2, 1]];
			q[[3, 1]] = -q[[3, 1]];
			d[[1, 1]] = -d[[1, 1]]];
		If[Det[r] < 0,
			r[[1]] = -r[[1]];
			d[[1, 1]] = -d[[1, 1]]];
		(* Return Q, D, R such that M = Q D R^T as usual. *)
		{q, d,Transpose[r]}];

(* Given phi1, phi2, phi3 ≥ 0, the eigenvalues of K, computes *)
(* the integral C of exp(K M^T R) with respect to the uniform *)
(* distribution. That is, C^-1 exp(K M^T R) is normalized. *)
(* Uses the one-dimensional integral of Wood (1993), with the *)
(* extraneous Pi^2 factor removed. *)
rotationFisherNormalizer[phi1_, phi2_, phi3_] :=
	Module[{p},
		p = Reverse[Sort[{phi1, phi2, phi3}]];
		NIntegrate[
			BesselI[0, (p[[1]] - p[[2]]) u] BesselI[0, (p[[1]] + p[[2]]) (1 - u)] Exp[p[[3]] (1 - 2 u)],
			{u, 0, 1}]];

(* Returns the maximum likelihood estimates Mhat and Khat *)
(* of the matrix Fisher distribution with parameters F = M K. *)
(* Mhat equals the projected mean of rotationProjectedMean. *)
rotationFisherMeanConcentration[rots_] :=
	Module[{q, g, v, f, phi1, phi2, phi3, sol, mHat, kHat},
		{q, g, v} = rotationSignedSVD[Mean[rots]];
		f[phi1_?NumericQ, phi2_?NumericQ, phi3_?NumericQ] :=
			phi1 g[[1, 1]] + phi2 g[[2, 2]] + phi3 g[[3, 3]] - Log[rotationFisherNormalizer[phi1, phi2, phi3]];
		sol = NMaximize[f[phi1, phi2, phi3], {phi1, phi2, phi3}];
		mHat = q.Transpose[v];
		kHat = v.DiagonalMatrix[{phi1, phi2, phi3} /. sol[[2]]].Transpose[v];
		{mHat, kHat}];

(* WARNING: MAY BE WRONG; R VERSION IS BETTER CHECKED. *)
(* Given 3x2 rotation matrices (rotations without their last columns). *)
(* And given a hypothesized mean 3x2 rotation matrix M0. *)
(* Returns the approximate p-value for M0 according to Downs (1972). *)
rotationFisher3By2PValue[rot3By2s_, m03By2_] :=
	Module[{n, yBar, hHat, mat, top, bottom, xi},
		n = Length[rot3By2s];
		yBar = Mean[rot3By2s];
		hHat = ellipsoidPower[Transpose[yBar].yBar, 1 / 2];
		top = 2 (IdentityMatrix[2] - hHat);
		bottom = 2 IdentityMatrix[2] - Transpose[yBar].m03By2 - Transpose[m03By2].yBar;
		xi = Sqrt[frobenius[top, top] / frobenius[bottom, bottom]];
		1 - CDF[FRatioDistribution[3, 3 n - 5], (xi^(-1 / 2) - 1) (n - 5 / 3)]];

(* Given rotation matrices, returns either a number B or a triple *)
(* {Mhat, g, F}. B is the bipolarity. If B ≤ 1/2, then the sample *)
(* is not like a cluster, so the mean is nonsensical, and we return *)
(* only B. If B > 1/2, then {M, g, F} is returned. M is the same *)
(* projected arithmetic mean from rotationProjectedMean, etc. g is *)
(* a metric on the space of quaternions. The (1 - alpha) confidence *)
(* region is those q such that q.g.q ≤ F. The p-value for a given q *)
(* is 1.0 - CDF[ChiSquareDistribution[3], q.g.q]. Based on Prentice *)
(* (1986). *)
rotationPrenticeMeanRegion[rots_, alpha_] :=
	Module[{n, xs, xxTs, i, eig, ls, aHats, aHat, mHat, mean, x2s, cHat, aHat4, fHat4, metric, f},
		n = Length[rots];
		(* Convert to quaternions in unusual order. *)
		xs = Map[quaternionFromRotation, rots];
		xs = Map[Function[x, {x[[2]], x[[3]], x[[4]], x[[1]]}], xs];
		(* Compute the lambda_i and A-hat. *)
		xxTs = Table[Outer[Times, xs[[i]], xs[[i]]], {i, 1, n}];
		eig = Transpose[Eigensystem[Mean[xxTs]]];
		{ls, aHats} = Transpose[Sort[eig, Function[{x, y}, x[[1]] < y[[1]]]]];
		If[ls[[1]] + ls[[4]] ≤ 1 / 2,
			ls[[1]] + ls[[4]],
			If[Det[aHats] < 0,
				aHats[[4]] = -aHats[[4]]];
			aHat = Transpose[aHats];
			(* Estimate the center M of the distribution. *)
			mHat = aHats[[4]];
			mean = rotationFromQuaternion[{mHat[[4]], mHat[[1]], mHat[[2]], mHat[[3]]}];
			(* Approximate the confidence region metric. *)
			x2s = Table[Transpose[aHat].xxTs[[i]].aHat, {i, 1, n}];
			x2s = Map[Function[m, {m[[1, 1]], m[[2, 2]], m[[3, 3]], m[[4, 4]]}], x2s];
			cHat = Sum[Outer[Times, x2s[[i]], x2s[[i]]], {i, 1, n}] / n;
			aHat4 = Transpose[aHats[[1;;3]]];
			fHat4 = aHat4.DiagonalMatrix[{
				(ls[[4]] - ls[[1]])^2 / cHat[[4, 1]],
				(ls[[4]] - ls[[2]])^2 / cHat[[4, 2]],
				(ls[[4]] - ls[[3]])^2 / cHat[[4, 3]]}].Transpose[aHat4];
			(* Put the metric back into my preferred quaternion order. *)
			metric = n {
				{fHat4[[4, 4]], fHat4[[4, 1]], fHat4[[4, 2]], fHat4[[4, 3]]},
				{fHat4[[1, 4]], fHat4[[1, 1]], fHat4[[1, 2]], fHat4[[1, 3]]},
				{fHat4[[2, 4]], fHat4[[2, 1]], fHat4[[2, 2]], fHat4[[2, 3]]},
				{fHat4[[3, 4]], fHat4[[3, 1]], fHat4[[3, 2]], fHat4[[3, 3]]}};
			f = InverseCDF[ChiSquareDistribution[3], 1 - alpha];
			{mean, metric, f}]];

(* Computes the Rayleigh statistic S (Sei et al., 2013). *)
(* Then compares it to the chi-squared distribution with *)
(* 9 degrees of freedom. Low P-values suggest that the *)
(* sample is not from the uniform distribution. *)
rotationUniformRayleigh[rots_] :=
	Module[{n, xBar, s},
		n = Length[rots];
		xBar = Mean[rots];
		s = 3 n Tr[Transpose[xBar].xBar];
		(*sStar = s (1 - (1 - s / 8) / (2 n));*)
		s];
rotationUniformRayleighPValue[rots_] :=
	1.0 - CDF[ChiSquareDistribution[9], rotationUniformRayleigh[rots]];

Tangent space methods

Clear[tangentVectorFromRotation, rotationFromTangentVector, rotationCovariance, rotationPrincipalComponents, rotationMahalanobisNorm, rotationBestGeodesic, rancourtVectorFromRotation, rotationFromRancourtVector, rotationRancourtMeanRegion];

(* Returns the rotation as a 3-vector v in the tangent space at *)
(* the center. That is, the vector is v = alpha * u expressing *)
(* Transpose[center].rotation. That is, rot = center.exp[matrix[v]]. *)
tangentVectorFromRotation[rot_, center_] :=
	Module[{log},
		log = rotationLog[Transpose[center].rot];
		{log[[3, 2]], log[[1, 3]], log[[2, 1]]}];
rotationFromTangentVector[v_, center_] :=
	center.rotationExp[{{0, -v[[3]], v[[2]]}, {v[[3]], 0, -v[[1]]}, {-v[[2]], v[[1]], 0}}];

(* Returns the MLE of the covariance matrix of the rotations. *)
(* center should be the rotations' mean, from rotationMean. *)
rotationCovariance[rots_, center_] :=
	Module[{n, vecs},
		n = Length[rots];
		vecs = Table[tangentVectorFromRotation[rots[[i]], center], {i, 1, n}];
		Mean[Table[Outer[Times, vecs[[i]], vecs[[i]]], {i, 1, n}]]];

(* Given a set of n rotations and their center (from rotationMean). *)
(* Returns a set of n decreasing nonnegative numbers and their *)
(* corresponding n unit vectors. The vectors are rotation axes, *)
(* describing how to move away from the center in each principal *)
(* direction of the sample. The numbers are like standard deviations. *)
rotationPrincipalComponents[rots_, center_]:=
	Module[{eig},
		eig = Eigensystem[rotationCovariance[rots, center]];
		{Sqrt[eig[[1]]], Map[unit, eig[[2]]]}];

(* Given a set of rotations and their center (from rotationMean). *)
(* Given inverse of covariance matrix (from rotationCovariance). *)
(* Returns a list of their Mahalanobis norms. For example, the *)
(* rotations at the ends of the principal geodesics have norm 1. *)
rotationMahalanobisNorm[rot_, center_, covarInv_] :=
	Module[{vec},
		vec = tangentVectorFromRotation[rot, center];
		Sqrt[vec.covarInv.vec]];

(* Given a set of rotations and their center (from rotationMean). *)
(* Returns a unit vector u, such that the best-fit geodesic through *)
(* the rotations is rotationFromAngleAxis[{t, u}].center. *)
rotationBestGeodesic[rots_, center_]:=
	rotationPrincipalComponents[rots, center][[2, 1]];

(* Given a rotation and a center (usually rotationProjectedMean), *)
(* converts to a vector in the tangent space at the center. *)
rancourtVectorFromRotation[rot_, center_] :=
	Module[{v},
		v = (Transpose[center].rot - Transpose[rot].center) / 2;
		{v[[3, 2]], v[[1, 3]], v[[2, 1]]}];
rotationFromRancourtVector[v_, center_] :=
	Module[{cosine, vv},
		cosine = Sqrt[1 - v.v];
		vv = {{0, -v[[3]], v[[2]]}, {v[[3]], 0, -v[[1]]}, {-v[[2]], v[[1]], 0}};
		center.(IdentityMatrix[3] + vv + vv.vv / (1 + cosine))];
(* Confidence region for the projected arithmetic mean. Assumes *)
(* data tightly clustered, then works in tangent space. Typically *)
(* alpha = 0.05 or 0.01. Returns triple {Mhat, g, f} such that *)
(* Mhat is the projected arithmetic mean and the (1 - alpha) *)
(* confidence region is all matrices *)
(* rotationFromRancourtVector[v, Mhat] such that v.g.v < f. *)
(* The p-value for a given v is *)
(* 1.0 - CDF[FRatioDistribution[3, n - 3], v.g.v]. *)
(* Based on Rancourt et al. (2000), but that paper is missing a *)
(* factor of n or n - 1 in its Eq. (2.6)? *)
rotationRancourtMeanRegion[rots_, alpha_] :=
	Module[{n, q, d, r,mHat, vs, s, g, f},
		n = Length[rots];
		{q, d, r} = rotationSignedSVD[Mean[rots]];
		mHat = q.Transpose[r];
		vs = Map[Function[rot, rancourtVectorFromRotation[rot, mHat]], rots];
		s = Sum[Outer[Times, vs[[i]], vs[[i]]], {i, 1,n}] / (n - 1);
		g = n (n - 3) / (3 (n - 1))Inverse[s];
		f = InverseCDF[FRatioDistribution[3, n - 3], 1 - alpha];
		{mHat, g, f}];

Sampling from distributions

Sampling from various distributions on SO(3). Testing the uniform distribution.

Clear[rotationUniform,  rotationWrappedTrivariateNormal, rotationIsotropicFisher, rotationFisher];

(* Returns a uniformly chosen rotation, following the subgroup *)
(* algorithm of Diaconis and Shahshahani (1987). *)
rotationUniform[] :=
	Module[{v, x, refl, a},
		v = randomSpherePoint[];
		x = unit[{1, 0, 0} - v];
		refl = IdentityMatrix[3] - 2 Outer[Times, x, x];
		a = RandomReal[{0, 2 Pi}];
		-refl.{{1, 0, 0}, {0, Cos[a], -Sin[a]}, {0, Sin[a], Cos[a]}}];
rotationUniform[n_] :=
	Table[rotationUniform[], {n}];

(* Returns a rotation drawn from the wrapped trivariate normal *)
(* distribution (e.g., Qiu et al., 2012) with center s and *)
(* concentration kappa. *)
rotationWrappedTrivariateNormal[s_, kappa_] :=
	Module[{v1, v2, v3},
		v1 = RandomVariate[NormalDistribution[0, kappa^-1]];
		v2 = RandomVariate[NormalDistribution[0, kappa^-1]];
		v3 = RandomVariate[NormalDistribution[0, kappa^-1]];
		s.rotationExp[{{0, -v3, v2}, {v3, 0, -v1}, {-v2, v1, 0}}]];
rotationWrappedTrivariateNormal[s_, kappa_, n_] :=
	Table[rotationWrappedTrivariateNormal[s, kappa], {n}];

(* Samples from the matrix Fisher distribution with center M *)
(* (3x3 rotation matrix) and concentration K (3x3 symmetric *)
(* positive-definite). Uses the rejection sampling algorithm of *)
(* Kent et al. (2013), "A new method to simulate the Bingham...". *)
(* The efficiency seems to be about 45%: 2.2 tries per success. *)
rotationFisher[lambda_, omega_, sigma_, bound_] :=
	Module[{x, w},
		x = unit[RandomVariate[MultinormalDistribution[{0, 0, 0, 0}, sigma]]];
		w = RandomReal[{0, 1}];
		(* Repeat until w < f*(x) / (m* g*(x)). *)
		While[bound w ≥ Exp[-x.lambda.x] (x.omega.x)^2,
			x = unit[RandomVariate[MultinormalDistribution[{0, 0, 0, 0}, sigma]]];
			w = RandomReal[{0, 1}]];
		rotationFromQuaternion[x]];
rotationFisher[m_, k_, n_] :=
	Module[{d, vecs, v, l, b, sol, bound, omega, sigma, xs},
		(* Diagonalize K as V Delta V^T. *)
		{d, vecs} = Eigensystem[k];
		v = Transpose[vecs];
		(* Build the Lambda matrix for the Bingham distribution. *)
		l = DiagonalMatrix[{0, 2 (d[[2]] + d[[3]]), 2 (d[[1]] + d[[3]]), 2 (d[[1]] + d[[2]])}];
		(* Compute b and the bound M*. *)
		sol= Solve[1 == 1 / (b) + 1 / (b + 2 l[[2, 2]]) + 1 / (b + 2 l[[3, 3]]) + 1 / (b + 2 l[[4, 4]]), {b}, Reals];
		b = Max[sol[[All, 1, 2]]];
		bound = Exp[(b - 4) / 2] (4 / b)^2;
		(* Compute Omega for ACG and Sigma = Omega^-1 for normal. *)
		omega = IdentityMatrix[4] + 2 l / b;
		sigma = Inverse[omega];
		xs = Table[m.v.rotationFisher[l, omega, sigma, bound].Transpose[v], {n}];
		xs];
rotationFisher[m_, k_] :=
	rotationFisher[m, k, 1][[1]];

(* Isotropic matrix Fisher distribution. See Qiu et al. (2012). *)
(* This function uses ordinary matrix Fisher with K = (kappa^2 / 2) I. *)
rotationIsotropicFisher[s_, kappa_, n_] :=
	rotationFisher[s, DiagonalMatrix[(kappa^2 / 2){1, 1, 1}], n];
rotationIsotropicFisher[s_, kappa_] :=
	rotationIsotropicFisher[s, kappa, 1][[1]];

Regression

Clear[rotationRegressionCurve, rotationRegressionSurface, rotationRegressionBandwidth, rotationRegressionBandwidthVerbose, rotationRegressionKernel];

(* Given pairs of the form {x, rot}, where x is a number, returns *)
(* {m, b, e}, where exp(x m).b is the best-fit geodesic and e is *)
(* the error squared in the fit. To get the R^2 statistic, compute *)
(* 1 - e/var, where var is the variance of the rots about their mean. *)
rotationRegressionCurve[xrs_] :=
	Module[{n, q, e, m32, m13, m21, w32, w13, w21, min, w, m},
		(* Let q be the r whose x is closest to zero. *)
		q = xrs[[argMin[Abs[xrs[[All, 1]]]][[1]], 2]];
		n = Length[xrs];
		e[m32_, m13_, m21_, w32_, w13_, w21_] :=
			Module[{mm, b, i},
				b = rotationExp[{{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}}].q;
				mm = {{0, -m21, m13}, {m21, 0, -m32}, {-m13, m32, 0}};
				Sum[rotationDistance[xrs[[i, 2]], rotationExp[xrs[[i, 1]] mm].b]^2, {i, 1, n}] / (2 n)];
		min = Minimize[
			{e[m32, m13, m21, w32, w13, w21],
			m32^2 + m13^2 + m21^2 ≤ Pi^2,
			w32^2 + w13^2 + w21^2 ≤ Pi^2},
			{m32, m13, m21, w32, w13, w21}];
		w = {{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}} /. min[[2]];
		m = {{0, -m21, m13}, {m21, 0, -m32}, {-m13, m32, 0}} /. min[[2]];
		{m, rotationExp[w].q, min[[1]]}];

(* Given triples of the form {x, y, rot}, where x and y are numbers, *)
(* returns {m, l, b, e}, where exp(x m + y l).b is the best-fit geodesic *)
(* and e is the error squared in the fit. To get the R^2 statistic, compute *)
(* 1 - e/var, where var is the variance of the rots about their mean. *)
rotationRegressionSurface[xyrs_] :=
	Module[{norms, n, q, e, m32, m13, m21, l32, l13, l21, w32, w13, w21, min, w, m, l},
		(* Let q be the r whose (x, y) is closest to (0, 0). *)
		norms = Map[Norm, xyrs[[All, 1;;2]]];
		q = xyrs[[argMin[norms][[1]], 3]];
		n = Length[xyrs];
		e[m32_, m13_, m21_, l32_, l13_, l21_, w32_, w13_, w21_] :=
			Module[{mm, ll, b, i},
				b = rotationExp[{{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}}].q;
				mm = {{0, -m21, m13}, {m21, 0, -m32}, {-m13, m32, 0}};
				ll = {{0, -l21, l13}, {l21, 0, -l32}, {-l13, l32, 0}};
				Sum[rotationDistance[
						xyrs[[i, 3]],
						rotationExp[xyrs[[i, 1]] mm + xyrs[[i, 2]] ll].b]^2,
					{i, 1, n}] / (2 n)];
		min = Minimize[
			{e[m32, m13, m21, l32, l13, l21, w32, w13, w21],
			m32^2 + m13^2 + m21^2 ≤ Pi^2,
			l32^2 + l13^2 + l21^2 ≤ Pi^2,
			w32^2 + w13^2 + w21^2 ≤ Pi^2},
			{m32, m13, m21, l32, l13, l21, w32, w13, w21}];
		w = {{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}} /. min[[2]];
		m = {{0, -m21, m13}, {m21, 0, -m32}, {-m13, m32, 0}} /. min[[2]];
		l = {{0, -l21, l13}, {l21, 0, -l32}, {-l13, l32, 0}} /. min[[2]];
		{m, l, rotationExp[w].q, min[[1]]}];

(* See rotationRegressionNonparametric for input types. *)
(* Returns the bandwidth determined by cross-validation. *)
(* Warning: This function takes orders of magnitude more time *)
(* than rotationRegressionKernel does. If you can't *)
(* wait, then manually pick a bandwidth, such as 0.09. *)
rotationRegressionBandwidth[xrs_, k_] :=
	Module[{f, h, sol},
		f[h_?NumericQ] :=
			Module[{rots, val},
				rots = Table[rotationRegressionKernel[
					Join[xrs[[;;(i - 1)]], xrs[[(i + 1);;]]],
					{xrs[[i, 1]]}, k, h][[1]], {i, 1, Length[xrs]}];
				val = Sum[rotationDistance[xrs[[i, 2]], rots[[i]]]^2, {i, 1, Length[xrs]}];
				val];
		sol = NMinimize[{f[h], h > 0}, {h}];
		h /. sol[[2]]];
rotationRegressionBandwidthVerbose[xrs_, k_] :=
	Module[{best, count, f, h, sol},
		best = 10^8;
		count = 0;
		f[h_?NumericQ] :=
			Module[{rots, val},
				count += 1;
				rots = Table[rotationRegressionKernel[
					Join[xrs[[;;(i - 1)]], xrs[[(i + 1);;]]],
					{xrs[[i, 1]]}, k, h][[1]], {i, 1, Length[xrs]}];
				val = Sum[rotationDistance[xrs[[i, 2]], rots[[i]]]^2, {i, 1, Length[xrs]}];
				If[val < best,
					best = val;
					Print["best = ", best, ", h = ", h, ", count = ", count, ", time = ", DateString[]]];
				val];
		sol = NMinimize[{f[h], h > 0}, {h}];
		h /. sol[[2]]];

(* Given a list of pairs, each of the form {x, rot}. *)
(* Given a list of x-values, at which to plot the regression curve. *)
(* Given a kernel, such as PDF[NormalDistribution[0, 1]]. *)
(* Given a bandwidth, as produced by rotationRegressionBandwidth. *)
(* Returns a list of rotations defining the regression curve. *)
rotationRegressionKernel[xrs_, xs_, k_, h_] :=
	Module[{n, sols, w32, w13, w21},
		n = Length[xrs];
		sols = Table[Minimize[{
			Sum[k[(xs[[j]] - xrs[[i, 1]]) / h] / h rotationDistance[
					xrs[[i, 2]],
					rotationExp[{{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}}]],
				{i, 1, n}] / Sum[k[(xs[[j]] - xrs[[i, 1]]) / h] / h, {i, 1, n}],
			w32^2 + w13^2 + w21^2 ≤ Pi}, {w32, w13, w21}],
			{j, 1, Length[xs]}];
		Table[
			rotationExp[{{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}} /. sols[[i, 2]]],
			{i, 1, Length[sols]}]];

Angle-axis and equal-volume plots

The angle-axis plot shows the angle-axis rotation (a, u) at the point a u. It does not respect volume. The equal-volume plot is a radially distorted version of the angle-axis plot, that respects volume. The equal-angle plot is not here. These functions include style arguments. Typical values are Automatic, PointSize[0.02], Thick, etc.

Clear[ballPlotCurveSegments, ballPlot, angleAxisPlot, volumetricPlot, rotationPlot];

ballPlotCurveSegments[c_, r_, t_]:=
	Module[{res, i, norm, bounce},
		res = {};
		For[i = 1, i < Length[c], i += 1,
			norm = Norm[c[[i]] - c[[i + 1]]];
			If[norm < 1.9 r,
				res = Append[res, c[[i]] + t (c[[i + 1]] - c[[i]])],
				bounce = (c[[i]] - c[[i + 1]]) r / norm;
				res = Append[res, c[[i]] + t (bounce - c[[i]])];
				res = Append[res, -bounce + t (c[[i + 1]] + bounce)]]];
		res];

ballPlot[points_, pointStyles_, curves_, curveStyles_, r_] :=
	Module[{x, y, z, e, container, axes, scatter, splits, t, curveStylesNew},
		e = 10^-4;
		container = ContourPlot3D[
			x^2 + y^2 + z^2 == r^2,
			{x, -r, r}, {y, -r, r}, {z, -r, r},
			Mesh->None, ContourStyle->Directive[Opacity[0.1]],
			BoxRatios->{1, 1, 1}, Boxed->False, Axes->False];
		axes = ParametricPlot3D[
			{{t, 0, 0}, {0, t, 0}, {0, 0, t}}, {t, 0, r},
			PlotStyle->{{Thick, Red}, {Thick, Green}, {Thick, Blue}}, Boxed->False, Axes->False];
		scatter = ListPointPlot3D[
			points,
			PlotRange->{{-r - e, r + e}, {-r - e, r + e}, {-r - e, r + e}},
			PlotStyle-> pointStyles,
			BoxRatios->{1, 1, 1}, Boxed->False, Axes->False];
		splits =Table[ballPlotCurveSegments[curves[[i]], r, t], {i, 1, Length[curves]}];
		curveStylesNew = Flatten[Table[
			Table[curveStyles[[i]], {Length[splits[[i]]]}],
			{i, 1, Length[curveStyles]}], 1];
		splits = Flatten[splits, 1];
		splits = ParametricPlot3D[
			splits, {t, 0, 1},
			PlotStyle->curveStylesNew, Boxed->False, Axes->False];
		Show[container, axes, scatter, splits, Boxed->False, Axes->False]];

(* Non-equal-volume plot of rotations as angle-axis vectors, *)
(* in a ball of radius pi. That is, (a, u) plots at point a u. *)
angleAxisPlot[aas_, pointStyles_, aaCurves_, curveStyles_] :=
	ballPlot[
		aas[[All, 1]] aas[[All, 2]], pointStyles,
		Table[aaCurves[[i, All, 1]] aaCurves[[i, All, 2]], {i, 1, Length[aaCurves]}],
		curveStyles, Pi];
(* Equal-volume plot of rotations and volumetrics. *)
(* Draws SO(3) as a ball of radius (3 / (4 pi))^(1 / 3). *)
volumetricPlot[vols_, pointStyles_, volCurves_, curveStyles_] :=
	ballPlot[
		vols, pointStyles,
		volCurves, curveStyles,
		(3 / (4 Pi))^(1 / 3)];
volumetricPlot[vols_] :=
	volumetricPlot[vols, Automatic, {}, Automatic];
rotationPlot[rots_, pointStyles_, rotCurves_, curveStyles_] :=
	ballPlot[
		(* This first Map call operates on rotation matrices, *)
		(* no matter how deeply nested in the list they are. *)
		Map[volumetricFromRotation, rots, {-3}], pointStyles,
		Map[volumetricFromRotation, rotCurves, {2}], curveStyles,
		(3 / (4 Pi))^(1 / 3)];
rotationPlot[rots_] :=
	rotationPlot[rots, Automatic, {}, Automatic];

Plotting functions

Equal-volume contour plots of arbitrary functions f : SO(3) -> R. For enhanced versatility near the boundaries of the equal-volume plot, the input to f is given as a volumetric vector, rather than as a rotation matrix.

Clear[volumetricLevelBase, volumetricMidpoint, volumetricLevelEdgesFacesCenter, 
	volumetricLevelRelevants, volumetricLevelBoxes, volumetricLevelCut, volumetricLevelPolygon, 
	volumetricLevelSurfaces, volumetricLevelPlot];

(* Returns a list of seven boxes, each a list of eight volumetric vectors. *)
volumetricLevelBase[] :=
	Module[{pnp, ppp, nnn, nnp, npn, npp, pnn, ppn, up, north, east, south, west, down, center},
		(* Make the eight outer vertices. *)
		pnp = volumetricFromRotation[
			{{-0.3333333333101354`,-0.6666666666724664`,0.6666666666724665`},
			{-0.6666666666724662`,-0.3333333333449323`,-0.6666666666550678`},
			{0.6666666666724662`,-0.6666666666550678`,-0.3333333333449322`}}];
		ppp = {1, -1, 1} pnp;
		nnn = {-1, -1, -1} ppp;
		nnp = {-1, -1, 1} ppp;
		npn = {-1, 1, -1} ppp;
		npp = {-1, 1, 1} ppp;
		pnn = {1, -1, -1} ppp;
		ppn = {1, 1, -1} ppp;
		(* Make the seven boxes. *)
		up = Join[{nnp, npp, ppp, pnp} / 2, {nnp, npp, ppp, pnp}];
		north = Join[{ppn, npn, npp, ppp} / 2, {ppn, npn, npp, ppp}];
		east = Join[{pnn, ppn, ppp, pnp} / 2, {pnn, ppn, ppp, pnp}];
		south = Join[{nnn, pnn, pnp, nnp} / 2, {nnn, pnn, pnp, nnp} ];
		west = Join[{npn, nnn, nnp, npp} / 2, {npn, nnn, nnp, npp}];
		down = Join[{ppn, pnn, nnn, npn} / 2, {ppn, pnn, nnn, npn}];
		center = Join[{nnn, pnn, ppn, npn, nnp, pnp, ppp, npp} / 2];
		{down, up, south, east, north, west, center}];

(* Returns the volumetric vector halfway between the two given. *)
volumetricMidpoint[vol1_, vol2_] :=
	Module[{plotRadius, rot1, rot2, rot},
		(* plotRadius = (3 / (4 Pi))^(1 / 3); *)
		plotRadius = 0.6203504908994`;
		If[Abs[Norm[vol1] - plotRadius] < 10^-8 && Abs[Norm[vol2] - plotRadius] < 10^-8,
			(* The rotations are at the boundary of the equal-volume plot. *)
			plotRadius (vol1 + vol2) / Norm[vol1 + vol2],
			(* The rotations are not at the boundary of the plot. *)
			rot1 = rotationFromVolumetric[vol1];
			rot2 = rotationFromVolumetric[vol2];
			rot = rotationExp[rotationLog[rot2.Transpose[rot1]] / 2].rot1;
			volumetricFromRotation[rot]]];

(* Given a box, returns the edge, face, and box midpoints. *)
(* Attaches the value of f to each of these points. *)
volumetricLevelEdgesFacesCenter[f_, box_] :=
	Module[{edges, faces, center, i},
		(* Find the point in the center of each edge. *)
		edges = Map[
			Function[edge,
				volumetricMidpoint[
					box[[edge[[1]], 1;;3]],
					box[[edge[[2]], 1;;3]]]],
			{{1, 2}, {2, 3}, {3, 4}, {4, 1}, {5, 6}, {6, 7},
				{7, 8}, {8, 5}, {1, 5}, {2, 6}, {3, 7}, {4, 8}}];
		(* Find the point in the center of each face. *)
		faces = Map[
			Function[face,
				volumetricMidpoint[edges[[face[[1]]]], edges[[face[[2]]]]]],
			{{1, 3}, {5, 7}, {1, 5}, {2, 6}, {3, 7}, {4, 8}}];
		(* Find the point in the center of the box. *)
		center = volumetricMidpoint[faces[[1]], faces[[2]]];
		(* Compute f at each of these new mesh points. *)
		edges = Table[Append[edges[[i]], f[edges[[i]]]], {i, 1, Length[edges]}];
		faces = Table[Append[faces[[i]], f[faces[[i]]]], {i, 1, Length[faces]}];
		center = Append[center, f[center]];
		{edges, faces, center}];

(* Given a box with f-values but not relevant contours, *)
(* returns a list of all contours crossed by the box. *)
(* If minDepth ≥ 1, then all contours are trivially relevant. *)
volumetricLevelRelevants[boxWithoutCont_, contours_, minDepth_] :=
	Module[{min, max},
		If[minDepth ≥ 1,
			contours,
			min = Min[boxWithoutCont[[All, 4]]];
			max = Max[boxWithoutCont[[All, 4]]];
			Select[contours, Function[c, min ≤ c ≤ max]]]];

(* Assumes that box should be subdivided into eight boxes. *)
(* Performs the subdivision, and records which contours hit which subboxes. *)
volumetricLevelBoxes[f_, minDepth_, maxDepth_, box_] :=
	Module[{edges, faces, center, boxesWithoutCont, boxes, i, rel},
		(* Form the eight new boxes. *)
		{edges, faces, center} = volumetricLevelEdgesFacesCenter[f, box];
		boxesWithoutCont = {{box[[1]], edges[[1]], faces[[1]], edges[[4]],
				edges[[9]], faces[[3]], center, faces[[6]]},
			{edges[[1]], box[[2]], edges[[2]], faces[[1]],
				faces[[3]], edges[[10]], faces[[4]], center},
			{faces[[1]], edges[[2]], box[[3]], edges[[3]],
				center, faces[[4]], edges[[11]], faces[[5]]},
			{edges[[4]], faces[[1]], edges[[3]], box[[4]],
				faces[[6]], center, faces[[5]], edges[[12]]},
			{edges[[9]], faces[[3]], center, faces[[6]],
				box[[5]], edges[[5]], faces[[2]], edges[[8]]},
			{faces[[3]], edges[[10]], faces[[4]], center,
				edges[[5]], box[[6]], edges[[6]], faces[[2]]},
			{center, faces[[4]], edges[[11]], faces[[5]],
				faces[[2]], edges[[6]], box[[7]], edges[[7]]},
			{faces[[6]], center, faces[[5]], edges[[12]],
				edges[[8]], faces[[2]], edges[[7]], box[[8]]}};
		(* Determine relevant contours; discard irrelevant boxes. *)
		boxes = {};
		For[i = 1, i ≤ Length[boxesWithoutCont], i += 1,
			rel = volumetricLevelRelevants[
				boxesWithoutCont[[i]], box[[9]], minDepth];
			If[Length[rel] ≥ 1,
				boxes = Append[boxes, Append[boxesWithoutCont[[i]], rel]]]];
		(* Recurse on the relevant boxes, unless this was the last depth. *)
		If[maxDepth == 1,
			boxes,
			Flatten[Table[
				volumetricLevelBoxes[f, minDepth - 1, maxDepth - 1, boxes[[i]]],
				{i, 1, Length[boxes]}], 1]]];

(* Helper for volumetricContourPolygon. Cuts an edge. *)
(* Returns a volumetric vector without f value. *)
volumetricLevelCut[a_, b_, contour_] :=
	Module[{t},
		If[a[[-1]]== b[[-1]],
			If[contour == a[[-1]], a[[1;;-2]], {}],
			t = (contour - a[[-1]]) / (b[[-1]] - a[[-1]]);
			If[t < 0 || t > 1,
				{},
				a[[1;;-2]] + t (b[[1;;-2]] - a[[1;;-2]])]]];

(* Returns the polygon for a single contour within a single box. *)
volumetricLevelPolygon[box_, contour_] :=
	Module[{cuts, ij, newCuts, i, j, found, v, c, mean},
		(* Find every point where the contour cuts the 12 edges. *)
		cuts = Map[
			Function[ij, 
				volumetricLevelCut[box[[ij[[1]]]], box[[ij[[2]]]], contour]],
			{{1, 2}, {2, 3}, {3, 4}, {4, 1}, {5, 6}, {6, 7}, {7, 8}, {8, 5},
				{1, 5}, {2, 6}, {3, 7}, {4, 8}}];
		newCuts = Select[cuts, Function[c, Length[c] == 3]];
		(* Remove duplicates. *)
		cuts = {};
		For[i = 1, i ≤ Length[newCuts], i += 1,
			found = False;
			j = 1;
			While[!found && j ≤ Length[cuts],
				found = (Norm[cuts[[j]] - newCuts[[i]]] < 10^-8);
				j += 1];
			If[!found,
				cuts = Append[cuts, newCuts[[i]]]]];
		If[Length[cuts] < 3,
			{},
			(* Use the first 3 cuts to determine a projection plane. *)
			v = Cross[cuts[[1]] - cuts[[2]], cuts[[3]] - cuts[[2]]];
			i = argMin[-v * v][[1]];
			(* Tag each cut with its projection to that plane. *)
			cuts = Map[Function[c, {Join[c[[1;;i - 1]], c[[i + 1;;3]]], c}], cuts];
			(* Sort the cuts by direction of displacement from the mean. *)
			mean = Mean[cuts[[All, 1]]];
			cuts = Map[Function[c, {c[[1]] - mean, c[[2]]}], cuts];
			cuts = SortBy[cuts,Function[c, ArcTan[c[[1, 1]], c[[1, 2]]]]];
			(* The cuts are now sorted counterclockwise w.r.t. the plane. *)
			cuts[[All, 2]]]];

(* The heart of volumetricContourPlot. Returns a list of surfaces. *)
(* Each surface is a list of polygons. Each polygon is a list of *)
(* volumetric vectors, in some decent order (CW or CCW). *)
volumetricLevelSurfaces[f_, contours_, minDepth_, maxDepth_] :=
	Module[{boxes, i, surfaces, surface, j, polygon},
		(* Attach values to the base mesh points. *)
		boxes = Map[
			Function[vol, Append[vol, f[vol]]],
			volumetricLevelBase[], {2}];
		(* Attach contours to the boxes. *)
		boxes = Table[Append[boxes[[i]], contours], {i, 1, Length[boxes]}];
		(* Recursively subdivide. *)
		boxes = Flatten[Table[
			volumetricLevelBoxes[f, minDepth, maxDepth, boxes[[i]]],
			{i, 1, Length[boxes]}], 1];
		(* Build a surface for each contour. *)
		surfaces = {};
		For[i = 1, i ≤ Length[contours], i += 1,
			surface = {};
			For[j = 1, j ≤ Length[boxes], j += 1,
				If[MemberQ[boxes[[j, 9]], contours[[i]]],
					polygon = volumetricLevelPolygon[boxes[[j]], contours[[i]]];
					If[Length[polygon] ≥ 3,
						surface = Append[surface, polygon]]]];
			surfaces = Append[surfaces, surface]];
		surfaces];

(* Returns a contour plot, showing certain level surfaces of the function f. *)
(* f takes volumetric vectors as input, and produces real numbers as output. *)
(* contours is a list of values of f, where level surfaces should be drawn. *)
(* 1 ≤ minDepth ≤ maxDepth control the quality of the plot. *)
volumetricLevelPlot[f_, contours_, minDepth_, maxDepth_] :=
	Module[{surfaces},
		surfaces = volumetricLevelSurfaces[f, contours, minDepth, maxDepth];
		Graphics3D[Join[
			{EdgeForm[], Opacity[0.5]},
			Map[Polygon, Join[surfaces]]], Boxed->False]];

Kamb contouring

Kamb contouring on SO(3).

Clear[rotationKambFunctionTerm, rotationKambFunction, rotationKambPlot, rotationKambPlotC];

(* A single term of the function f used in Kamb contouring. *)
rotationKambFunctionTerm[dataRot_, rot_, r_, coeffs_] :=
	Module[{alpha},
		alpha = rotationDistance[dataRot, rot];
		If[alpha > r,
			0,
			Total[Table[coeffs[[i]] alpha^(i - 1), {i, 1, Length[coeffs]}]]]];

(* Returns the function f, such that contouring f is Kamb contouring the data. *)
rotationKambFunction[dataRots_, r_, degree_] :=
	Module[{coeffs, i},
		(* Determine the coefficients of the weighting polynomial. *)
		coeffs = If[degree == 3,
			(* This is the cubic case. *)
			{(2 r^4 - 2 r^3 Sin[r]) / (12 r Sin[r] + 24 Cos[r] + r^4 - 24), 0,
				(6 r Sin[r] - 6 r^2) / (12 r Sin[r] + 24 Cos[r] + r^4 - 24),
				(4 r - 4 Sin[r]) / (12 r Sin[r] + 24 Cos[r] + r^4 - 24)},
			If[degree == 1,
				(* This is the linear case. *)
				{(2 r^2 - 2 r Sin[r]) / (2 Cos[r] + r^2 - 2), (2 Sin[r] - 2 r) / (2 Cos[r] + r^2 - 2)},
				(* Default to no weighting. *)
				{1}]];
		Function[vol,
			Total[Table[rotationKambFunctionTerm[
					dataRots[[i]], rotationFromVolumetric[vol], r, coeffs],
			{i, 1, Length[dataRots]}]]]];

(* Returns a Kamb-contoured density plot of the given rotations. *)
(* k is the smoothing factor; 3 is a typical value. mults specifies *)
(* which multiples of the stddev should be shown as level surfaces. *)
(* degree is 0, 1, or 3 for the degree of the weighting function. *)
(* 1 ≤ minDepth ≤ maxDepth control the quality of the plot. *)
rotationKambPlot[rots_, k_, mults_, degree_, minDepth_, maxDepth_, pointStyles_] :=
	Module[{n, r, p, sigma, f, levels},
		(* Determine Kamb contouring parameters. *)
		n = Length[rots];
		r = NSolve[{rad - Sin[rad] == Pi k^2 / (n + k^2), 0 ≤ rad ≤ Pi}, rad, Reals];
		r = r[[1, 1, 2]];
		p = (r - Sin[r]) / Pi;
		sigma = Sqrt[n p (1 - p)];
		levels = sigma mults;
		f = rotationKambFunction[rots, r, degree];
		Show[
			volumetricLevelPlot[f, levels, minDepth, maxDepth],
			rotationPlot[rots, pointStyles, {}, Automatic]]];

(* Like rotationKambPlot, but just one multiple, and much much faster. *)
(* Requires georotationsKamb.c to be compiled to georotationsKamb.o *)
(* and placed in same directory as this notebook. *)
rotationKambPlotC[rots_, k_, mult_, degree_, depth_, pointStyles_] :=
	Module[{stream, r, fileToC, fileFromC, error, polygons},
		(* Write meta-data to the configuration file. *)
		stream = OpenWrite[];
		Write[stream, Length[rots]];
		Write[stream, Chop[k]];
		Write[stream, Chop[mult]];
		Write[stream, degree];
		Write[stream, depth];
		(* Write rotations to the configuration file in column-major order. *)
		Map[Function[r, 
			Write[stream,  Chop[r[[1, 1]]]];
			Write[stream,  Chop[r[[2, 1]]]];
			Write[stream,  Chop[r[[3, 1]]]];
			Write[stream,  Chop[r[[1, 2]]]];
			Write[stream,  Chop[r[[2, 2]]]];
			Write[stream,  Chop[r[[3, 2]]]];
			Write[stream,  Chop[r[[1, 3]]]];
			Write[stream,  Chop[r[[2, 3]]]];
			Write[stream,  Chop[r[[3, 3]]]]], rots];
		(* Invoke the C program on the configuration file. *)
		fileToC = Close[stream];
		stream = OpenWrite[];
		fileFromC = Close[stream];
		SetDirectory[NotebookDirectory[]];
		error = Run["./georotationsKamb.o " <> fileToC <> " " <> fileFromC];
		polygons = Get[fileFromC];
		DeleteFile[fileToC];
		DeleteFile[fileFromC];
		(* Return error string or graphics object. *)
		If[StringQ[polygons],
			polygons,
			Show[
				Graphics3D[Join[
					{EdgeForm[], Opacity[0.5]},
					Map[Polygon, polygons]], Boxed->False],
				rotationPlot[rots, pointStyles, {}, Automatic]]]];

Experimental I

Here we check that Wood (1993)’s matrix Fisher normalizing constant on SO(3) is correct. It seems to have an extraneous Pi^2.

Clear[a, p, t, p1, p2, p3];
p1 = 9;
p2 = 5;
p3 = 4;
f = Tr[DiagonalMatrix[{p1, p2, p3}].rotationFromAngleAxis[{a, cartesianFromSpherical[{1, p, t}]}]]
g = Sin[p] (1 - Cos[a]) / (4 Pi^2);
NIntegrate[Exp[f] g, {a, 0, Pi}, {p, 0, Pi}, {t, 0, 2 Pi}]
rotationFisherNormalizer[p1, p2, p3]

Clear[gs];
gs[rots_] :=
	Module[{g},
		g = rotationSignedSVD[Mean[rots]][[2]];
		{g[[1, 1]], g[[2, 2]], g[[3, 3]]}];
gs[rots]

Clear[optPhis];
optPhis[g1_?NumericQ, g2_?NumericQ, g3_?NumericQ] :=
	Module[{phi1, phi2, phi3, f, sol},
		f[phi1_?NumericQ, phi2_?NumericQ, phi3_?NumericQ] :=
			phi1 g1 + phi2 g2 + phi3 g3  - Log[rotationFisherNormalizer[phi1, phi2, phi3]];
		sol = NMaximize[f[phi1, phi2, phi3], {phi1, phi2, phi3}];
		{phi1, phi2, phi3} /. sol[[2]]];
Timing[optPhis[0.9934600947599254`,0.9902656323143567`,0.9853273716174588`]]
Timing[optPhis[0.3, 0.2, 0.1]]

phis = {40, 6, 5};
rots = rotationFisher[IdentityMatrix[3], DiagonalMatrix[phis], 1000];
gss = gs[rots]
Apply[optPhis, gss]

rots = rotationUniform[1000];
gss = gs[rots]
Apply[optPhis, gss]

optPhis[0.99, 0.5, 0.090]
optPhis[0.98, 0.5, 0.090]
optPhis[0.97, 0.5, 0.090]
optPhis[0.96, 0.5, 0.090]
optPhis[0.95, 0.5, 0.090]

delta = 0.055;
SetDirectory[NotebookDirectory[]];
DateString[]
posTable = Table[optPhis[g1, g2, g3],
	{g1, 0, 0.99, delta}, {g2, 0, g1, delta}, {g3, 0, g2, delta}];
Put[CForm[posTable], "posTable.R"]
DateString[]
negTable = Table[optPhis[-g1, g2, g3],
	{g1, 0, 0.99, delta}, {g2, 0, g1, delta}, {g3, 0, g2, delta}];
Put[CForm[negTable], "negTable.R"]
DateString[]

posTable[[1]]
posTable[[2]]
posTable[[3]]
posTable[[4]]

negTable[[1]]
negTable[[2]]
negTable[[3]]
negTable[[4]]

Experimental II

This Rancourt et al. (2000) plot is great, except that it doesn’t clip the confidence region against the boundary of the plot on those rare occasions when it has to.

Clear[trianglesSphere, trianglesEllipsoid, rotationRancourtPlot];
trianglesSphere[m_, n_] :=
	Module[{mesh, i, j},
		mesh = {};
		For[j = 0, j < n, j += 1,
			For[i = 0, i < m - 1, i += 1,
				mesh = Append[mesh,
					{cartesianFromSpherical[{1, i Pi / m, j 2 Pi / n}],
					cartesianFromSpherical[{1, (i + 1) Pi / m, j 2 Pi / n}],
					cartesianFromSpherical[{1, (i + 1) Pi / m, (j + 1) 2 Pi / n}]}]];
			For[i = 1, i < m, i += 1,
				mesh = Append[mesh,
					{cartesianFromSpherical[{1, i Pi / m, j 2 Pi / n}],
					cartesianFromSpherical[{1, (i + 1) Pi / m, (j + 1) 2 Pi / n}],
					cartesianFromSpherical[{1, i Pi / m, (j + 1) 2 Pi / n}]}]]];
		mesh];
trianglesEllipsoid[m_, n_, ell_] :=
	Module[{q, vals, d, ts},
		{vals, q} = Eigensystem[ell];
		q = Transpose[Map[unit, q]];
		d = DiagonalMatrix[vals^(-1 / 2)];
		ts = trianglesSphere[m, n];
		ts = Map[Function[p, q.d.p], ts, {2}];
		ts];
rotationRancourtPlot[rots_, alpha_, m_] :=
	Module[{mHat, g, f, tris},
		{mHat, g, f} = rotationRancourtMeanRegion[rots, alpha];
		tris = trianglesEllipsoid[m, 2 m, g / f];
		tris = Map[Function[v,
			volumetricFromRotation[rotationRancourtRotation[v, mHat]]], tris, {2}];
		Show[
			rotationPlot[rots, Automatic, {}, Automatic],
			Graphics3D[Join[{EdgeForm[], Opacity[0.5]}, Map[Polygon, tris]], Boxed->False]]];

data = rotationFisher[rotationFromVolumetric[{0, 0, 0.615}], DiagonalMatrix[{9, 8, 7}], 50];
rotationRancourtPlot[data, 0.05, 8]

Clear[volumetricLevelPlotC];
volumetricLevelPlotC[f_, contours_, minDepth_, maxDepth_] :=
	Module[{stream, r, fileToC, fileFromC, error, polygons},
		(* Generate the boilerplate preamble. *)
		!!
		(* Insert the custom data. *)
		!!
		(* Generate the rest of the boilerplate. *)
		!!
		(* Compile and run the C program. *)
		fileToC = Close[stream];
		stream = OpenWrite[];
		fileFromC = Close[stream];
		SetDirectory[NotebookDirectory[]];
		!!error = Run["./georotationsKamb.o " <> fileToC <> " " <> fileFromC];
		polygons = Get[fileFromC];
		DeleteFile[fileToC];
		DeleteFile[fileFromC];
		(* Return error string or graphics object. *)
		If[StringQ[polygons],
			polygons,
			Graphics3D[Join[
				{EdgeForm[], Opacity[0.5]},
				Map[Polygon, polygons]], Boxed->False]];

Geodesic methods modulo a subgroup

rotationDistance[r_, q_, g_] :=
	Min[Map[Function[a, rotationDistance[a.r, q]], g]];

rotationDiameter[rots_, g_] :=
	Module[{diameter, i, j},
		diameter = 0;
		For[i = 1, i ≤ Length[rots], i += 1,
			For[j = i + 1, j ≤ Length[rots], j += 1,
				diameter = Max[diameter, rotationDistance[rots[[i]], rots[[j]], g]]]];
		diameter];

rotationVariance[rots_, center_, g_] :=
	Mean[Table[rotationDistance[center, rots[[i]], g]^2, {i, 1, Length[rots]}]] / 2;

rotationClusters[rots_, gamma0_, epsilon_, bound_, g_] :=
	Module[{dd, t, gamma, wgd0, wgd1, error},
		(* Build a table of distances squared between the rotations *)
		dd = Table[rotationDistance[rots[[i]], rots[[j]], g]^2, {i, 1, Length[rots]}, {j, 1, Length[rots]}];
		(* Initialize the loop. *)
		gamma = gamma0;
		wgd1 = rotationClustersWGD[dd, gamma];
		error = epsilon + 1;
		t = 0;
		While[t < bound && error ≥ epsilon,
			(* Prepare for this iteration. *)
			t += 1;
			wgd0 = wgd1;
			(* Update the clustering. *)
			gamma = rotationClustersGamma[dd, gamma];
			(* Test the loop termination condition. *)
			wgd1 = rotationClustersWGD[dd, gamma];
			error = Abs[wgd1 - wgd0]];
		{gamma, error, t}];

rotationRegressionCurve[xrs_, g_] :=
	Module[{n, q, e, m32, m13, m21, w32, w13, w21, min, w, m},
		(* Let q be the r whose x is closest to zero. *)
		q = xrs[[argMin[Abs[xrs[[All, 1]]]][[1]], 2]];
		n = Length[xrs];
		e[m32_, m13_, m21_, w32_, w13_, w21_] :=
			Module[{mm, b, i},
				b = rotationExp[{{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}}].q;
				mm = {{0, -m21, m13}, {m21, 0, -m32}, {-m13, m32, 0}};
				Sum[rotationDistance[xrs[[i, 2]], rotationExp[xrs[[i, 1]] mm].b, g]^2, {i, 1, n}] / (2 n)];
		min = Minimize[
			{e[m32, m13, m21, w32, w13, w21],
			m32^2 + m13^2 + m21^2 ≤ Pi^2,
			w32^2 + w13^2 + w21^2 ≤ Pi^2},
			{m32, m13, m21, w32, w13, w21}];
		w = {{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}} /. min[[2]];
		m = {{0, -m21, m13}, {m21, 0, -m32}, {-m13, m32, 0}} /. min[[2]];
		{m, rotationExp[w].q, min[[1]]}];

rotationRegressionSurface[xyrs_, g_] :=
	Module[{norms, n, q, e, m32, m13, m21, l32, l13, l21, w32, w13, w21, min, w, m, l},
		(* Let q be the r whose (x, y) is closest to (0, 0). *)
		norms = Map[Norm, xyrs[[All, 1;;2]]];
		q = xyrs[[argMin[norms][[1]], 3]];
		n = Length[xyrs];
		e[m32_, m13_, m21_, l32_, l13_, l21_, w32_, w13_, w21_] :=
			Module[{mm, ll, b, i},
				b = rotationExp[{{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}}].q;
				mm = {{0, -m21, m13}, {m21, 0, -m32}, {-m13, m32, 0}};
				ll = {{0, -l21, l13}, {l21, 0, -l32}, {-l13, l32, 0}};
				Sum[rotationDistance[
						xyrs[[i, 3]],
						rotationExp[xyrs[[i, 1]] mm + xyrs[[i, 2]] ll].b,
						g]^2,
					{i, 1, n}] / (2 n)];
		min = Minimize[
			{e[m32, m13, m21, l32, l13, l21, w32, w13, w21],
			m32^2 + m13^2 + m21^2 ≤ Pi^2,
			l32^2 + l13^2 + l21^2 ≤ Pi^2,
			w32^2 + w13^2 + w21^2 ≤ Pi^2},
			{m32, m13, m21, l32, l13, l21, w32, w13, w21}];
		w = {{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}} /. min[[2]];
		m = {{0, -m21, m13}, {m21, 0, -m32}, {-m13, m32, 0}} /. min[[2]];
		l = {{0, -l21, l13}, {l21, 0, -l32}, {-l13, l32, 0}} /. min[[2]];
		{m, l, rotationExp[w].q, min[[1]]}];

rotationRegressionBandwidth[xrs_, k_, g_] :=
	Module[{f, h, sol},
		f[h_?NumericQ] :=
			Module[{rots, val},
				rots = Table[rotationRegressionKernel[
					Join[xrs[[;;(i - 1)]], xrs[[(i + 1);;]]],
					{xrs[[i, 1]]}, k, h][[1]], {i, 1, Length[xrs]}];
				val = Sum[rotationDistance[xrs[[i, 2]], rots[[i]], g]^2, {i, 1, Length[xrs]}];
				val];
		sol = NMinimize[{f[h], h > 0}, {h}];
		h /. sol[[2]]];
rotationRegressionBandwidthVerbose[xrs_, k_, g_] :=
	Module[{best, count, f, h, sol},
		best = 10^8;
		count = 0;
		f[h_?NumericQ] :=
			Module[{rots, val},
				count += 1;
				rots = Table[rotationRegressionKernel[
					Join[xrs[[;;(i - 1)]], xrs[[(i + 1);;]]],
					{xrs[[i, 1]]}, k, h][[1]], {i, 1, Length[xrs]}];
				val = Sum[rotationDistance[xrs[[i, 2]], rots[[i]], g]^2, {i, 1, Length[xrs]}];
				If[val < best,
					best = val;
					Print["best = ", best, ", h = ", h, ", count = ", count, ", time = ", DateString[]]];
				val];
		sol = NMinimize[{f[h], h > 0}, {h}];
		h /. sol[[2]]];
rotationRegressionKernel[xrs_, xs_, k_, h_, g_] :=
	Module[{n, sols, w32, w13, w21},
		n = Length[xrs];
		sols = Table[Minimize[{
			Sum[k[(xs[[j]] - xrs[[i, 1]]) / h] / h rotationDistance[
					xrs[[i, 2]],
					rotationExp[{{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}}],
					g],
				{i, 1, n}] / Sum[k[(xs[[j]] - xrs[[i, 1]]) / h] / h, {i, 1, n}],
			w32^2 + w13^2 + w21^2 ≤ Pi}, {w32, w13, w21}],
			{j, 1, Length[xs]}];
		Table[
			rotationExp[{{0, -w21, w13}, {w21, 0, -w32}, {-w13, w32, 0}} /. sols[[i, 2]]],
			{i, 1, Length[sols]}]];

rotationMean[rots_, epsilon_, t_, seed_, g_] :=
	Module[{mean, errorSquared, k, i, j, ws, ds, w},
		mean = seed;
		k = 0;
		errorSquared = epsilon^2 + 1;
		While[errorSquared ≥ epsilon^2 && k < t,
			w = {{0, 0, 0} ,{0, 0, 0}, {0, 0, 0}};
			For[i = 1, i ≤ Length[rots], i += 1,
				ws = Table[
					rotationLog[Transpose[mean].g[[j]].rots[[i]]],
					{j, 1, Length[g]}];
				ds = Table[frobenius[ws[[j]], ws[[j]]], {j, 1, Length[ws]}];
				j = argMin[ds][[1]];
				w += ws[[j]]];
			w = w / Length[rots];
			mean = mean.rotationExp[w];
			errorSquared = frobenius[w, w];
			k += 1];
		{mean, Sqrt[errorSquared], k}];