DifferentialRiccatiEquations.ShiftsModule

This module groups all pre-defined shift strategies.

Extended help

To define a custom shift strategy, create a mutable subtype of Shifts.Strategy, define a method for Shifts.init and, optionally, methods for Shifts.update! and Shifts.take!.

struct FortyTwo <: Shifts.Strategy end

Shifts.init(::FortyTwo, _...) = FortyTwo()
Shifts.take!(::FortyTwo) = 42

If it is customary to generate multiple shift parameters at once, that are then to be used one-by-one, define a method for Shifts.take_many! and have Shifts.init return a Shifts.BufferedIterator.

struct FibonacciShifts <: Shifts.Strategy
"Number of shifts to generate at a time"
n::Int

function FibonacciShifts(n::Int)
n >= 2 || error("batch size is too small")
new(n)
end
end

struct FibonacciShiftsIterator
n::Int
f1::Int
f2::Int
end

function Shifts.init(f::FibonacciShifts)
Shifts.BufferedIterator(FibonacciShiftsIterator(f.n, 0, 1))
end

function Shifts.take_many!(it::FibonacciShiftsIterator)
n = it.n

# Generate n shifts at once:
f = Vector{Int}(undef, n)
f[1] = it.f1
f[2] = it.f2
for i in 3:n
f[i] = f[i-1] + f[i-2]
end

# Prepare next batch:
it.f1 = f[end-1] + f[end]
it.f2 = f[end] + it.f1
return f
end
DifferentialRiccatiEquations.Shifts.CyclicType
Cyclic(::Shifts.Strategy)
Cyclic(values)

Cycle through precomputed values or the shifts produced by the inner strategy. That is, continue with the first parameter once the last one has been consumed.

Examples:

Cyclic(Heuristic(10, 20, 20))
Cyclic(Float64[-1, -2, -3])
DifferentialRiccatiEquations.Shifts.HeuristicType
Shifts.Heuristic(nshifts, k₊, k₋)

Compute heuristic or sub-optimal shift parameters following Algorithm 5.1 of

Penzl: A cyclic low rank Smith method for large sparse Lyapunov equations, SIAM J. Sci. Comput., 21 (1999), pp. 1401-1418. DOI: 10.1137/S1064827598347666

DifferentialRiccatiEquations.Shifts.ProjectionType
Shifts.Projection(u::Int)

Compute shift parameters based on the u most recent increments comprising the solution candidate.

Only even u > 1 are allowed, such that an ADI double-step can properly be accounted for.

See section 5.3.1 of

Kürschner: Efficient low-rank solution of large-scale matrix equations. Otto-von-Guericke-Universität Magdeburg (2016).

The strategy has first been presented in

Benner, Kürschner, Saak: Self-generating and efficient shift parameters in ADI methods for large Lyapunov and Sylvester equations, Electronic Transactions on Numerical Analysis, 43 (2014), pp. 142-162. https://etna.math.kent.edu/volumes/2011-2020/vol43/abstract.php?vol=43&pages=142-162

DifferentialRiccatiEquations.Shifts.WrappedType
Wrapped(func!, ::Shifts.Strategy)

Apply func! to the set of shifts produced by the inner strategy via Shifts.take_many!. This strategy may be used, e.g., to filter or reorder the shifts. Complex-valued shifts must occur in conjugated pairs.

Examples:

Wrapped(reverse, Projection(2))
Wrapped(Projection(4)) do shifts
filter(s -> real(s) < -1, shifts)
end
DifferentialRiccatiEquations.Shifts.take!Method
Shifts.take!(shifts)

Return the next shift parameter from shift generator shifts.

This operation may be expensive. Compute new shift parameters, if needed.

Default: popfirst!(shifts)

DifferentialRiccatiEquations.Shifts.update!Method
Shifts.update!(shifts, X, R, Vs...)

Pass most recent solution update to shift generator shifts.

• X: current solution candidate
• R: outer factor of residual corresponding to X
• Vs: outer factors of most recent updates comprising X

This operation must be cheap. Defer the computation of new shift parameters to Shifts.take! or Shifts.take_many!.

