Geodesics and the Matrix Exponential.

Jason Cantarella, January 2013.

This demonstration shows some code for the matrix exponential and for geodesics in some particular Grassmann and Stiefel manifolds. First, we note that Mathematica has a matrix exponential function built in:

A = {{0, -0.25},{0.25,0}};
A // MatrixForm

MatrixExp[A]  // MatrixForm

We can check immediately that the matrix resulting from a skew-symmetric input is orthogonal:

MatrixExp[A].Transpose[MatrixExp[A]] // MatrixForm

We now want to define entries in some Grassmann and Stiefel manifolds and compute geodesics from them. We first note that we can get a random entry in O(n) or U(n) by taking the QR decomposition of a matrix with Gaussian entries.

RandomOrthogonal[n_] := QRDecomposition[RandomVariate[NormalDistribution[],{n,n}]][[1]];

A = RandomOrthogonal[3];
A // MatrixForm

A.Transpose[A] // MatrixForm

ToComplex[{a_,b_}] := a + b I;
RandomUnitary[n_] := QRDecomposition[Map[ToComplex,RandomVariate[BinormalDistribution[0],{n,n}],{2}]][[1]];

A = RandomUnitary[3];
A // MatrixForm

A.ConjugateTranspose[A] // MatrixForm

Now we’d like some code to generate “horizontal” directions in the Stiefel and Grassmannian sense. Notice that Mathematica has a syntax for address blocks of a matrix directly:

A[[i;;j,k;;l]]

is the (j-i) x (l-k) submatrix of A. Using this, we can write down elements in our horizontal spaces just by making large skew-symmetric matrices and then zeroing certain blocks.

RandomStiefelHorizontal[n_,k_] := Module[{A,i,j},
A = RandomVariate[NormalDistribution[],{n,n}];
A = A - Transpose[A];
A[[k+1;;n, k+1;;n]] = Table[0,{i,1,n-k},{j,1,n-k}];
A
];

RandomStiefelHorizontal[5,3] // MatrixForm

RandomGrassmannHorizontal[n_,k_] := Module[{A,i,j},
A = RandomVariate[NormalDistribution[],{n,n}];
A = A - Transpose[A];
A[[k+1;;n, k+1;;n]] = Table[0,{i,1,n-k},{j,1,n-k}];
A[[1;;k,1;;k]] = Table[0,{i,1,k},{j,1,k}];
A
];

RandomGrassmannHorizontal[5,3] // MatrixForm

Doing the same thing for complex matrices just means we have to generate skew-Hermitian matrices first:

RandomComplexStiefel[n_,k_] := Module[{A,i,j},
A = Map[ToComplex,RandomVariate[BinormalDistribution[0],{n,n}],{2}];
A = A - ConjugateTranspose[A];
A[[k+1;;n, k+1;;n]] = Table[0,{i,1,n-k},{j,1,n-k}];
A
];
RandomComplexGrassmannHorizontal[n_,k_] := Module[{A,i,j},
A = Map[ToComplex,RandomVariate[BinormalDistribution[0],{n,n}],{2}];
A = A - ConjugateTranspose[A];
A[[k+1;;n, k+1;;n]] = Table[0,{i,1,n-k},{j,1,n-k}];
A[[1;;k,1;;k]] = Table[0,{i,1,k},{j,1,k}];
A
];

Finally, we’ve got orthogonal and unitary matrices, and all kinds of tangent information ready to go. So we’ll start with the Stiefel manifold V_2(R^3).

StiefelVis[A_] := Module[{x,y},
x = A[[All,1]];
y = A[[All,2]];
Graphics3D[{Arrowheads[{.03}],Thickness[0.005],Arrow[{{0,0,0},x}],Arrow[{{0,0,0},y}],
{FaceForm[Opacity[0.2]],Polygon[{-x - y,x - y,x + y,-x + y}]}},Boxed->False]
];

MyAxes = Graphics3D[{Arrowheads[{-.03,.03}],Arrow[{{-1,0,0},{1,0,0}}],Arrow[{{0,-1,0},{0,1,0}}],Arrow[{{0,0,-1},{0,0,1}}]}];

A = RandomOrthogonal[3];

Show[StiefelVis[A],MyAxes]

We now need to go ahead and compute the geodesics. We know that the general formula is given by 

Q(t) = Q e^(t A),

where A is the (skew-symmetric or skew-hermitian) tangent vector to O(n) or U(n) so that the tangent vector to the geodesic at Q(0) is Q(0) A.

Geodesic[Base_,Dir_,t_] := Base. MatrixExp[t Dir];

We now try some demos.

1. Geodesic in a high-dimensional Stiefel manifold.

We apply the matrix exponential to see how the geodesic formula works in practice in the Stiefel manifold V_7(R^23) of 7-frames in R^23.

Base = RandomOrthogonal[23];
MatrixPlot[Base]

MatrixPlot[Base.Transpose[Base]]

