"OrthogonalProjection" Method for NDSolve

Introduction

Consider the matrix differential equation

where the initial value (Subscript[y, 0])==y(0)∈(ℝ^(m×p)) is given. Assume that ((Subscript[y, 0])^T)(Subscript[y, 0])==I, that the solution has the property of preserving orthonormality, ((y(t))^T)y(t)==I, and that it has full rank for all t≥0.

From a numerical perspective, a key issue is how to numerically integrate an orthogonal matrix differential system in such a way that the numerical solution remains orthogonal. There are several strategies that are possible. One approach, suggested in [DRV94], is to use an implicit Runge–Kutta method (such as the Gauss scheme). Some alternative strategies are described in [DV99] and [DL01].

The approach taken here is to use any reasonable numerical integration method and then postprocess using a projective procedure at the end of each integration step.

An important feature of this implementation is that the basic integration method can be any built-in numerical method, or even a user-defined procedure. In the following examples an explicit Runge–Kutta method is used for the basic time stepping. However, if greater accuracy is required an extrapolation method could easily be used, for example, by simply setting the appropriate (RefLink[Method, paclet:ref/Method]) option.

Projection Step

At the end of each numerical integration step you need to transform the approximate solution matrix of the differential system to obtain an orthogonal matrix. This can be carried out in several ways (see for example [DRV94] and [H97]):

The Newton and Schulz methods are quadratically convergent, and the number of iterations may vary depending on the error tolerances used in the numerical integration. One or two iterations are usually sufficient for convergence to the orthonormal polar factor (see the following) in IEEE double-precision arithmetic.

QR decomposition is cheaper than singular value decomposition (roughly by a factor of two), but it does not give the closest possible projection.

Definition (Thin singular value decomposition [GVL96]): Given a matrix A∈(ℝ^(m×p)) with m≥p there exist two matrices U∈(ℝ^(m×p)) and V∈(ℝ^(p×p)) such that (U^T) A V is the diagonal matrix of singular values of A, Σ==diag((Subscript[σ, 1]),…,(Subscript[σ, p]))∈(ℝ^(p×p)), where (Subscript[σ, 1])≥⋯≥(Subscript[σ, p])≥0. U has orthonormal columns and V is orthogonal.

Definition (Polar decomposition): Given a matrix A and its singular value decomposition U Σ (V^T), the polar decomposition of A is given by the product of two matrices Z and P where Z==U (V^T) and P==V Σ (V^T). Z has orthonormal columns and P is symmetric positive semidefinite.

The orthonormal polar factor Z of A is the matrix that solves

for the 2 and Frobenius norms [H96].

Schulz Iteration

The approach chosen is based on the Schulz iteration, which works directly for m≥p. In contrast, Newton iteration for m>p needs to be preceded by QR decomposition.

Comparison with direct computation based on the singular value decomposition is also given.

The Schulz iteration is given by

(Subscript[Y, i+1])==(Subscript[Y, i])+(Subscript[Y, i])(I-(Subscript[Y, i]^T) (Subscript[Y, i]))/2,(Subscript[Y, 0])==A.

The Schulz iteration has an arithmetic operation count per iteration of 2(m^2)p+2m (p^2) floating-point operations, but is rich in matrix multiplication [H97].

In a practical implementation, GEMM-based level 3 BLAS of LAPACK [LAPACK99] can be used in conjunction with architecture-specific optimizations via the Automatically Tuned Linear Algebra Software [ATLAS00]. Such considerations mean that the arithmetic operation count of the Schulz iteration is not necessarily an accurate reflection of the observed computational cost. A useful bound on the departure from orthonormality of A is in [H89]: ||(A^T) A-I(Subscript[||, F]). Comparison with the Schulz iteration gives the stopping criterion ||(A^T) A-I(Subscript[||, F])<τ for some tolerance τ.

Standard Formulation

Assume that an initial value (Subscript[y, n]) for the current solution of the ODE is given, together with a solution (Subscript[y, n+1])==(Subscript[y, n])+Δ (Subscript[y, n]) from a one-step numerical integration method. Assume that an absolute tolerance τ for controlling the Schulz iteration is also prescribed.

The following algorithm can be used for implementation.

Step 1. Set (Subscript[Y, 0])==(Subscript[y, n+1]) and i==0.

Step 2. Compute E==I-(Subscript[Y, i]^T) (Subscript[Y, i]).

Step 3. Compute (Subscript[Y, i+1])==(Subscript[Y, i])+(Subscript[Y, i])E/2.

Step 4. If ||E(Subscript[||, F])≤τ or i==(Subscript[i, max]), then return (Subscript[Y, i+1]).

