API

Orbit(a, p)

Unknowns of the search: the space-time field a, in any discretisation, and the parameter vector p = [log ω, c₁, c₂, …], the logarithm of the frequency ω = 2π/T followed by one drift speed per drift direction of the System. The frequency is searched through its logarithm, which keeps it positive and makes a step a relative change of ω; frequency returns ω. Orbits broadcast like vectors; dot adds the inner product of the fields to that of the parameters.

x = Orbit(a, [log(2π/T), c])

frequency(x::Orbit) -> ω

The frequency ω = 2π/T of the orbit, from its logarithm x.p[1].

System(nl, lin, adj, dds!, x::Orbit; linearise!, ddi=(), B=I, project=identity)

The system with nonlinear, linearised and adjoint operators nl, lin and adj, rescaled-time derivative dds! and function linearise!(op, a) setting the linearisation point of lin or adj, on orbits shaped like x, e.g. for the Kuramoto–Sivashinsky example

System(KSNonlinear(g), KSLinearised(g), KSLinearised(g; adjoint=true), dds!, x;
       linearise!, ddi=(ddx!,))

ddi is a tuple of functions dd!(out, a) -> out, the derivatives along the directions of translational symmetry in which the orbit drifts with an unknown speed, one per drift speed in x.p. The default () searches for periodic orbits. Leave out the directions in which a symmetry of the problem forbids drifting.

B is the preconditioner, an invertible operator on orbits for which precondition! and precondition_adjoint! are defined; the identity I by default.

project(a) -> a restricts the search to a subspace of fields invariant under the operators, the fields fixed by a symmetry of the problem: it projects the field a in place, orthogonally, onto the subspace. It is applied to the initial orbit, and to every residual, gradient, Jacobian action and application of the preconditioner, so that the search stays in the subspace despite rounding errors. The default identity does nothing.

solve!(x::Orbit, F::System, method) -> x

Search for an orbit of the system F, starting from x and overwriting it with the result. method is an LBFGS object, which minimises ½‖r‖² and is robust far from a solution, or a NewtonHookstep object, which solves r = 0 and converges fast close to one. Both hold their own options. A typical search runs L-BFGS first and refines with the hookstep:

solve!(x, F, LBFGS(maxiter=300))
solve!(x, F, NewtonHookstep(krylov_dim=150))

LBFGS(; memory=10, maxiter=100, tol=1e-10, verbose=true, io=stdout,
        callback=info -> false)

L-BFGS minimisation of ½‖r‖², keeping memory curvature pairs, for solve!. Stops when ‖r‖ < tol, after maxiter iterations, or when callback(info), called at the start and after every iteration, returns true. info is a named tuple with the iteration iter, the orbit x, the residual norm res, the cumulative evaluations of the system and krylov, empty here. A Trace records them. Robust far from a solution.

NewtonHookstep(; krylov_dim=50, krylov_tol=1e-3, Δ=1e-2, Δmax=1, eta=1e-3,
                 maxiter=20, tol=1e-10, verbose=true, io=stdout,
                 callback=info -> false)

Newton–Krylov hookstep solution of r = 0, for solve!. Krylov spaces grow up to dimension krylov_dim, or until the GMRES solution of the Newton system has a residual below krylov_tol times ‖r‖, and the hookstep is then computed in that space; the trust-region radius starts at Δ and never exceeds Δmax; steps are accepted when the actual reduction of ½‖r‖² is at least eta times the predicted one. Stops when ‖r‖ < tol, after maxiter iterations, or when callback(info), called at the start and after every iteration, returns true. info is a named tuple with the iteration iter, the orbit x, the residual norm res, the cumulative evaluations of the system, krylov, the relative residual of the GMRES solution of the Newton system after each Arnoldi step of the iteration, and step, "newton" or "hook", the kind of the accepted step. A Trace records them. Fast near a solution.

Trace()

History of a search, filled when passed as callback to LBFGS or NewtonHookstep. At the start and after every iteration it records:

• iter, the iteration number of the method;
• res, the residual norm ‖r‖;
• p, the parameters [log ω, c₁, …];
• time, the seconds elapsed since the first record;
• evaluations, the cumulative numbers of residuals, gradients, Jacobian actions and preconditioner applications of the system;
• krylov, for the hookstep, the relative residual of the linear model after each Arnoldi step of the iteration; empty for L-BFGS and at the start;
• step, for the hookstep, "newton" if the accepted step is the full Newton step, inside the trust region, or "hook" if it lies on its boundary; empty for L-BFGS and at the start.

The same trace can follow several calls to solve!, e.g. L-BFGS and then the hookstep: the iteration numbers restart at zero with each call.

trace = Trace()
solve!(x, F, LBFGS(callback=trace))
solve!(x, F, NewtonHookstep(callback=trace))

trace.res # residual norms, in order

precondition!(out::Orbit, B, p::Orbit) -> out

Apply the inverse of the preconditioner B to the orbit p: out = B⁻¹ p. Extend it, together with precondition_adjoint!, for each preconditioner type.

precondition_adjoint!(out::Orbit, B, p::Orbit) -> out

Apply the adjoint of the inverse of the preconditioner B to the orbit p: out = B⁻⁺ p, in the inner product of the orbits. For a self-adjoint B it is precondition!.

residual!(r::Orbit, F::System, x::Orbit) -> r

The augmented residual at x: the field r.a = ω ∂s a - Σᵢ cᵢ ∂ᵢ a - N(a), and zero parameters r.p, the phase conditions, which vanish at the current point. Minus r is the right-hand side of the Newton system.

objective(F::System, x::Orbit) -> R

R = ½‖r‖². The residual is left in F.residual.

gradient!(g::Orbit, F::System, x::Orbit) -> R

Write the gradient of R = ½‖r‖² into g and return R. Moves the linearisation point of the operators to x.

jacobian!(out::Orbit, F::System, x::Orbit, δx::Orbit) -> out

The Jacobian of the residual at x, with the phase conditions as extra rows, applied to δx. Assumes the linearisation point was set to x with F.linearise!(F.lin, x.a).

jacobian_adjoint!(out::Orbit, F::System, x::Orbit, w::Orbit) -> out

The adjoint of jacobian! in the inner product of the orbits, applied to w: with w = (w_a, [w_s, w₁, …]), the field and the multipliers of the phase conditions,

out.a = -ω ∂s w_a + Σᵢ cᵢ ∂ᵢ w_a - L⁺ w_a + w_s ∂s a + Σᵢ wᵢ ∂ᵢ a
out.p = [ω ⟨∂s a, w_a⟩, -⟨∂₁ a, w_a⟩, …]

so that ⟨w, 𝒥 v⟩ = ⟨𝒥⁺ w, v⟩ for all orbits v. Assumes the linearisation point of the adjoint operator was set to x with F.linearise!(F.adj, x.a). Needed by Krylov solvers of the least-squares problem min ‖𝒥 δx + r‖, such as LSQR; the gradient of R = ½‖r‖² is 𝒥⁺ (r, 0).