Default: no-op.

DifferentialRiccatiEquations.CallbacksModule

This module groups all callback functions that are called at various points during the solution process. See their respective docstrings for more information. Every callback takes an observer as its first argument, which is passed via an optional keyword argument to solve.

During solve(::GALEProblem, ::LyapunovSolver; observer):

During solve(::GAREProblem, ::AlgebraicRiccatiSolver; observer):

During solve(::GDREProblem, ::Algorithm; observer):

Extended help

Hook into above callbacks by first defining a custom observer type.

mutable struct ResidualObserver
norms::Vector{Float64}
abstol::Float64

ResidualObserver() = new(Float64[], -1.0)
end

Then, create custom methods to above callbacks. If the observer needs to store any information, use some global variables (not recommended), have the observer be mutable. Note that Callbacks has to be imported manually; this is a deliberate choice.

import DifferentialRiccatiEquations.Callbacks

function Callbacks.observe_gale_step!(o::ResidualObserver, _prob, _alg, abstol::Float64, _reltol)
o.abstol = abstol
end

function Callbacks.observe_gale_step!(o::ResidualObserver, _iter, _sol, _residual, residual_norm::Float64)
push!(o.norms, residual_norm)
end

The observer is passed into the solution procedure as follows.

prob = GALEProblem(E, A, C)
obs = ResidualObserver()

solve(prob, alg; observer=obs)

@show obs.norms[end] <= obs.abstol
Todo

Update extended help to use doctests.

DifferentialRiccatiEquations.Callbacks.observe_gale_done!Method
observe_gale_done!(observer, iters::Int, X, residual, residual_norm)

Notify observer at the end of solving the GALE. The observer may compute ans store any metrics of the subsequent arguments, but it must not modify any of them.

• X: solution candidate
• residual: residual corresponding to X; usually of the same data type as X
• residual_norm: internal approximation of the norm of residual
DifferentialRiccatiEquations.Callbacks.observe_gale_metadata!Method
observe_gale_metadata!(observer, desc::String, metadata)

Notify observer on some metadata the algorithm has computed. desc gives a brief description of the metadata.

Example

The ADI calls observe_gale_metadata!(observer, "ADI shifts", μ), where μ are the (newly) computed ADI shift parameters.

DifferentialRiccatiEquations.Callbacks.observe_gale_step!Method
observe_gale_step!(observer, iter::Int, X, residual, residual_norm)

Notify observer for an iterative GALE algorithm, that iteration number iter has been completed. The observer may compute ans store any metrics of the subsequent arguments, but it must not modify any of them.

• X: solution candidate
• residual: residual corresponding to X; usually of the same data type as X
• residual_norm: internal approximation of the norm of residual
Note

The iterations iter may not be consequtive. If an algorithm computes multiple steps at once and has no (cheap) representation of the intermediate solution candidates, the difference between the values of iter of subsequent calls to observe_gale_step! may differ by more than one.

DifferentialRiccatiEquations.DRESolutionType

Solution to a generalized differential Riccati equation (DRE) as returned by solve(::GDREProblem, alg; kwargs...). The solution has three fields:

• X::Vector{T}: state X(t); T may be a Matrix or LDLᵀ
• K::Vector{<:Matrix}: feedback K(t) := B' * X(t) * E
• t::Vector{<:Real}: discretization time

By default, the state X is only stored at the boundaries of the time span, as one is mostly interested only in the feedback matrices K. To store the full state trajectory, pass save_state=true to solve.

DifferentialRiccatiEquations.GDREProblemType

Generalized differential Riccati equation

E'ẊE = C'C + A'XE + E'XA - E'XBB'XE
X(t0) = X0

having the fields E, A, C, X0, and tspan=(t0, tf).

DifferentialRiccatiEquations.LDLᵀType
LDLᵀ{TL,TD}(L::TL, D::TD)
LDLᵀ{TL,TD}(Ls::Vector{TL}, Ds::Vector{TD})

A lazy representation of L * D * L' that supports the following functions:

