Framing
An animated version of a graphic I produced for my paper “Stiefel manifolds and polygons”, which appears in the Proceedings of Bridges 2019: Mathematics, Art, Music, Architecture, Education, Culture. The actual conference, which is an annual math art conference, happens in two weeks in Linz, Austria; if anybody is planning to attend, I will be giving a talk on Thursday, and also have two pieces (A to Z and Truncation) in the art exhibit.
To explain the animation a bit more: as I’ve mentioned before , there’s a correpondence between polygons in space and points in the Grassmann manifold G_2(\mathbb{C}^n) of 2-dimensional linear subspaces of \mathbb{C}^n (by polygon in space I mean a closed, piecewise-linear curve in \mathbb{R}^3).
See that post for more details, but what I didn’t mention there is that the correspondence is really between points in the Grassmannian and framed polygons, meaning that at each edge I have a triple of mutually perpendicular vectors: one pointing in the direction of the edge and the other two giving a perpendicular basis for the orthogonal complement of the edge direction (i.e., the normal space).
While we don’t require these triples to be othonormal, we do require them to be mutually orthogonal and the same length, so we can think of the matrix whose columns are the vectors in the triple as a scaled element of SO(3), the group of rotations of \mathbb{R}^3. Then the process of going from a framed polygon to the Grassmannian is the following: convert each scaled element of SO(3) to a quaternion (formally, this means lifting to a scaled element of \mathrm{Spin}(3)), interpret the quaternion as u + v j, where u and v are complex numbers, let \vec{u} and \vec{v} be the complex n-vectors of the u’s and v’s, respectively, and then take the span of \vec{u} and \vec{v}.
The virtue of doing this is that geodesics in the Grassmannian are easy to compute explicitly, which gives us a nice way of morphing between different framed polygons. The animation shows the geodesic in the Grassmannian between the +3-framed regular 200-gon (meaning that the frame vectors spin three times around the edge vectors as we travel all the way around the polygon), and the 0-framed 200-gon (in this case the frame vectors just always point straight up and straight out). The curves on the surface are all parallel to the paths traced out by the frame vectors.
It’s probably not a coincidence that there are 3 singularities and we’re going from a +3-framed polygon to a 0-framed polygon: getting rid of twist in the framing is expensive, and it seems to optimal to scale some edges down near length 0, then twist, then scale back up, so one needs to do this three times to get rid of three twists in the framing.
There’s a fair amount of code that goes into this. First of all, we need to get the Quaternions package:
<< Quaternions`
Also, a couple of helper functions for converting between real and complex matrices:
Real4nToComplex2n[M_] :=
Transpose[Complex @@ # & /@ # & /@ (Partition[#, 2] & /@ Transpose[M])];
Complex2nToReal4n[M_] := Transpose[Flatten /@ (ReIm /@ Transpose[M])];
For lifting from a scaled element of SO(3) to a quaternion, I’m using a slight modification of the function rotm2q from Reichlin’s quaternion package for Octave:
MatrixToQuaternion[A_] := Module[{M, T, l, s, w, x, y, z},
l = Norm[A[[1]]];
M = 1/l A;
T = Tr[M];
Which[T > 0,
s = 1/(2 Sqrt[T + 1]);
w = 1/(4 s);
x = s (M[[2, 3]] - M[[3, 2]]);
y = s (M[[3, 1]] - M[[1, 3]]);
z = s (M[[1, 2]] - M[[2, 1]]);,
M[[1, 1]] > M[[2, 2]] && M[[1, 1]] > M[[3, 3]],
s = 2 Sqrt[1 + M[[1, 1]] - M[[2, 2]] - M[[3, 3]]];
w = (M[[2, 3]] - M[[3, 2]])/s;
x = s/4;
y = (M[[1, 2]] + M[[2, 1]])/s;
z = (M[[1, 3]] + M[[3, 1]])/s;,
M[[2, 2]] > M[[3, 3]],
s = 2 Sqrt[1 + M[[2, 2]] - M[[1, 1]] - M[[3, 3]]];
w = (M[[3, 1]] - M[[1, 3]])/s;
x = (M[[1, 2]] + M[[2, 1]])/s;
y = s/4;
z = (M[[2, 3]] + M[[3, 2]])/s;,
True,
s = 2 Sqrt[1 + M[[3, 3]] - M[[1, 1]] - M[[2, 2]]];
w = (M[[1, 2]] - M[[2, 1]])/s;
x = (M[[1, 3]] + M[[3, 1]])/s;
y = (M[[2, 3]] + M[[3, 2]])/s;
z = s/4;
];
Sqrt[l] {w, x, y, z}
];
On the other hand, to go from a point in the Grassmannian (or Stiefel manifold) to a framed polygon, we use these functions:
QuaternionFramedEdges[
StiefelForm_] := {Conjugate[#] ** Quaternion[0, 1, 0, 0] ** #,
Conjugate[#] ** Quaternion[0, 0, 1, 0] ** #,
Conjugate[#] **
Quaternion[0, 0, 0, 1] ** #} & /@ (Quaternion @@ # & /@
Transpose[StiefelForm]);
QFramedEdges[StiefelForm_] :=
List @@ #[[2 ;;]] & /@ # & /@ QuaternionFramedEdges[StiefelForm];
Next, geodesics in the Grassmannian are computed very straightforwardly: suppose P_1 and P_2 are our two 2-dimensional linear subspaces in \mathbb{C}^n. Then the idea is to find the orthonormal bases for P_1 and P_2 which are as close to each other as possible, and then simply directly rotate one to the other. More precisely, we want to find unit vectors a_i \in P_i so that |\langle a_1, a_2\rangle | is maximal, then find unit b_i \in a_i^\bot \cap P_i that maximize |\langle b_1, b_2\rangle |. Then (a_i,b_i) is an orthonormal basis for P_i, and we realize the geodesic from P_1 to P_2 while directly rotating a_1 to a_2 and simultaneously rotating b_1 to b_2 (with speeds chosen so that both rotation finish at the same time). In turn, an easy way to find a_1 and b_1 is that they are the leading eigenvectors of the linear map which orthogonally projects from P_1 to P_2, and then back from P_2 to P_1.
So here’s some code for computing geodesics:
ComplexProjectionBasis[{A_, B_}, {C_, D_}] :=
Eigenvectors[
Transpose[{A, B}].Conjugate[{A, B}].Transpose[{C, D}].Conjugate[{C, D}], 2];
ComplexGrassmannianGeodesic[{A_, B_}, {C_, D_}, t_] :=
Module[{a, b, c, d, cPerp, dPerp, dist1, dist2},
{a, b} = Orthogonalize[ComplexProjectionBasis[{C, D}, {A, B}]];
{c, d} = Orthogonalize[ComplexProjectionBasis[{A, B}, {C, D}]];
{cPerp, dPerp} = {Normalize[c - Re[Conjugate[c].a]*a],
Normalize[d - Re[Conjugate[d].b]*b]};
dist1 = ArcCos[Re[Conjugate[c].a]];
dist2 = ArcCos[Re[Conjugate[d].b]];
{Cos[t*dist1]*a + Sin[t*dist1]*cPerp,
Cos[t*dist2]*b + Sin[t*dist2]*dPerp}
];
Finally, then, we just need to define our two framed 200-gons:
circleframe = With[{n = 200},
1/10 Real4nToComplex2n[
Transpose[
MatrixToQuaternion /@ (Table[
Transpose@
Prepend[{{-#[[i, 2]], #[[i, 1]], 0}, {0, 0, 1}}, #[[i]]],
{i, 1, n}]
&[Append[#, 0] & /@ N[CirclePoints[n]]])]]
];
circleframe2 = With[{n = 200}, 1/10 Real4nToComplex2n[
Transpose[
MatrixToQuaternion /@ (Table[
Transpose@
Prepend[
RotationMatrix[
6 ? i/n].{{-#[[i, 2]], #[[i, 1]], 0}, {0, 0, 1}}, #[[i]]],
{i, 1, n}]
&[Append[#, 0] & /@ N[CirclePoints[n]]])]]
];
And put everything together:
DynamicModule[{frame, n, f, g, h,
cols = RGBColor /@ {"#f6490d", "#1dced8", "#161c2e"}},
Manipulate[
frame =
Transpose[
QFramedEdges[
Complex2nToReal4n[
ComplexGrassmannianGeodesic[circleframe, circleframe2, t]]]];
f = Interpolation[
Transpose[{N@Range[0, 1, 1/(Length[#] - 1)], #}]] &[
Prepend[#, #[[-1]]] &[# - ConstantArray[Mean[#], Length[#]] &[
Accumulate[frame[[1]]]]]];
g = Interpolation[
Transpose[{N@Range[0, 1, 1/(Length[#] - 1)], #}]] &[
Prepend[#, #[[-1]]] &[frame[[2]]]];
h = Interpolation[
Transpose[{N@Range[0, 1, 1/(Length[#] - 1)], #}]] &[
Prepend[#, #[[-1]]] &[frame[[3]]]];
n = Length[frame[[1]]];
ParametricPlot3D[
f[t] + 3 (Cos[?] g[t] + Sin[?] h[t]), {t, 0, 1}, {?, 0, 2 ?},
MeshFunctions -> {#5 &}, Mesh -> 5,
MeshShading -> {cols[[1]], cols[[2]]}, MeshStyle -> None,
PlotPoints -> 50, Boxed -> False,
ViewPoint -> -100 Normalize[Cross[f[0.], f[.25]]],
ViewVertical -> f[0.], ViewAngle -> ?/300, PlotRange -> 1/2,
ImageSize -> {540, 540}, Axes -> False, Lighting -> "Neutral",
Background -> cols[[-1]]],
{t, 0., 1}]
]
