MATH590 DIF: Information geometry

In this final module homework we combine concepts of differential geometry with information theory to investigate the utility of the new field of information geometry to data analysis. We’ll use the same model system of viscoelastic flow in a convergent channel. Recall that this simple physical system corresponds to a generic situation encountered across the sciences: a system subjected to background stochastic  forcing and experiences a change in external forcing (i.e., the channel wall geometry in this case).

Viscoelastic flow in a convergent channel

Consider the flow of polymers in a convergent channel -1<x<1. The x-velocity of the background flow in the channel is assumed to be known:

u[x_]=(2/π)ArcTan[20x]+2;
Plot[u[x],{x,-1,1},Axes->False,Frame->True,FrameLabel->{"x","u(x)"},GridLines->Automatic,PlotLabel->"Background flow velocity"]

Assuming a dilute aqueous solution, the polymers respond to changes in the liquid velocity by changing their conformation. This process is used, for instance, to stretch DNA molecules prior to identifying base pairs. The simplest model for the complicated internal conformation is a vector Q from one end of the molecule to the other end, the so-called dumbbell model. 

We are interested in extracting the probability distribution functions for Q from data. In practice, this would be obtained by measurements. Here, however, we generate synthetic data from a stochastic process that models the stretching of the polymer by shear within the flow according to the stochastic differential equation

dQ = (η u’(x) Q - k Q) dt + σ dW 

with η,k,σ, model constants that correspond to viscosity, dumbbell stiffness, temperature, respectively. The dW term is a Wiener process modeling the Brownian motion of the dumbbell ends. Note that Q=0 is taken to signify that the “spring” between the two ends of the dumbbell is at its equilibrium, zero-force length

Thermal effects

Taking η=0 leads to a model that is only influenced by temperature dQ = - k Q dt +  σ dW, in which thermal effects force the dumbbell away from its Q=0 equilibrium length (expressed through the mean value μ=0). This is known as an Ornstein-Uhlenbeck process. Generate the synthetic data for 5 instances (i.e., 5 polymers) over 100 time steps of 0.01s.