• +(::LDLᵀ, ::LDLᵀ) and +(::LDLᵀ{TL,TD}, ::Tuple{TL,TD})
• *(::Real, ::LDLᵀ)
• size
• rank which yields the length of the inner dimension, i.e. size(D, 1)
• zero which yields a rank 0 representation
• concatenate! (expert use only)
• compress! (expert use only)

Iterating the structure yields L::TL and D::TD. This calls compress!, if necessary.

For convenience, the structure might be converted to a matrix via Matrix. It is recommended to use this only for testing.

DifferentialRiccatiEquations.LowRankUpdateType
LowRankUpdate{TA,T,TU,TV}(A::TA, α::T, U::TU, V::TV)

Lazy representation of A + inv(α)*U*V that supports the following functions:

• \ via the Sherman-Morrison-Woodbury formula
• +(::LowRankUpdate, ::AbstractMatrix) to update A
• adjoint which returns a LowRankUpdate
• size

Iterating the structure produces the components A, α, U and V.

It is recommended to use lr_update to create a suitable representation of A + inv(α)*U*V.

DifferentialRiccatiEquations.NewtonADIType
NewtonADI()

Kleinman-Newton method to solve algebraic Riccati equations. The algebraic Lyapunov equations arizing at every Newton steps are solved using the ADI.

solve(prob::GAREProblem, NewtonADI(); kwargs...)

Supported keyword arguments:

• reltol = size(prob.A, 1) * eps(): relative Riccati residual tolerance
• maxiters = 5: maximum number of Newton steps
• observer: see Callbacks
• adi_initprev = false: whether to use previous Newton iterate as the initial guess for the ADI. If false, the default initial value of zero is used.
• adi_kwargs::NamedTuple: keyword arguments to pass to solve(_, ::ADI; adi_kwargs...)
• inexact = true: whether to allow (more) inexact Lyapunov solutions
• inexact_forcing = quadratic_forcing: compute the forcing parameter η = inexact_forcing(i, residual_norm) as described by Dembo et al. (1982), where i::Int is the Newton step and residual_norm::Float64 is the norm of the Riccati residual. See quadratic_forcing, and superlinear_forcing.
• inexact_hybrid = true: whether to switch to the classical Newton method, if the absolute Lyapunov tolerance of the classical Newton method is less strict (i.e. larger) than the tolerance η * residual_norm.
• linesearch = true: whether to perform an Armijo line search if the Riccati residual did not decrease sufficiently, see e.g. Benner et al. (2015).

Default arguments to Lyapunov solver, which can all be overwritten by adi_kwargs:

• maxiters = 100: maximum number of ADI iterations
• observer = observer
• initial_guess: see adi_initprev above
• reltol: defaults a fraction of the Riccati tolerance, reltol/10
• abstol: controlled by inexact* above, if inexact = true.

References:

• Dembo, Eisenstat, Steihaug: Inexact Newton Methods. 1982. https://doi.org/10.1137/0719025
• Benner, Heinkenschloss, Saak, Weichelt: Inexact low-rank Newton-ADI method for large-scale algebraic Riccati equations. 2015. http://www.mpi-magdeburg.mpg.de/preprints/
DifferentialRiccatiEquations.concatenate!Method
concatenate!(X::LDLᵀ)

Concatenate the internal components such that L and D may be obtained via L, D = X. This function is roughly equivalent to L = foldl(hcat, X.Ls) and D = foldl(dcat, Ds), where dcat is pseudo-code for "diagonal concatenation".

This is a somewhat cheap operation.

See also: compress!

LinearAlgebra.normMethod
norm(::LDLᵀ)

Compute the Frobenius norm of a LDLᵀ factorization. The technique is similar to the one described in

Benner, Li, Penzl. Numerical solution of large-scale Lyapunov equations, Riccati equations, and linear-quadratic optimal control problems. Numerical Linear Algebra with Applications 2008. DOI: 10.1002/nla.622

• Lang2015N Lang, H Mena, and J Saak, "On the benefits of the LDLT factorization for large-scale differential matrix equation solvers" Linear Algebra and its Applications 480 (2015): 44-71. doi:10.1016/j.laa.2015.04.006