Dir = RandomStiefelHorizontal[23,7];
MatrixPlot[Dir]

PlotLabel->Style["Endpoint Q(t) of Geodesic",Bold,FontFamily->"Helvetica",FontSize->14]

This plot shows the geodesic from Base in Direction “dir” in the orthogonal group (which lies over the corresponding geodesic in the Stiefel manifold) as a pair of MatrixPlots.

Manipulate[
MyA =Geodesic[Base,Dir,t];
GraphicsGrid[{{MatrixPlot[MyA,PlotLabel->Style["Endpoint Q(t) of Geodesic",Bold,FontFamily->"Helvetica",FontSize->14]],MatrixPlot[MyA.Transpose[MyA],PlotLabel->Style["Q(t).(Q(t))^T",Bold,FontFamily->"Helvetica",FontSize->14]]}},ImageSize->Large],{t,0,5}]

Of course, we have to be worried about the error in our computation: is the matrix “drifting” away from the orthogonal group? To check, we can plot the Frobenius norm distance between Q(t).(Q(t))^T and the identity matrix.

DistanceToI[Base_,Dir_,t_] := Module[{A},
A = Geodesic[Base,Dir,t];
Norm[A.Transpose[A] - IdentityMatrix[{Length[A],Length[A]}]]
];

DistanceToI[Base,Dir,0.2]

Plot[DistanceToI[Base,Dir,t],{t,0,5},ImageSize->Large,PlotLabel->Style["Distance from Q(t).Q(t)^T to I",Bold,FontFamily->"Helvetica",FontSize->14]]

Plot[DistanceToI[Base,Dir,t],{t,0,50},ImageSize->Large,PlotLabel->Style["Distance from Q(t).Q(t)^T to I",Bold,FontFamily->"Helvetica",FontSize->14]]

2. Geodesic in a low-dimensional Stiefel Manifold

Of course, the really nice picture will be of a geodesic in the manifold where we can see what’s going on: V_2(R^3). We try a random geodesic:

LowBase = RandomOrthogonal[3];
LowDir    = RandomStiefelHorizontal[3,2];

Manipulate[Show[{StiefelVis[Geodesic[LowBase,LowDir,t]],MyAxes},ImageSize->Large,PlotRange->{{-1.5,1.5},{-1.5,1.5},{-1.5,1.5}}],{t,0,10}]

3. Geodesic in a Grassmann manifold.

We now look at a geodesic in a Grassmann manifold.

GBase = RandomOrthogonal[12];
GDir = RandomGrassmannHorizontal[12,2];
MatrixPlot[GDir]

Manipulate[
GrA =Geodesic[GBase,GDir,t];
GraphicsGrid[{{MatrixPlot[GrA,PlotLabel->Style["Endpoint Q(t) of Geodesic",Bold,FontFamily->"Helvetica",FontSize->14]],MatrixPlot[GrA.Transpose[GrA],PlotLabel->Style["Q(t).(Q(t))^T",Bold,FontFamily->"Helvetica",FontSize->14]]}},ImageSize->Large],{t,0,5}]

4. Grassmannians and Polygons.

One of the applications of Grassmann and Stiefel manifolds that we’ll return to over and over again in the course will be the identification of Grassmannians and Stiefel manifolds with covers of closed polygon spaces of fixed total length. It’s worth writing down exactly how this is supposed to work.

MatrixToPlanePolygon[A_] := Accumulate[{#[[1]]^2 - #[[2]]^2,2 #[[1]] #[[2]]} & /@ A[[All,1;;2]]];

What we’re doing here is isolating the first two columns of the matrix, transposing it, and then viewing the ith column of the resulting 2 x n matrix as giving the real and imaginary components of a complex number z_i. If we take z_i^2, we get the complex number 

((z_i)_1^2 - (z_i)_2^2)  + (2 (z_i)_1 (z_i)_2) I

We note that the sum of these complex number is 0 iff the original pair of columns in A were orthogonal vectors of the same norm in R^n. 

Now the norm of z_i^2 is the square of the norm of z_i, which is exactly (z_i)_1^2 + (z_i)_2^2. Summing these over the coordinates yields the sum of the squares of the norms of the first two columns of A, which is always 2, since these are orthonormal.

DisplayPlanePoly[P_] := Graphics[{Arrowheads[0.02],Arrow /@ Transpose[{P,RotateRight[P]}]},PlotRange->{{-1,1},{-1,1}},Axes->True,AxesOrigin->{0,0},ImageSize->Large];

We can now take a look at some geodesics in the space of polygons:

PGBase = RandomOrthogonal[4];
PGDir = RandomGrassmannHorizontal[4,2];

Manipulate[DisplayPlanePoly[MatrixToPlanePolygon[Geodesic[PGBase,PGDir,t]]],{t,0,10}]

Questions: 

What is the distance on polygon space that is being minimized here? 
How can we understand the geodesic from the variation which generates it? 
Which geodesics are closed?