Step 5. Set i==i+1 and go to step 2.

Increment Formulation

(RefLink[NDSolve, paclet:ref/NDSolve]) uses compensated summation to reduce the effect of rounding errors made by repeatedly adding the contribution of small quantities Δ (Subscript[y, n]) to (Subscript[y, n]) at each integration step [H96]. Therefore, the increment Δ (Subscript[y, n]) is returned by the base integrator.

An appropriate orthogonal correction Δ (Subscript[Y, i]) for the projective iteration can be determined using the following algorithm.

Step 1. Set Δ (Subscript[Y, 0])==0 and i==0.

Step 2. Set (Subscript[Y, i])==Δ (Subscript[Y, i])+(Subscript[y, n+1]).

Step 3. Compute E==I-(Subscript[Y, i]^T) (Subscript[Y, i]).

Step 4. Compute Δ (Subscript[Y, i+1])==Δ (Subscript[Y, i])+(Subscript[Y, i])E/2.

Step 5. If ||E(Subscript[||, F])≤τ or i==(Subscript[i, max]), then return Δ (Subscript[Y, i+1])+Δ (Subscript[y, n]).

Step 6. Set i==i+1 and go to step 2.

This modified algorithm is used in "OrthogonalProjection" and shows an advantage of using an iterative process over a direct process, since it is not obvious how an orthogonal correction can be derived for direct methods.

Examples

Orthogonal Error Measurement

A function to compute the Frobenius norm ||A(Subscript[||, F]) of a matrix A can be defined in terms of the (RefLink[Norm, paclet:ref/Norm]) function as follows:

FrobeniusNorm[a_?MatrixQ]:= Norm[a, Frobenius];

An upper bound on the departure from orthonormality of A can then be measured using this function [H97]:

OrthogonalError[a_?MatrixQ]:= FrobeniusNorm[Transpose[a].a - IdentityMatrix[Last[Dimensions[a]]]];

This defines the utility function for visualizing the orthogonal error during a numerical integration:

(* Utility function for extracting a list of values of the independent variable at which the integration method has sampled *)

TimeData[{v_?VectorQ,___?VectorQ}]:=TimeData[v];

TimeData[{if:(InterpolatingFunction[__])[_],___}]:=
Part[if,0,3,1];

(* Utility function for plotting the orthogonal error in a numerical integration *)

