This module groups all pre-defined shift strategies.

• Shifts.Heuristic
• Shifts.Projection
• Shifts.Cyclic
• Shifts.Wrapped

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

BufferedIterator(generator)

Initialize an internal buffer of type Vector{<:Number} from Shifts.take_many!(generator) and return shifts one-by-one using popfirst!. Refill the buffer once it is depleated.

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])

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

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

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

Shifts.init(::Shifts.Strategy, prob)

Create and initialize a shift generator from problem data. The returned iterator will immediately be Shifts.update!ed with initial guess and residual of the iteration.

safe_sort!(shifts)

Ensure that complex conjugated values are located adjacent to one another.

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)

Shifts.take_many!(generator)

Return a Vector{<:Number} of shift parameters to be used within a Shifts.BufferedIterator.

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.

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):

• observe_gale_start!
• observe_gale_step!
• observe_gale_done!
• observe_gale_failed!
• observe_gale_metadata!

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

• observe_gare_start!
• observe_gare_step!
• observe_gare_done!
• observe_gare_failed!

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

• observe_gdre_start!
• observe_gdre_step!
• observe_gdre_done!

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)
alg = ADI()
obs = ResidualObserver()

solve(prob, alg; observer=obs)

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

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

observe_gale_failed!(observer)

Notify observer that the algorithm has failed to solve the GALE. observe_gale_done! will be called regardless.

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.

observe_gale_start!(observer, prob::GALEProblem, alg::LyapunovSolver)

Notify observer at the start of solving the GALE.

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

observe_gdre_done!(observer)

Notify observer at the end of solving the GDRE.

observe_gdre_start!(observer, ::GDREProblem, ::Algorithm)

Notify observer at the start of solving the GDRE.

observe_gdre_step!(observer, t::Float64, X, K)

Notify observer that the step to time point t has been completed.

• X: solution at time t
• K: feedback matrix K = B' * X * E where B denotes the input map of the associated GDREProblem

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.

Generalized algebraic Lyapunov equation

A'XE + E'XA = -C

having the fields A, E and C.

Generalized algebraic (continuous time) algebraic Riccati equation

Q + A'XE + E'XA - E'XGXE = 0

Generalized differential Riccati equation

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

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

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.

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.

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/

compress!(X::LDLᵀ)

Concatenate the internal components and perform a column compression following ^[Lang2015].

This is an expensive operation.

See also: concatenate!

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!

lr_update(A::Matrix, α, U, V)
lr_update(A::AbstractSparseMatrixCSC, α, U, V)

Return a suitable representation of A + inv(α)*U*V. For dense A, compute A + inv(α)*U*V directly. For sparse A, return a LowRankUpdate.

quadratic_forcing(_, residual_norm) = min(0.1, 0.9 * residual_norm)

Exemplary forcing term to obtain quadratic convergence in the inexact Newton method. residual_norm::Float64 refers to the norm of the previous Newton residual. See NewtonADI.

superlinear_forcing(i, _) = 1 / (i^3 + 1)

Exemplary forcing term to obtain superlinear convergence in the inexact Newton method. i::Int refers to the current Newton step. See NewtonADI.

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