`DifferentialRiccatiEquations.Shifts`

— ModuleThis 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.BufferedIterator`

— Type`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.

`DifferentialRiccatiEquations.Shifts.Cyclic`

— Type```
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.Heuristic`

— Type`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.Projection`

— Type`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.Wrapped`

— Type`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.init`

— Function`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.

`DifferentialRiccatiEquations.Shifts.safe_sort!`

— Method`safe_sort!(shifts)`

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

`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.take_many!`

— Function`Shifts.take_many!(generator)`

Return a `Vector{<:Number}`

of shift parameters to be used within a `Shifts.BufferedIterator`

.

`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.Callbacks`

— ModuleThis 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)
alg = ADI()
obs = ResidualObserver()
solve(prob, alg; observer=obs)
@show obs.norms[end] <= obs.abstol
```

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_failed!`

— Method`observe_gale_failed!(observer)`

Notify `observer`

that the algorithm has failed to solve the GALE. `observe_gale_done!`

will be called regardless.

`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_start!`

— Method`observe_gale_start!(observer, prob::GALEProblem, alg::LyapunovSolver)`

Notify `observer`

at the start of solving the GALE.

`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`

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.Callbacks.observe_gdre_done!`

— Method`observe_gdre_done!(observer)`

Notify `observer`

at the end of solving the GDRE.

`DifferentialRiccatiEquations.Callbacks.observe_gdre_start!`

— Method`observe_gdre_start!(observer, ::GDREProblem, ::Algorithm)`

Notify `observer`

at the start of solving the GDRE.

`DifferentialRiccatiEquations.Callbacks.observe_gdre_step!`

— Method`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`

`DifferentialRiccatiEquations.DRESolution`

— TypeSolution 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.GALEProblem`

— TypeGeneralized algebraic Lyapunov equation

`A'XE + E'XA = -C`

having the fields `A`

, `E`

and `C`

.

`DifferentialRiccatiEquations.GAREProblem`

— TypeGeneralized algebraic (continuous time) algebraic Riccati equation

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

`DifferentialRiccatiEquations.GDREProblem`

— TypeGeneralized 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)`

.

`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.LowRankUpdate`

— Type`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.NewtonADI`

— Type`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.compress!`

— Method`compress!(X::LDLᵀ)`

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

This is an expensive operation.

See also: `concatenate!`

`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!`

`DifferentialRiccatiEquations.lr_update`

— Function```
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`

.

`DifferentialRiccatiEquations.quadratic_forcing`

— Method`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`

.

`DifferentialRiccatiEquations.superlinear_forcing`

— Method`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`

.

`LinearAlgebra.norm`

— Method`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