η=0; k=1; σ=10; μ=0; x=0; nPoly=5;
proc=OrnsteinUhlenbeckProcess[μ,σ,η u'[x]+ k];
db=Table[RandomFunction[proc,{0,1,0.01},2]["ValueList"]\[Transpose],{nPoly}];

Represent the time-evolution of the two components of Q=(Q1,Q2)

ListPlot[db,Joined->True,AspectRatio->Automatic,Frame->True,Axes->False,FrameLabel->{"Q1","Q2"}]

Also construct a direct representation of the evolution of the dumbbell vector itself

For[i=1; dbVec={}, i≤nPoly, i++,AppendTo[dbVec, Graphics[{ColorData[i,"ColorList"][[i]],Map[Line[{{0,0},#}]& ,db[[i]]]}]] ];

Show[dbVec,Frame->True,GridLines->Automatic]

The above can be interpreted as depicting the evolution over 100 time increments of the end-to-end (dumbbell) vector for nPoly=5 different polymers.

Flow effects

Now turn on viscosity, by setting η>0

η=.1; k=1.; σ=10.; μ=0; x=0; nPoly=5;
proc=OrnsteinUhlenbeckProcess[μ,σ,η u'[x]+ k];
db=Table[RandomFunction[proc,{0,1,0.01},2]["ValueList"]\[Transpose],{nPoly}];

Represent the time-evolution of the two components of Q=(Q1,Q2)

ListPlot[db,Joined->True,AspectRatio->Automatic,Frame->True,Axes->False,FrameLabel->{"Q1","Q2"}]

Karhunen-Loeve, SVD

In the previous section, a simplified physical model was used to generate synthetic data as opposed to measured data (e.g., accelerometer measurements, images). Consider now the main problem of interest in this module of extracting an approximation of the PDF of the data. The PDF itself is unknown beforehand, but can be approximated as a multivariate Gaussian distribution through the singular value decomposition (SVD). The goal in this module is to build a point cloud of such multivariate Gaussian approximations. In subsequent modules, the point cloud is used to construct non-Gaussian PDF approximants.

Generate data

η=.1; k=1.; σ=10.; μ=0; x=0; nPoly=50;
proc=OrnsteinUhlenbeckProcess[μ,σ,η u'[x]+ k];
data=Flatten[Table[RandomFunction[proc,{0,1,0.01},2]["ValueList"]\[Transpose],{nPoly}],1];

m=Mean[data];
X=Map[Plus[#,-m]&,data];
{Ns,Nd}=Dimensions[X]

Compute the SVD, limited to two modes. The singular values specify the data variance

{U,S,V}=SingularValueDecomposition[X,2];
MatrixForm[N[(S^2)/Ns]]

Geometry of the 2-sphere

Metric tensor

Relate traditional notation for parameters to indexed notation

ξ={ϕ,θ}

g={{1,0},{0,((Sin[ϕ])^2)}}

ginv=Inverse[g]

Connection coefficients from metric tensor (Christoffel symbols)

Γ is a rank 3 tensor. First, we compute (Subscript[((Γ^1)), ij]), (Subscript[((Γ^2)), ij]) separately

Γ1=Table[  Sum[ (1/2) ginv[[1,l]] (D[g[[l ,i]],ξ[[ j]]] + D[g[[ l,j]],ξ[[i ]]] - D[g[[i ,j]],ξ[[l ]]]) ,{l,2}] ,{i,2},{j,2}]

MatrixForm[Γ1]

Γ2=Table[  Sum[ (1/2) ginv[[2,l]] (D[g[[l ,i]],ξ[[ j]]] + D[g[[ l,j]],ξ[[i ]]] - D[g[[i ,j]],ξ[[l ]]]) ,{l,2}] ,{i,2},{j,2}]

Now, introduce a general formulation for computation of the full rank-3 tensor

Γ=Table[  Sum[ (1/2) ginv[[k,l]] (D[g[[l ,i]],ξ[[ j]]] + D[g[[ l,j]],ξ[[i ]]] - D[g[[i ,j]],ξ[[l ]]]) ,{l,2}] ,{k,2},{i,2},{j,2}]

Compare the two computational methods to check they give the same result

{MatrixForm[Γ[[1]]],MatrixForm[Γ1]}

{MatrixForm[Γ[[2]]],MatrixForm[Γ2]}

Also define connection coefficient functions along a curve

Γt = Γ /. {ϕ->ϕ[t],θ->θ[t]}

Geodesics

Define the system of geodesic equations

GeodesicSystem =Table[D[ξ[[k]][t],{t,2}] + Sum[Γt[[k,i,j]] D[ξ[[i]][t],t]  D[ξ[[j]][t],t],{i,2},{j,2}]==0 ,{k,2}]

The above is a system of two ordinary differential equations of second order. Four scalar initial or boundary value conditions are required. First consider four initial conditions

InitialConditions = {ϕ[0]==(π/2),ϕ'[0]==-1,θ[0]==0,θ'[0]==1}

First, try to find an analytical solution. Note that finding an analytical solution to this nonlinear system is not guaranteed.

DSolve[ Flatten[{GeodesicSystem,InitialConditions}],{ϕ[t],θ[t]},t]

Next try to find a numerical solution

sol=NDSolve[ Flatten[{GeodesicSystem,InitialConditions}],{ϕ[t],θ[t]},{t,0,1}][[1]]

ParametricPlot3D[{Sin[ϕ[t]] Cos[θ[t]],Sin[ϕ[t]] Sin[θ[t]],Cos[ϕ[t]]} /. sol,{t,0,1}]

GeodesicRays[system_, nG_,tf_]:=Module[{α,geodesics,sol,g,InitialConditions},
α=Table[2π(k-1)/nG,{k,nG}];
InitialConditions = Table[
{ϕ[0]==(π/2),ϕ'[0]==Cos[ α[[k]]],θ[0]==0,θ'[0]==Sin[α[[k]]]},{k,nG}];
geodesics={};
For[k=1,k≤nG,k++,
sol=NDSolve[ Flatten[{system,InitialConditions[[k]]}],{ϕ[t],θ[t]},{t,0,tf}][[1]];
g=ParametricPlot3D[{Sin[ϕ[t]] Cos[θ[t]],Sin[ϕ[t]] Sin[θ[t]],Cos[ϕ[t]]} /. sol,{t,0,tf}];
AppendTo[geodesics,g];
];
Return[Show[geodesics]]
]

GeodesicRays[GeodesicSystem,18,π]

Homework questions

Question 1

Compute the Fisher metric for univariate Gaussian distributions.

Question 2

Formulate and solve the geodesic equation for univariate Gaussian distributions. Plot the geodesics in the (μ,σ) plane, the Poincare hyperbolic plane.

Question 3

Plot the trajectory of the Gaussian PDFs for Qx,Qy obtained by the SVD in the Poincare hyperbolic plane along the channel flow for viscosities η=10^(-k), k=1,2,3,4,5.

Question 4

Arbitrarily choose one of the viscosities. Recompute the Qx,Qy PDFs for a larger number of sample paths, e.g., nPoly=250 instead of nPoly=50.  Plot the trajectory in the Poincare plane of the PDFs for both the nPoly=250 and the nPoly=50 cases.

Question 5

Presumably the means of Qx, Qy obtained for the larger number of samples (nPoly=250) is more accurate. Correct the trajectory for the nPoly=50 case by geodesic transport onto the means obtained from the nPoly=250 case. Plot the initial and corrected trajectories in the Poincare plane.