OrthogonalErrorPlot[sol_]:=
Module[{errdata, samples, soldata},
(* Form a list of times at which the method is invoked *)
samples = TimeData[sol];
(* Form a list of solutions at the integration times *)
soldata=Map[(sol/.t->#)&,samples];
(* Form a list of the orthogonal errors *)
errdata = Map[OrthogonalError, soldata];
ListLinePlot[Transpose[{samples,errdata}], Frame->True,PlotLabel->"Orthogonal error ||(Y^T)Y - I\!(Subscript[||, F]) vs time"
]
];

Square Systems

This example concerns the solution of a matrix differential system on the orthogonal group (Subscript[O, 3])(ℝ) (see [Z98]).

The matrix differential system is given by

with

and

The solution evolves as

The eigenvalues of Y(t) are (Subscript[λ, 1])==1, (Subscript[λ, 2])==exp(t i (Sqrt[3])), (Subscript[λ, 3])==exp(-ti(Sqrt[3])). Thus as t approaches π/(Sqrt[3]), two of the eigenvalues of Y(t) approach -1. The numerical integration is carried out on the interval [0,2]:

n = 3;

A = (({{0, -1, 1}, {1, 0, 1}, {-1, -1, 0}}));

Y=Table[y[i,j][t],{i,n},{j,n}];

F =A +( IdentityMatrix[n] -Transpose[Y].Y);

(* Vector differential system *)

system=Thread[Flatten[D[Y,t]]==Flatten[F.Y]];

(* Vector initial conditions *)

ics=Thread[Flatten[(Y/.t->0)]==Flatten[IdentityMatrix[Length[Y]]]];

eqs={system,ics};

vars = Flatten[Y];

time={t,0,2};

This computes the solution using an explicit Runge–Kutta method. The appropriate initial step size and method order are selected automatically, and the step size may vary throughout the integration interval, which is chosen in order to satisfy local relative and absolute error tolerances. Alternatively, the order of the method could be specified by using a (RefLink[Method, paclet:ref/Method]) option:

solerk=NDSolve[eqs,vars,time,Method->"ExplicitRungeKutta"];

This computes the orthogonal error, or absolute deviation from the orthogonal manifold, as the integration progresses. The error is of the order of the local accuracy of the numerical method:

solerk = Y/.First[solerk];

OrthogonalErrorPlot[solerk]

This computes the solution using an orthogonal projection method with an explicit Runge–Kutta method used for the basic integration step. The initial step size and method order are the same as earlier, but the step size sequence in the integration may differ:

solop=NDSolve[eqs,vars,time,Method->{"OrthogonalProjection",
Method->"ExplicitRungeKutta", Dimensions->Dimensions[Y]}];

Using the orthogonal projection method, the orthogonal error is reduced to approximately the level of roundoff in IEEE double-precision arithmetic:

solop = Y/.First[solop];

OrthogonalErrorPlot[solop]

The Schulz iteration, using the incremental formulation, generally yields smaller errors than the direct singular value decomposition.

Rectangular Systems

In the following example it is shown how the implementation of the orthogonal projection method also works for rectangular matrix differential systems. Formally stated, the interest is in solving ordinary differential equations on the Stiefel manifold, the set of n×p orthogonal matrices with p<n.

Definition The Stiefel manifold of n×p orthogonal matrices is the set (Subscript[V, n,p])(ℝ)=={Y∈(ℝ^(n×p))|(Y^T) Y==(Subscript[I, p])}, 1≤p<n, where (Subscript[I, p]) is the p×p identity matrix.

Solutions that evolve on the Stiefel manifold find numerous applications such as eigenvalue problems in numerical linear algebra, computation of Lyapunov exponents for dynamical systems and signal processing.

Consider an example adapted from [DL01]:

where (Subscript[q, 0])==1/(((Sqrt[n]) [1,…,1])^T), A==diag[(Subscript[a, 1]),…,(Subscript[a, n])]∈(ℝ^(n×n)), with (Subscript[a, i])==(((-1))^i)α,i==1,…,n and α>0.

The exact solution is given by

Normalizing q(t) as

it follows that Y(t) satisfies the following weak skew-symmetric system on (Subscript[V, n,1])(ℝ):

In the following example, the system is solved on the interval [0,5] with α==9/10 and dimension n==2:

p = 1;

n = 2;

α = (9/10);

ics = (1/(Sqrt[n]))Table[1,{n}];

avec = Table[(((-1))^i) α,{i, n}];

A = DiagonalMatrix[avec];

Y=Table[y[i,1][t],{i, n},{j,p}];

F = (IdentityMatrix[Length[Y]] -Y.Transpose[Y]).A;

system=Thread[Flatten[D[Y,t]]==Flatten[F.Y]];

ics=Thread[Flatten[(Y/.t->0)]==ics];

eqs={system,ics};

vars = Flatten[Y];

tfinal = 5.;

time={t,0,tfinal};

This computes the exact solution which can be evaluated throughout the integration interval:

solexact = Transpose[{(((#)/(Norm[#, 2])))}]& @ (((Exp[avec t])/(Sqrt[n])) );

This computes the solution using an explicit Runge–Kutta method:

solerk=NDSolve[eqs,vars,time,Method->"ExplicitRungeKutta"];

solerk = Y /. First[solerk];

This computes the componentwise absolute global error at the end of the integration interval:

(solexact -solerk) /. t->tfinal

This computes the orthogonal error—a measure of the deviation from the Stiefel manifold:

OrthogonalErrorPlot[solerk]

This computes the solution using an orthogonal projection method with an explicit Runge–Kutta method as the basic numerical integration scheme:

solop=NDSolve[eqs,vars,time,Method->{"OrthogonalProjection",
Method->"ExplicitRungeKutta", Dimensions->Dimensions[Y]}];

solop= Y /. First[solop];

The componentwise absolute global error at the end of the integration interval is roughly the same as before since the absolute and relative tolerances used in the numerical integration are the same:

(solexact - solop) /. t->tfinal

Using the orthogonal projection method, however, the deviation from the Stiefel manifold is reduced to the level of roundoff:

OrthogonalErrorPlot[solop]

Implementation

The implementation of the method "OrthogonalProjection" has three basic components:

Options can be used to modify the stopping criteria for the Schulz iteration. One option provided by the code is "IterationSafetyFactor" which allows control over the tolerance τ of the iteration. The factor is combined with a Unit in the Last Place, determined according to the working precision used in the integration (ULP≈2.22045×(10^(-16)) for IEEE double precision).

The Frobenius norm used for the stopping criterion can be computed efficiently using the LAPACK LANGE functions [LAPACK99].

The option (RefLink[MaxIterations, paclet:ref/MaxIterations]) controls the maximum number of iterations that should be carried out.

Option Summary

Options of the method "OrthogonalProjection".