API Reference

Adam and the methods that build their direction from a pair of bias-corrected moment averages as it does. What the union selects for is that shared character and not a shared implementation: the direction is a moving average that is deliberately allowed not to descend on an individual step, which is what default_linesearch's fixed Static and the ensure_descent! exemption in FirstOrderMethodWithState are about.

The implementation is shared by only some of them, so do not read the union as a promise about the cache:

• AdamWithEuclideanDecay differs only in how the direction is finished off, so it reuses AdamCache and AdamState outright and OptimizerState(AdamWithEuclideanDecay(), x) returns an AdamState.
• ScalarMomentAdam has neither, because its second moment is a number rather than an element of 𝔤ʰᵒʳ and its recursion is therefore not Adam's. It carries ScalarMomentAdamCache and ScalarMomentAdamState of its own.

CURVATURE_TOLERANCE

How positive $\delta^T\gamma$ has to be, relative to $\|\delta\|\|\gamma\|$, for the secant pair to be used. See curvature_is_usable.

DEFAULT_STEP_CEILING

How far along its direction a manifold solve may step, in multiples of $2\pi$: the c of step_αmax(c, δ), and hence the params.αmax solver_step! hands the line search through linesearch_parameters. Set per solve with Optimizer(x, F; step_ceiling = …); Inf switches the ceiling off.

Why a ceiling is needed at all, and why the caller has to set it

A line search bounds the step it returns by the merit — it stops when $\varphi$ stops falling. On a compact manifold $\varphi$ is bounded, so that test never fires: at $\alpha = 10^9$, $\varphi$ can be genuinely lower than at $\alpha = 0$, and a search that reports a decrease there is telling the truth. Measured on the SVD problem of test/optimizer_convergence/svd_optim.jl, SimpleSolvers.Quadratic returned $\alpha = 4.3\times10^7$ on a direction of norm 5.54 — a step of $\|\alpha\delta\| = 2.4\times10^8$ — and did it again two steps later on a steepest-descent direction, so the behaviour belonged to the search and not to what it was handed. Retracting a lift of that norm loses the manifold, and the solve then reported convergence from a point that was no longer on $St(20, 3)$. That was issue A1b.

The merit is not a bound on the step, and the bound that does exist is not a property of $\varphi$ at all: it is the $2\pi$ of a rotation divided by $\|\delta\|$, it changes at every solver step, and nothing a line search can measure reveals it. SimpleSolvers 0.12 therefore splits the ceiling in two — a method's own DEFAULT_LINESEARCH_αmax backstop, and a per-call params.αmax — and leaves this half here deliberately. Upstream measured why the backstop alone cannot close it: at $\|\delta\| \approx 5.5$ its 65536 still permits $\|\alpha\delta\| = 3.6\times10^5$, five orders above the $2\pi$ past which a rotation stops going anywhere, and four of the eight starting points still diverge.

Why 1

c is the ceiling in multiples of $2\pi$, so 1 says "never step further along the direction than one full turn". That is the scale the geometry gives; the value is not tuned to the SVD problem beyond confirming that it leaves it alone. Over the converging solve next to the diverging one — seed 2 at 90 iterations — the largest $\|\alpha\delta\|$ is 2.03, well inside $2\pi$, which is what makes a ceiling that forbids the $10^8$ step cost the ordinary ones nothing.

The methods whose update! needs the method rather than a SimpleSolvers.Hessian, because they carry state of their own out of which the direction is built — a momentum or a pair of moments. The (quasi-)Newton methods are the complement: their direction comes from the Hessian and the method object holds nothing.

GradientMethod is in neither group: it has no Hessian and no state, and takes the Hessian branch because NoHessian is all that branch needs from it.

solver_step! also uses this to decide whether ensure_descent! applies. It must not: a momentum and a moment average are deliberately allowed not to descend on an individual step, and the decay of AdamWithEuclideanDecay tilts the direction further away from the gradient still.

That exemption is about the direction before the search, and it used to be extended to the restart after it as well. It no longer is — see linesearch_rejected and issue A7. Declining to overrule the direction is not the same as taking the longest step along it, which is what a rejected search hands back.

OptimizerSolution

A type alias for the solution of an optimizer, which can be either an AbstractVector or a Manifold.

VectorStorageMatrix

The structured matrix types that keep their free parameters in one vector, reachable as parent(A): SkewSymMatrix, SymmetricMatrix and the two AbstractTriangulars.

Each of them stores $n(n\pm1)/2$ numbers behind an $n \times n$ interface, so none of the generic AbstractArray methods, which either broadcast or reshape, is usable for them:

• three of the four have no setindex! for a broadcast to write through. SkewSymMatrix and the two triangulars cannot represent an arbitrary matrix at all, so there is no entry to assign. SymmetricMatrix is the exception — it does have one (writing $[A]_{ij}$ writes $[A]_{ji}$ too, which is what a symmetric matrix means), so the broadcast would give the right answer for that one type, at twice the work and only as long as it stays the exception.
• similar has to preserve the type, because the optimizer caches allocate their scratch with it and then require every array to have the same type as the parameter. The AbstractArray fallback returns a dense Matrix.
• ParameterHandling.flatten has to round trip through the free parameters. Its AbstractMatrix method reshapes the flattened vector back to $n \times n$, and $n(n\pm1)/2$ numbers do not reshape to that.

For every one of these types the free parameters are the coordinates the optimizer should work in, so each primitive is the corresponding operation on parent. _add!, _rac!, _square!, _div!, _rmul!, _difference! and flatten live next to the rest of their family in optimizers/named_tuple_wrapper.jl, l2norm in optimizers/optimizer_status.jl; update_section! in global_sections/global_sections.jl splits on this alias for the same reason.

This is what lets a SymmetricMatrix or a triangular matrix be an optimizer parameter — which is what GeometricMachineLearning's SympNet, symplectic-attention and volume-preserving layers need, and what this package could not do before.

The alias lives in a file of its own because it is needed by global_sections.jl, which is included long before the primitives are.

Matrix(λY::GlobalSection)

Put λY into matrix form.

This is not recommended if speed is important!

Use apply_section and global_rep instead!

AbstractExponentialAlgorithm

Supertype of the algorithms that Geodesic can use to evaluate the matrix exponential.

A horizontal lift factors as $\bar{B} = B'(B'')^T$ into two $N\times{}2n$ matrices, and the exponential of a product in that order is

\[\exp(B'(B'')^T) = \mathbb{I} + B'\,\mathfrak{A}((B'')^TB')\,(B'')^T, \qquad \mathfrak{A}(X) = \sum_{n=1}^\infty \frac{X^{n-1}}{n!},\]

so the whole computation reduces to one $2n\times{}2n$ matrix function — the argument $X = (B'')^TB'$ is small even when $N$ is large. $\mathfrak{A}$ is the function usually written $\varphi_1(X) = (\exp(X) - \mathbb{I})X^{-1}$, though it is defined by the series and is perfectly regular at a singular $X$.

The subtypes differ only in how that $2n\times{}2n$ function is evaluated, and they are mathematically identical — every one of them returns the exponential map, so the one-parameter subgroup property Geodesic relies on holds for all of them. They differ in accuracy at a large lift, in cost, and in which backends they run on. On a random StiefelLieAlgHorMatrix(20, 3) scaled to $\|\bar{B}\| = 361$, against the cost of one retraction at $N = 200$, $n = 10$:

"Forward error" is the relative distance to exp(Matrix(B)); both it and check come from the same lift, so the columns are comparable row by row. The cost column is one whole retraction, most of which is the $N\times{}N$ assembly they share; the $\mathfrak{A}$ call on its own, which is what actually separates the three algorithms that evaluate it, is 0.021 ms and 201 KiB for ScaledSquaring, 0.037 ms and 330 KiB for NativePade, and 0.053 ms and 114 KiB for AugmentedPade — note that the delegating one, which discards three quarters of its work, is nonetheless the lightest allocator, because Base.exp reuses buffers where these two build a fresh $2n\times{}2n$ temporary per operation.

ScaledSquaring remains the default because it is the cheapest portable algorithm. NativePade is the independent portable cross-check, ProjectedSkew if check matters more than the last digit of the exponential, and AugmentedPade as the CPU reference that delegates its numerics to Base.exp. TaylorSeries is the pre-0.2.0 behaviour and is retained only so the regression is reproducible; it is not a usable retraction.

A new one has to supply 𝔄(X::AbstractMatrix, ::NewAlgorithm); geodesic and everything above it then follow. An algorithm that does not go through $\mathfrak{A}$ at all — ProjectedSkew is the one such — supplies geodesic(::AbstractLieAlgHorMatrix, ::NewAlgorithm) instead.

See GeometricOptimizers.𝔄 for the implementations.

AbstractLieAlgHorMatrix <: AbstractMatrix

AbstractLieAlgHorMatrix is a supertype for various horizontal components of Lie algebras. We usually call this $\mathfrak{g}^\mathrm{hor}$.

See StiefelLieAlgHorMatrix and GrassmannLieAlgHorMatrix for concrete examples.

AbstractOptimizerProblem

See SimpleSolvers.LinesearchProblem and OptimizerProblem.

AbstractRetraction

Supertype of the retractions, i.e. of the maps that turn an element of a horizontal Lie algebra back into a point of the manifold.

A retraction is how every step in this package is taken: an OptimizerMethod produces a direction $B$ in the horizontal component of $\mathfrak{g}$, and $\Lambda \mapsto \Lambda\cdot\mathrm{retract}(B)$ moves the global section — and with it the point — without ever leaving the manifold. An instance is passed to Optimizer as retraction = Cayley(); the default is Cayley.

Concrete subtypes are Geodesic and Cayley. A new one has to supply retraction(::NewRetraction, x); the callable form R(x) and everything in update_section! then follow.

For every manifold layer one has to specify a retraction that takes the layer and elements of the (global) tangent space.

AbstractTriangular

See UpperTriangular and LowerTriangular.

Adam(T; β₁, β₂, δ)

The Adam optimizer, with the defaults suggested in [24, page 301].

The cache consists of a first and a second moment, stored in bias-corrected form, i.e. updated in the $t$-th iteration (counted from $t = 1$) as

\[m_1 \gets \frac{\beta_1 - \beta_1^t}{1 - \beta_1^t}m_1 + \frac{1 - \beta_1}{1 - \beta_1^t}\nabla{}L,\]

\[m_2 \gets \frac{\beta_2 - \beta_2^t}{1 - \beta_2^t}m_2 + \frac{1 - \beta_2}{1 - \beta_2^t}\nabla{}L\odot\nabla{}L,\]

from which the direction is computed as $-m_1/(\sqrt{m_2} + \delta)$.

T is the element type of the parameters that are to be optimized; unlike MomentumMethod, Adam is not converted by Optimizer, so Adam(Float32) is needed for Float32 parameters.

AdamCache <: OptimizerCache

Cache for the gradient optimizer.

Fields

• x::Manifold: the solution,
• g: the gradient (for the manifold case this is in AbstractLieAlgHorMatrix form),
• δ: the direction,
• Δg: difference in gradients, needed for OptimizerStatus,
• g̃: scratch for latest_gradient; see GradientCache,
• g̃_is_current: whether g̃ is the gradient at x; see store_gradient!,
• section: the GlobalSection.

AdamState <: OptimizerState

State for the gradient optimizer.

AdamWithEuclideanDecay(T; β₁, β₂, δ, λ)

Adam with the decoupled weight decay of [28], which is the Euclidean $\lambda{}x$ and therefore acts on the unconstrained weights only.

The moments are the ones of Adam — this method shares its cache and its state — and the decay is applied to the direction instead of to the gradient, i.e. the direction is

\[-\left(\frac{m_1}{\sqrt{m_2} + \delta} + \lambda{}x\right)\]

so that a step with learning rate $\eta$ (see below) is

\[x \gets x - \eta\frac{m_1}{\sqrt{m_2} + \delta} - \eta\lambda{}x.\]

That is what decoupled means: $\lambda{}x$ never enters $m_1$ or $m_2$, so the amount by which a weight is shrunk does not depend on the size of its gradient. Adding $\lambda{}x$ to the gradient instead — $L^2$ regularization, which is what plain Adam on a penalized objective does — is not the same thing, because the second moment then rescales the penalty away.

Why the name says Euclidean

$\lambda{}x$ is the gradient of $\frac{\lambda}{2}||x||^2$, which is constant on the StiefelManifold and the GrassmannManifold — both are compact, with $||Y||_F^2 = \mathrm{tr}(Y^TY) = n$ — so its Riemannian gradient $\mathtt{rgrad}(Y, \lambda{}Y)$ vanishes identically and the decay does nothing at all to a manifold weight. The method therefore decays the ordinary arrays of a NamedTuple of parameters and leaves the Manifolds in it alone — which is the case it exists for, a network that keeps Stiefel weights next to unconstrained ones — and on a bare Manifold it is Adam, for every $\lambda$. Passing a nonzero λ together with parameters that are entirely manifolds is therefore warned about rather than silently ignored.

The derivation, the Grassmann case and why the name AdamW is held in reserve for a Riemannian decay instead (see issue #28) are in the Weight Decay on Manifolds page.

Arguments

As for Adam, T is the element type of the parameters (AdamWithEuclideanDecay(Float32) for Float32 parameters) and the learning rate $\eta$ is not stored here: it is the α of the line search, i.e. linesearch = Static(η), which is also the default (see default_linesearch). λ is hence multiplied by $\eta$, exactly as in [28], where the two are decoupled from each other but the decay is still scaled by the schedule.

AugmentedPade <: AbstractExponentialAlgorithm

Evaluate $\mathfrak{A}$ as a block of a larger ordinary exponential.

\[\exp\begin{pmatrix} X & \mathbb{I} \\ \mathbb{O} & \mathbb{O} \end{pmatrix} = \begin{pmatrix} \exp(X) & \mathfrak{A}(X) \\ \mathbb{O} & \mathbb{I} \end{pmatrix}\]

so one call to Base.exp on a $4n\times{}4n$ matrix returns $\mathfrak{A}(X)$ in its upper-right block. That hands the numerics to Julia's own exponential — a degree-13 Padé approximant with its own scaling and squaring, and the most heavily exercised implementation available — at the cost of exponentiating a matrix four times the size and discarding three quarters of it.

Accurate to the same order as ScaledSquaring, and about twice as expensive in the $\mathfrak{A}$ call itself though much less than that once the $N\times{}N$ assembly is counted. Its value is that it introduces no new numerics: it is the reference the other algorithms are tested against in test/retractions/exponential_accuracy.jl.

See AbstractExponentialAlgorithm for the alternatives.

BFGS

The Broyden–Fletcher–Goldfarb–Shanno method, from [27], and the default algorithm of Optimizer. The direction is $-Q\nabla{}f(x)$ for an approximate inverse Hessian $Q$, updated from the secant pair $\delta = x^{(k)} - x^{(k-1)}$, $\gamma = \nabla{}f^{(k)} - \nabla{}f^{(k-1)}$ by

\[Q \gets Q - \frac{\delta\gamma^TQ + Q\gamma\delta^T - \left(1 + \frac{\gamma^TQ\gamma}{\delta^T\gamma}\right)\delta\delta^T}{\delta^T\gamma}\]

Its state is BFGSState and the update itself lives in BFGSCache, which skips the update when the secant pair fails curvature_is_usable — the condition that keeps $Q$ positive definite. restart!(::BFGSState) discards $Q$ when a line search reports that it could not decrease the merit.

Unlike Newton this builds its curvature up over the iterations rather than evaluating a Hessian, so it takes more of them; unlike DFP it produces a direction already scaled like a Newton step, and accepts $\alpha = 1$ on most of them. Both run on an AbstractVector, a NamedTuple of parameters and a bare Manifold alike: $Q$ is sized by the intrinsic dimension and the secant pair is taken in the horizontal lift, so nothing in the update needs a vector-valued point.

Takes no parameters: like every OptimizerMethod it only produces a direction, and how far the optimizer travels along it is the line search's business.

BFGSCache

The OptimizerCache for the BFGS algorithm. Also see update!(::BFGSCache, ::OptimizerState, ::AbstractVector, ::AbstractVector).

g̃ is the scratch for latest_gradient and g̃_is_current says whether it is the gradient at x; see GradientCache, which carries the same pair for the same reason, and store_gradient!.

BFGSState <: OptimizerState

The OptimizerState corresponding to the BFGS method.

Keys

• x̄
• s: stores the previous direction. This needs to be stored in addition to the previous solution because of the manifold case.
• ḡ
• f̄
• Q

Cayley <: AbstractRetraction

The Cayley transform as a retraction, and the default one.

\[\mathrm{Cayley}(B) = \left(\mathbb{I} - \frac{1}{2}B\right)^{-1}\left(\mathbb{I} + \frac{1}{2}B\right)\]

For $B$ in the horizontal component this maps into the Lie group, so the retracted point stays on the manifold to round-off. cayley never forms the $N\times{}N$ inverse: it factors $B = B'(B'')^T$ into two $N\times{}2n$ matrices with lift_factors and inverts a $2n\times{}2n$ matrix instead.

Examples

using GeometricOptimizers
using GeometricOptimizers: Cayley, check

Y = rand(StiefelManifold, 5, 3)
B = GeometricOptimizers.global_rep(GlobalSection(Y), rand(5, 3))
check(Cayley()(B)) < 1e-14        # the retracted element is still on the manifold

# output

true

See cayley for the implementation and Geodesic for the alternative.

DFP

The Davidon–Fletcher–Powell method, from [27]: the direction is $-Q\nabla{}f(x)$ for an approximate inverse Hessian $Q$, updated from the secant pair $\delta = x^{(k)} - x^{(k-1)}$, $\gamma = \nabla{}f^{(k)} - \nabla{}f^{(k-1)}$ by

\[Q \gets Q - \frac{Q\gamma\gamma^TQ}{\gamma^TQ\gamma} + \frac{\delta\delta^T}{\delta^T\gamma}\]

(equation 6.15 there). Its state is BFGSState, under the alias DFPState, and the update itself lives in DFPCache, which symmetrizes $Q\gamma\gamma^TQ$ explicitly and skips the update when the secant pair fails curvature_is_usable.

Where BFGS is the better-conditioned choice and the default, DFP is here because it is the other classical member of the family and because the two differ in exactly the way that shows what the line search is doing: DFP's direction is systematically under-scaled, so it wants a median $\alpha$ well above 1 and is starved by a search that cannot exceed it. See default_linesearch, which measures this, converges DFP under the default expanding Backtracking, and gives StrongWolfe(T; c₂ = 0.1) as the steadier explicit choice for a DFP-heavy workload.

Takes no parameters: like every OptimizerMethod it only produces a direction, and how far the optimizer travels along it is the line search's business.

DFPCache <: OptimizerCache

The OptimizerCache corresponding to the DFP method.

g̃ is the scratch for latest_gradient and g̃_is_current says whether it is the gradient at x; see GradientCache, which carries the same pair for the same reason, and store_gradient!.

DFPState

This is equivalent to BFGSState.

DecayingStatic(T; η₁, η₂, n)

A SimpleSolvers.LinesearchMethod that takes no search at all, like SimpleSolvers.Static, but whose step size decays geometrically with the iteration number:

\[\alpha(t) = \gamma^t\eta_1, \qquad \gamma = \exp(\log(\eta_2/\eta_1)/n),\]

so that $\alpha(0) = \eta_1$ and $\alpha(n) = \eta_2$. It keeps decaying past n; the decay is not floored at $\eta_2$, because a floor is exactly what stops the iteration from converging (see below).

$\alpha(0)$ is a limit of the schedule and not a step a solve takes: solve! calls increase_iteration_number! before solver_step!, so the first step is $\alpha(1) = \gamma\eta_1$, the $k$-th is $\alpha(k)$, and the horizon is reached exactly at iteration n rather than one before it. $\eta_1$ is therefore an upper bound on the step rather than the first one — by a factor of $\gamma$, which for the defaults is 0.9908. The offset is deliberate: it is also how GeometricMachineLearning's AdamOptimizerWithDecay counts, which is what makes AdamOptimizerWithDecay reproduce it step for step (test/adam_optimizer_with_decay.jl asserts this against a live solve).

Why this exists

Adam produces a direction of magnitude $\approx{}1$ per component whatever the gradient is — that scale-freeness is the point of it. With a constant step $\alpha$ the iterate therefore does not converge to the minimizer but circles it at a distance of order $\alpha$, so no gradient-based stopping criterion is ever met and the solve runs until max_iterations. A step size that goes to zero is what turns it back into a convergent iteration.

This is what the AdamWithDecay method of v0.1.0 did with its own η₁, η₂ and n_epochs fields. It was removed when the step size moved out of the OptimizerMethods and into the line search; this is the replacement, and being a LinesearchMethod it composes with GradientMethod and MomentumMethod just as well.

Examples

ls = DecayingStatic(; η₁ = 1e-2, η₂ = 1e-6, n = 1000)

# output

DecayingStatic from α = 0.01 to α = 1.0e-6 over 1000 iterations.

Geodesic(algorithm = ScaledSquaring()) <: AbstractRetraction

The exponential map as a retraction, i.e. the true geodesic of the manifold.

\[\mathrm{Geodesic}(B) = \exp(B)\]

Because this is the matrix exponential, $\alpha \mapsto \mathrm{Geodesic}(\alpha{}B)$ is a one-parameter subgroup: it follows the geodesic through the point in the direction $B$ exactly, and $\mathrm{Geodesic}((\alpha+\beta)B) = \mathrm{Geodesic}(\alpha{}B)\mathrm{Geodesic}(\beta{}B)$. That is the property Cayley lacks, and the reason a derivative-based line search is exact here.

geodesic exploits the sparsity of a horizontal lift rather than exponentiating the full $N\times{}N$ matrix: the only matrix function it evaluates is on a $2n\times{}2n$ argument. Since 0.2.0 that also makes it the cheaper of the two for $N \gtrsim 50$ — see the note on Cayley.

algorithm selects how the exponential is evaluated. All of them compute the same map — the choice is one of accuracy at a large lift, cost, and backend support — and the default ScaledSquaring is the one to use unless there is a reason not to. See AbstractExponentialAlgorithm for the comparison.

Examples

using GeometricOptimizers
using GeometricOptimizers: Geodesic, check

Y = rand(StiefelManifold, 5, 3)
B = GeometricOptimizers.global_rep(GlobalSection(Y), rand(5, 3))
check(Geodesic()(B)) < 1e-14      # the retracted element is still on the manifold

# output

true

A large lift is where the algorithms part company, and the reason the default is what it is:

using GeometricOptimizers
using GeometricOptimizers: Geodesic, ScaledSquaring, TaylorSeries, check
import Random
Random.seed!(1234)

B = 60 * rand(StiefelLieAlgHorMatrix, 20, 3)      # ‖B̄‖ ≈ 393

check(Geodesic()(B)) < 1e-12, check(Geodesic(TaylorSeries())(B)) < 1e-12

# output

(true, false)

See geodesic for the implementation and Cayley for the cheaper alternative.

GlobalSection(Y)

Construct a global section for Y.

A global section $\lambda$ is a mapping from a homogeneous space $\mathcal{M}$ to the corresponding Lie group $G$ such that

\[\lambda(Y)E = Y,\]

Also see apply_section and global_rep.

Implementation

For an implementation of GlobalSection for a custom array (especially manifolds), the function global_section has to be generalized.

GradientCache <: OptimizerCache

Cache for the gradient optimizer.

Fields

• x::Manifold: the solution,
• g: the gradient (for the manifold case this is in AbstractLieAlgHorMatrix form),
• δ: the direction,
• Δg: difference in gradients,
• g̃: scratch for latest_gradient,
• g̃_is_current: whether g̃ is the gradient at x; see store_gradient!,
• section: the GlobalSection.

Implementation

The gradient descent algorithm.

GradientState <: OptimizerState

State for the gradient optimizer.

GrassmannLieAlgHorMatrix(B::AbstractMatrix, N::Integer, n::Integer)

Build an instance of GrassmannLieAlgHorMatrix based on an arbitrary matrix B of size $(N-n)\times{}n$.

GrassmannLieAlgHorMatrix is the horizontal component of the Lie algebra of skew-symmetric matrices (with respect to the canonical metric).

Extended help

The projection here is: $\pi:S \to SE/\sim$ where

\[E = \begin{bmatrix} \mathbb{I}_{n} \\ \mathbb{O}_{(N-n)\times{}n} \end{bmatrix},\]

and the equivalence relation is

\[V_1 \sim V_2 \iff \exists A\in\mathcal{S}_\mathrm{skew}(n) \text{ such that } V_2 = V_1 + \begin{bmatrix} A \\ \mathbb{O} \end{bmatrix}\]

An element of GrassmannLieAlgMatrix takes the form:

\[\begin{pmatrix} \bar{\mathbb{O}} & B^T \\ B & \mathbb{O} \end{pmatrix},\]

where $\bar{\mathbb{O}}\in\mathbb{R}^{n\times{}n}$ and $\mathbb{O}\in\mathbb{R}^{(N - n)\times(N-n)}.$

GrassmannLieAlgHorMatrix(D::AbstractMatrix, n::Integer)

Take a big matrix as input and build an instance of GrassmannLieAlgHorMatrix.

The integer $N$ in $Gr(n, N)$ here is the number of rows of D.

Extended help

If the constructor is called with a big $N\times{}N$ matrix, then the projection is performed the following way:

\[\begin{pmatrix} A & B_1 \\ B_2 & D \end{pmatrix} \mapsto \begin{pmatrix} \bar{\mathbb{O}} & -B_2^T \\ B_2 & \mathbb{O} \end{pmatrix}.\]

This can also be seen as the operation:

\[D \mapsto \Omega(E, DE - EE^TDE),\]

where $\Omega$ is the horizontal lift GeometricOptimizers.Ω.

GrassmannManifold <: Manifold

The GrassmannManifold is based on the StiefelManifold.

HessianBFGS <: Hessian

A struct derived from SimpleSolvers.Hessian to be used for an Optimizer.

HessianDFP <: Hessian

The SimpleSolvers.Hessian corresponding to the DFP method.

IterativeHessian <: Hessian

An abstract type derived from SimpleSolvers.Hessian. Its main purpose is defining a supertype that encompasses HessianBFGS and HessianDFP for dispatch.

LowerTriangular(S::AbstractVector, n::Int)

Build a lower-triangular matrix from a vector.

A lower-triangular matrix is an $n\times{}n$ matrix that has zeros on the diagonal and on the upper triangular.

The data are stored in a vector $S$ similarly to other matrices. See UpperTriangular, SkewSymMatrix and SymmetricMatrix.

The struct two fields: S and n. The first stores all the entries of the matrix in a sparse fashion (in a vector) and the second is the dimension $n$ for $A\in\mathbb{R}^{n\times{}n}$.

Examples

using GeometricOptimizers
S = [1, 2, 3, 4, 5, 6]
LowerTriangular(S, 4)

# output

4×4 LowerTriangular{Int64, Vector{Int64}}:
 0  0  0  0
 1  0  0  0
 2  3  0  0
 4  5  6  0

LowerTriangular(A::AbstractMatrix)

Build a lower-triangular matrix from a matrix.

This is done by taking the lower left of that matrix.

Examples

using GeometricOptimizers
M = [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
LowerTriangular(M)

# output

4×4 LowerTriangular{Int64, Vector{Int64}}:
  0   0   0  0
  5   0   0  0
  9  10   0  0
 13  14  15  0

Manifold <: AbstractMatrix

A manifold in GeometricOptimizers is a sutype of AbstractMatrix. All manifolds are matrix manifolds and therefore stored as matrices. More details can be found in the docstrings for the StiefelManifold and the GrassmannManifold.

MomentumCache <: OptimizerCache

Cache for the gradient optimizer.

Fields

• x::Manifold: the solution,
• g: the gradient (for the manifold case this is in AbstractLieAlgHorMatrix form),
• δ: the direction,
• Δg: difference in gradients (used in OptimizerStatus),
• g̃: scratch for latest_gradient; see GradientCache for why it is separate from g, and note that this is the cache the separation was measured on,
• g̃_is_current: whether g̃ is the gradient at x; see store_gradient!,
• section: the GlobalSection.

MomentumMethod(α)

The gradient descent algorithm with momentum, i.e. the heavy ball method.

Stores the momentum coefficient α. The momentum is accumulated as

\[ p \gets \alpha{}p + \nabla{}L\]

and the direction is $-p$.

MomentumState <: OptimizerState

State for the gradient optimizer.

NativePade(θ = 0.5) <: AbstractExponentialAlgorithm

Evaluate $\mathfrak{A} = \varphi_1$ with a native degree-6 diagonal Padé approximant.

The argument is first divided by $2^s$ until its induced one-norm is at most θ. On that small argument the method evaluates

\[\mathfrak{A}(X) \approx q_6(X)^{-1}p_6(X)\]

directly at $2n\times{}2n$. The denominator inverse is not a dense solve: five effective Newton–Schulz refinements start from the identity and use only matrix products. Scaling is undone with the same low-rank squaring recursion as ScaledSquaring, so no matrix larger than the input is formed.

The degree and the threshold are paired deliberately, and unlike ScaledSquaring's, 0.5 here is a ceiling and not a preference. The Newton–Schulz count is fixed at five, so the inverse residual is exactly $(\mathbb{I} - q_6)^{32}$, and $\|\mathbb{I} - q_6\|_1 \leq \sum_{k\geq1}|q_k|\theta^k = 0.256$ at $\theta = 1/2$, putting it at 2e-19 — below round-off, where the degree-6 Padé error already sits. Nothing bounds it above that. Worst relative error against AugmentedPade over 400 random $6\times6$ arguments of one-norm exactly $\theta$:

It fails silently: a fixed number of refinements simply stops converging, and nothing in the result says so. ScaledSquaring has no such limit — it sums its series until the terms vanish, which is why its own docstring can sweep $\theta$ over [0.125, 4] and find nothing to choose between — so the two thresholds are not interchangeable and this constructor rejects $\theta > 1/2$. Lowering it is safe and merely adds squarings, which is why it is a parameter at all.

The coefficients are not invented here. $q_6$ is the denominator of the $[7/6]$ Padé approximant of $\exp$, whose closed form is standard [19, 20], and $p_6$ is $(N - D)/x$ for that approximant's numerator $N$ — which divides exactly, both having constant term one, and which inherits the $O(x^{13})$ order. The Newton–Schulz iteration for the inverse is likewise classical [20]. What is assembled here is the pairing: a portable solve-free inverse in place of the dense LU a rational approximant normally needs.

Like ScaledSquaring, this uses GeometricOptimizers.opnorm₁, GeometricOptimizers.unit_matrix, reductions and matrix products only, and therefore runs without scalar indexing on a KernelAbstractions backend. It is an independent portable cross-check rather than the default, and three measurements say why. Its fixed rational evaluation is 0.037 ms against 0.021 ms for ScaledSquaring on the isolated $\mathfrak{A}$ call at $N = 200$, $n = 10$. It allocates 330 KiB there against 201 KiB — 1.6×, and exactly the figure that does not stay a mere constant factor on a backend where an allocation costs a synchronisation. And in Float32 its check is the worst of the three at the top of the norm sweep, 1.0e-4 against 4.0e-5 and 3.4e-5; the forward error is not — 1.2e-5 against 1.1e-5 and 9.6e-6 — so what degrades is the orthogonality of the retracted point, not the agreement with the exponential. AugmentedPade remains the CPU reference that delegates all numerics to Base.exp.

See AbstractExponentialAlgorithm for the comparison.

Newton

Newton's method: the direction solves $\nabla^2f(x)\delta = -\nabla{}f(x)$ with the exact Hessian, which SimpleSolvers.HessianAutodiff supplies.

Unlike BFGS and DFP this needs no approximation to build up, so it converges in few iterations, but it also inherits the Hessian's indefiniteness: where $\nabla^2f$ is not positive definite the direction ascends, and ensure_descent! substitutes the steepest-descent direction for it.

NewtonOptimizerCache <: OptimizerCache

Keys

• x: current iterate (this stores the guess called by the functions generated with linesearch_problem),
• Δx: direction of optimization step (difference between x and x̄); this is obtained by multiplying rhs with the inverse of the Hessian,
• g: gradient value (this stores the gradient associated with x),
• g̃: scratch for latest_gradient — the gradient at a line-search trial point, and at the accepted iterate once solver_step! has refreshed it; see GradientCache for why this is a field of its own and not an alias for g,
• g̃_is_current: whether g̃ is the gradient at x; see store_gradient!,
• Δg: gradient difference (difference between g̃ and g); this is used for computing the OptimizerStatus,
• rhs: the right hand side used to compute the update,
• H: the Hessian matrix evaluated at x,

Also compare this to SimpleSolvers.NonlinearSolverCache.

NewtonOptimizerState <: OptimizerState

The optimizer state is needed to update the Optimizer. This is different from OptimizerStatus and OptimizerResult which serve as diagnostic tools.

We note that this is also used for the BFGS and the DFP optimizer.

Keys

• x
• x̄
• g
• ḡ
• f̄
• f̄

Optimizer

The optimizer that stores all the information needed for an optimization problem.

This problem can be solved by calling solve!(::AbstractVector, ::Optimizer).

Keys

• algorithm::OptimizerState,
• problem::OptimizerProblem,
• gradient::SimpleSolvers.Gradient,
• hessian::SimpleSolvers.Hessian,
• config::SimpleSolvers.Options,
• cache::OptimizerCache,
• linesearch::SimpleSolvers.Linesearch.

Examples

F(x) = sum(sin.(x) .^ 2)
x = ones(3)
algorithm = Newton()
state = OptimizerState(algorithm, x)
optimizer = Optimizer(x, F; algorithm = algorithm, linesearch = Bisection())

solve!(x, state, optimizer)
x

# output

3-element Vector{Float64}:
 0.0
 0.0
 0.0

The choice of line search does not change which minimum this one finds, only how close to it the stopping criterion lets the solve get — Backtracking stops at 1.4e-9 per component, where Bisection above happens to land on 0 exactly:

x = ones(3)
algorithm = Newton()
state = OptimizerState(algorithm, x)
optimizer = Optimizer(x, F; algorithm = algorithm, linesearch = Backtracking())

solve!(x, state, optimizer)
F(x) < 1e-15

# output

true

Optimizer(x, F; ∇F!, mode, algorithm, linesearch, retraction, options_kwargs...)

Build an Optimizer for the objective F at the parameters x.

Implementation

OptimizerCache

See e.g. NewtonOptimizerCache and BFGSCache.

Extended help

OptimizerMethod <: SolverMethod

The OptimizerMethod is used in Optimizer and determines the algorithm that is used.

OptimizerProblem <: AbstractOptimizerProblem

Used in Optimizer. Also compare this to SimpleSolvers.NonlinearProblem.

Examples

julia> x = ones(3); F(x) = sum(sin.(x) .^ 2)
F (generic function with 1 method)

julia> OptimizerProblem(F, x)
OptimizerProblem{Float64, typeof(F), Missing, Missing}(F, missing, missing)

OptimizerResult

Serves as a diagnostic tool for the Optimizer and is the return argument of solve!.

Keys

• status::OptimizerStatus: current status of the optimization,
• x: solution,
• f: function value at solution,
• trace: one OptimizerTraceEntry per iteration if Options.store_trace was set, and empty otherwise. See trace.

An OptimizerState is a data structure that is used to dispatch on different algorithms.

It needs to implement three methods,

initialize!(alg::OptimizerState, ::AbstractVector)
update!(alg::OptimizerState, ::AbstractVector)
solver_step!(::AbstractVector, alg::OptimizerState)

that initialize and update the state of the algorithm and perform an actual optimization step.

Further the following convenience methods should be implemented,

problem(alg::OptimizerState)
gradient(alg::OptimizerState)
hessian(alg::OptimizerState)
linesearch(alg::OptimizerState)

which return the problem to optimize, its gradient and (approximate) Hessian as well as the linesearch algorithm used in conjunction with the optimization algorithm if any.

See NewtonOptimizerState for a struct that was derived from OptimizerState.

OptimizerStatus

Contains residuals (relative and absolute) and various convergence properties.

This is also used in OptimizerResult.

Examples

x = ones(3)
state = NewtonOptimizerState(x)
cache = NewtonOptimizerCache(x)
f = 1.
config = Options()
OptimizerStatus(state, cache, f; config = config)

# output

 * Convergence measures

    |x - x'|               = NaN
    |x - x'|/|x'|          = NaN
    |f(x) - f(x')|         = NaN
    |f(x) - f(x')|/|f(x')| = NaN
    |g(x) - g(x')|         = NaN
    |g(x)|                 = NaN

OptimizerTraceEntry

One iteration's worth of the trace solve! records when SimpleSolvers.Options.store_trace is set.

Keys

• iteration: the iteration number the entry was taken at,
• f: the objective there,
• rg: the gradient residual $\|\nabla{}f\|$ there.

See trace.

ProjectedSkew <: AbstractExponentialAlgorithm

Exponentiate the lift in a basis of its own range, where it is a small skew-symmetric matrix.

$\bar{B}$ is skew-symmetric of rank at most $2n$, so its range and its row space coincide and a thin QR of $B'$ gives an $N\times{}2n$ orthonormal Q with $\bar{B} = QMQ^T$ for $M = Q^T\bar{B}Q$ skew-symmetric and $2n\times{}2n$. Then

\[\exp(\bar{B}) = \mathbb{I} + Q\left(\exp(M) - \mathbb{I}\right)Q^T,\]

and $\exp(M)$ is formed from an eigendecomposition: $iM$ is Hermitian for real skew $M$, so $M = -iV\Lambda{}V^*$ and $\exp(M) = \Re\left(V e^{-i\Lambda} V^*\right)$, which is orthogonal by construction rather than by cancellation.

This is the only algorithm whose check does not degrade with the size of the lift: it stays between 2.1e-15 and 5.3e-15 from $\|\bar{B}\| = 5.8$ to $\|\bar{B}\| = 767$, where the other three drift from 1e-15 to 7e-14. Orthogonality is structural here — it comes from the eigenvector matrix, not from the accuracy of a series. The trade is the forward error against exp(Matrix(B)), which is the largest of the four at all but the very largest lifts and up to 4.4× ScaledSquaring's, and a QR plus an eigendecomposition instead of matrix products, costing 1.1×–1.6× one whole retraction over the sizes measured.

The gap is widest in Float32, where the other three are at the mercy of the format: over the same sweep check climbs to 4.0e-5 for ScaledSquaring, 3.4e-5 for AugmentedPade and 1.0e-4 for NativePade, which of them is worst depending on the lift below the top of that sweep, while this stays between 1.0e-6 and 3.1e-6 from one end to the other. Choose it when staying on the manifold matters more than agreeing with the exponential to the last bit — a long Float32 run, for instance, where check accumulates over thousands of steps.

See AbstractExponentialAlgorithm for the alternatives.

QuasiNewtonOptimizerMethod <: OptimizerMethod

Includes BFGS and DFP.

ScalarMomentAdam(T; β₁, β₂, δ, ambient_norm)

Cayley ADAM, [26, Algorithm 2], as an experimental Stiefel-only baseline for Adam.

The moments are stored in bias-corrected form, as Adam's are, but the second one is a scalar:

\[m_1 \gets \frac{\beta_1 - \beta_1^t}{1 - \beta_1^t}m_1 + \frac{1 - \beta_1}{1 - \beta_1^t}\bar{G},\]

\[m_2 \gets \frac{\beta_2 - \beta_2^t}{1 - \beta_2^t}m_2 + \frac{1 - \beta_2}{1 - \beta_2^t}\lVert\bar{G}\rVert^2,\]

where $\bar{G}\in\mathfrak{g}^\mathrm{hor}$ is the gradient in the global tangent space representation, and the direction is $-m_1/\sqrt{m_2 + \delta}$.

That one scalar is the entire difference from Adam, and it is what the name says: $\lVert\bar{G}\rVert^2$ is a squared gradient norm where Adam's $m_2$ is a squared gradient, so the second moment carries no direction and the method assigns one adaptive learning rate to the whole matrix instead of one per coordinate. That is a reproduction of a published algorithm and not a straw man — on some objectives it will beat Adam.

Only a single StiefelManifold{T} solution is supported; ordinary arrays, NamedTuples, Grassmann solutions and mixed parameter trees throw an ArgumentError. As for Adam, T is the element type of the parameters and is not converted by Optimizer.

Arguments

β₁, β₂ and δ are the source's $\beta_1$, $\beta_2$ and $\varepsilon$, all converted to T, with $0 \le \beta_1, \beta_2 < 1$ and $\delta \ge 0$ validated. ambient_norm selects which $\lVert\cdot\rVert^2$ the second moment accumulates and defaults to false; see GeometricOptimizers._squared_gradient_norm, which is also where the reason the two choices are not interchangeable up to a constant is written down.

Extended help

The source's Algorithm 2, a table mapping each of its symbols onto the code, the derivation that makes its lines 8–10 be global_rep, and the two departures above with their reasons are on the Optimizer Methods page. The third departure is the momentum's transport: the source re-projects $M_k$ onto the tangent space at the new iterate (its equation (6)) where update_section! carries the lift here, which agrees for a Lie group and differs on a proper homogeneous space — the one departure that changes the iterates against the authors' implementation on the same objective and seed.

What is recorded only here is three discrepancies between the source and its own implementation, which are not departures of this port:

• $v_1$. Its line 2 initializes it to $1$, stiefel_optimizer.py to $0$ (v_buffer = torch.zeros([1])), and line 6's $1 - \beta_2^k$ is $0$ at the $k = 0$ its loop starts from. $0$ is followed, which is also Adam's convention here and which makes the first direction a normalized gradient step.
• $\varepsilon$ inside or outside the root. Its line 7 and its implementation (vnew_hat.add(epsilon).sqrt()) both put it inside; Adam here puts it outside, $\sqrt{m_2}+\delta$. The source's placement is used, so this is the one convention ScalarMomentAdam does not share with Adam.
• $s$, the number of fixed-point iterations. Its line 1 gives $s = 2$; Cayley_loop in gutils.py runs for i in range(5). Neither number reaches this method, because the iteration they count is the approximation of the Cayley transform that the admonition above says is not ported — Cayley is evaluated exactly. It is recorded because a reader comparing this port against the authors' code will meet it.

ScalarMomentAdamCache <: OptimizerCache

Cache for ScalarMomentAdam.

The fields are AdamCache's, except that m₂ is a single number rather than an element of 𝔤ʰᵒʳ:

• x::StiefelManifold: the solution,
• g: the gradient, in StiefelLieAlgHorMatrix form,
• δ: the direction,
• Δg: difference in gradients, needed for OptimizerStatus,
• g̃: scratch for latest_gradient; see GradientCache,
• g̃_is_current: whether g̃ is the gradient at x; see store_gradient!,
• m₁: the first moment, a horizontal lift, stored bias-corrected,
• m₂::T: the second moment, the scalar of [26], stored bias-corrected,
• m̃₂::T: √(m₂ + method.δ), the scalar the direction is divided by. The δ in it is ScalarMomentAdam's regularizer and not the δ field above, which is the direction,
• section: the GlobalSection.

ScalarMomentAdamState <: OptimizerState

State for ScalarMomentAdam.

As for ScalarMomentAdamCache, m₂ is a scalar and not a horizontal lift.

ScaledSquaring(θ = 0.5) <: AbstractExponentialAlgorithm

Evaluate $\mathfrak{A}$ by scaling and squaring, and the default.

The series for $\mathfrak{A}$ converges for every argument but is only accurate for a small one: at $\|X\| \gg 1$ its terms cancel catastrophically, and the partial sum reaches $2.5\cdot10^{18}$ where the result is of order one. So halve the argument until it is small, sum the series there, and undo the halving by squaring.

The squaring is done on the $2n\times{}2n$ factor rather than on the assembled $N\times{}N$ matrix, which is possible because the low-rank form is closed under squaring:

\[\left(\mathbb{I} + B'W(B'')^T\right)^2 = \mathbb{I} + B'\left(2W + WXW\right)(B'')^T, \qquad X = (B'')^TB',\]

so one squaring of the exponential is one application of $W \mapsto 2W + WXW$. With $s$ chosen so that $\|X\|_1/2^s \leq θ$, the whole algorithm is s small matrix products on top of a series that now converges in a handful of terms — cheaper than summing the unscaled series, not just more accurate.

θ is the norm below which the series is summed. It barely matters: at $\|\bar{B}\| = 155$ every $θ \in [0.125, 4]$ — a 32-fold range — gives a check between 9.9e-15 and 5.0e-14 and a forward error between 6.4e-15 and 8.2e-15, and neither column is monotone in θ. Nothing in the measurement singles out 0.5; it needs no tuning because no value in that range does better.

Like TaylorSeries and NativePade, this uses nothing but matrix products and norms, so it runs unchanged on a KernelAbstractions backend — including the identities it needs, which come from GeometricOptimizers.unit_matrix and not from Base.one. It remains the default because it is the cheaper of the two usable portable algorithms: the isolated $\mathfrak{A}$ call at $N = 200$, $n = 10$ is 0.021 ms and 201 KiB against NativePade's 0.037 ms and 330 KiB. The norm is taken by GeometricOptimizers.opnorm₁ rather than by LinearAlgebra.opnorm(X, 1), which is a scalar-indexing loop and would give up exactly the property this paragraph claims.

See AbstractExponentialAlgorithm for the alternatives.

SkewSymMatrix(S::AbstractVector, n::Integer)

Instantiate a skew-symmetric matrix with information stored in vector S.

A skew-symmetric matrix $A$ is a matrix $A^T = -A$.

Internally the struct saves a vector $S$ of size $n(n-1)\div2$. The conversion is done the following way:

\[[A]_{ij} = \begin{cases} 0 & \text{if $i=j$} \\ S[( (i-2) (i-1) ) \div 2 + j] & \text{if $i>j$}\\ S[( (j-2) (j-1) ) \div 2 + i] & \text{else}. \end{cases}\]

So $S$ stores a string of vectors taken from $A$: $S = [\tilde{a}_1, \tilde{a}_2, \ldots, \tilde{a}_n]$ with $\tilde{a}_i = [[A]_{i1},[A]_{i2},\ldots,[A]_{i(i-1)}]$.

Also see SymmetricMatrix, LowerTriangular and UpperTriangular.

Examples

using GeometricOptimizers
S = [1, 2, 3, 4, 5, 6]
SkewSymMatrix(S, 4)

# output

4×4 SkewSymMatrix{Int64, Vector{Int64}}:
 0  -1  -2  -4
 1   0  -3  -5
 2   3   0  -6
 4   5   6   0

SkewSymMatrix(A::AbstractMatrix)

Perform 0.5 * (A - A') and store the matrix in an efficient way (as a vector with $n(n-1)/2$ entries).

If the constructor is called with a matrix as input it returns a skew-symmetric matrix via the projection:

\[A \mapsto \frac{1}{2}(A - A^T).\]

Examples

using GeometricOptimizers
M = [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
SkewSymMatrix(M)

# output

4×4 SkewSymMatrix{Float64, Vector{Float64}}:
 0.0  -1.5  -3.0  -4.5
 1.5   0.0  -1.5  -3.0
 3.0   1.5   0.0  -1.5
 4.5   3.0   1.5   0.0

Extended help

Note that the constructor is designed in such a way that it always returns matrices of type SkewSymMatrix{<:AbstractFloat} when called with a matrix, even if this matrix is of type AbstractMatrix{<:Integer}.

If the user wishes to allocate a matrix SkewSymMatrix{<:Integer} then call:

SkewSymMatrix(::AbstractVector, n::Integer)

Note that this is different from LowerTriangular and UpperTriangular as no porjection takes place there.

StiefelLieAlgHorMatrix(A::SkewSymMatrix, B::AbstractMatrix, N::Integer, n::Integer)

Build an instance of StiefelLieAlgHorMatrix based on a skew-symmetric matrix A and an arbitrary matrix B.

An element of StiefelLieAlgMatrix takes the form:

\[\begin{pmatrix} A & B^T \\ B & \mathbb{O} \end{pmatrix},\]

where $A$ is skew-symmetric (this is SkewSymMatrix in GeometricOptimizers).

Also see GrassmannLieAlgHorMatrix.

Extended help

StiefelLieAlgHorMatrix is the horizontal component of the Lie algebra of skew-symmetric matrices (with respect to the canonical metric).

The projection here is: $\pi:S \to SE$ where

\[E = \begin{bmatrix} \mathbb{I}_{n} \\ \mathbb{O}_{(N-n)\times{}n} \end{bmatrix}.\]

The matrix $E$ is implemented under StiefelProjection in GeometricOptimizers.

StiefelLieAlgHorMatrix(D::AbstractMatrix, n::Integer)

Take a big matrix as input and build an instance of StiefelLieAlgHorMatrix.

The integer $N$ in $St(n, N)$ is the number of rows of D.

Extended help

If the constructor is called with a big $N\times{}N$ matrix, then the projection is performed the following way:

\[\begin{pmatrix} A & B_1 \\ B_2 & D \end{pmatrix} \mapsto \begin{pmatrix} \mathrm{skew}(A) & -B_2^T \\ B_2 & \mathbb{O} \end{pmatrix}.\]

The operation $\mathrm{skew}:\mathbb{R}^{n\times{}n}\to\mathcal{S}_\mathrm{skew}(n)$ is the skew-symmetrization operation. This is equivalent to calling of SkewSymMatrix with an $n\times{}n$ matrix.

This can also be seen as the operation:

\[D \mapsto \Omega(E, DE) = \mathrm{skew}\left(2 \left(\mathbb{I} - \frac{1}{2} E E^T \right) DE E^T\right).\]

Also see GeometricOptimizers.Ω.

StiefelManifold <: Manifold

An implementation of the Stiefel manifold [11]. The Stiefel manifold is the collection of all matrices $Y\in\mathbb{R}^{N\times{}n}$ whose columns are orthonormal, i.e.

\[ St(n, N) = \{Y: Y^TY = \mathbb{I}_n \}.\]

The Stiefel manifold can be shown to have manifold structure (as the name suggests) and this is heavily used in GeometricOptimizers. It is further a compact space. More information can be found in the docstrings for rgrad(::StiefelManifold, ::AbstractMatrix) and metric(::StiefelManifold, ::AbstractMatrix, ::AbstractMatrix).

StiefelProjection(backend, T, N, n)

Make a matrix of the form $\begin{bmatrix} \mathbb{I} & \mathbb{O} \end{bmatrix}^T$ for a specific backend and data type.

An array that essentially does vcat(I(n), zeros(N-n, n)) with GPU support.

Extended help

An instance of StiefelProjection should technically also belong to StiefelManifold.

StiefelProjection(B::AbstractLieAlgHorMatrix)

Extract necessary information from B and build an instance of StiefelProjection.

Necessary information here referes to the backend, the data type and the size of the matrix.

The size is queried through B.N and B.n.

Examples

using GeometricOptimizers
using GeometricOptimizers: StiefelProjection

B₁ = rand(StiefelLieAlgHorMatrix, 5, 2)
B₂ = rand(GrassmannLieAlgHorMatrix, 5, 2)
E = [1. 0.; 0. 1.; 0. 0.; 0. 0.; 0. 0.]

StiefelProjection(B₁) ≈ StiefelProjection(B₂) ≈ E 

# output

true

StiefelProjection(A::AbstractMatrix)

Extract necessary information from A and build an instance of StiefelProjection.

Necessary information here referes to the backend, the data type and the size of the matrix.

SymmetricMatrix(S::AbstractVector, n::Integer)

Instantiate a symmetric matrix with information stored in vector S.

A SymmetricMatrix $A$ is a matrix $A^T = A$.

Internally the struct saves a vector $S$ of size $n(n+1)\div2$. The conversion is done the following way:

\[[A]_{ij} = \begin{cases} S[( (i-1) i ) \div 2 + j] & \text{if $i\geq{}j$}\\ S[( (j-1) j ) \div 2 + i] & \text{else}. \end{cases}\]

So $S$ stores a string of vectors taken from $A$: $S = [\tilde{a}_1, \tilde{a}_2, \ldots, \tilde{a}_n]$ with $\tilde{a}_i = [[A]_{i1},[A]_{i2},\ldots,[A]_{ii}]$.

Also see SkewSymMatrix, LowerTriangular and UpperTriangular.

Examples

using GeometricOptimizers
S = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
SymmetricMatrix(S, 4)

# output

4×4 SymmetricMatrix{Int64, Vector{Int64}}:
 1  2  4   7
 2  3  5   8
 4  5  6   9
 7  8  9  10

SymmetricMatrix(A::AbstractMatrix)

Perform a projection and store the matrix in an efficient way (as a vector with $n(n+1)/2$ entries).

If the constructor is called with a matrix as input it returns a symmetric matrix via the projection:

\[A \mapsto \frac{1}{2}(A + A^T).\]

Examples

using GeometricOptimizers
M = [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
SymmetricMatrix(M)

# output

4×4 SymmetricMatrix{Float64, Vector{Float64}}:
 1.0   3.5   6.0   8.5
 3.5   6.0   8.5  11.0
 6.0   8.5  11.0  13.5
 8.5  11.0  13.5  16.0

Extended help

Note that the constructor is designed in such a way that it always returns matrices of type SymmetricMatrix{<:AbstractFloat} when called with a matrix, even if this matrix is of type AbstractMatrix{<:Integer}.

If the user wishes to allocate a matrix SymmetricMatrix{<:Integer} then call

$julia SymmetricMatrix(::AbstractVector, n::Integer)$`

Note that this is different from LowerTriangular and UpperTriangular as no porjection takes place there.

TaylorSeries <: AbstractExponentialAlgorithm

Sum the series for $\mathfrak{A}$ directly, without scaling. This is not a usable retraction.

It is the behaviour of every version of this package up to 0.2.0, retained only so that the regression is reproducible from the test suite and so the working algorithms have a baseline to be compared against. check(geodesic(B, TaylorSeries())), on a random StiefelLieAlgHorMatrix(20, 3) scaled up:

At $\|\bar{B}\| = 79$ the "retracted" point is not on the Stiefel manifold in any sense, and by $767$ the series has overflowed. The series is summed on $X = (B'')^TB'$ whose norm is $\approx\|\bar{B}\|^2/4$, and its terms cancel: the partial sum reaches $2.5\cdot10^{18}$ where the result is of order one, so stopping when a term falls below eps leaves a relative error of $\varepsilon\|\mathfrak{A}(X)\|$ rather than $\varepsilon$. Making the termination test relative to the partial sum instead of absolute does not help, and was measured not to change any of the numbers above — the loss is in the cancellation, not in when the summation stops.

Use ScaledSquaring, which fixes this and is also the cheaper of the two — 1.7× at $N = 200$, $n = 10$ and 4.6× at $N = 500$, $n = 50$, because the scaled series converges in a handful of terms where the unscaled one grinds through hundreds.

See AbstractExponentialAlgorithm for the alternatives.

UpperTriangular(S::AbstractVector, n::Int)

Build an upper-triangular matrix from a vector.

An upper-triangular matrix is an $n\times{}n$ matrix that has zeros on the diagonal and on the lower triangular.

The data are stored in a vector $S$ similarly to other matrices. See LowerTriangular, SkewSymMatrix and SymmetricMatrix.

The struct two fields: S and n. The first stores all the entries of the matrix in a sparse fashion (in a vector) and the second is the dimension $n$ for $A\in\mathbb{R}^{n\times{}n}$.

Examples

using GeometricOptimizers
S = [1, 2, 3, 4, 5, 6]
UpperTriangular(S, 4)

# output

4×4 UpperTriangular{Int64, Vector{Int64}}:
 0  1  2  4
 0  0  3  5
 0  0  0  6
 0  0  0  0

UpperTriangular(A::AbstractMatrix)

Build an upper-triangular matrix from a matrix.

This is done by taking the upper right of that matrix.

Examples

using GeometricOptimizers
M = [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
UpperTriangular(M)

# output

4×4 UpperTriangular{Int64, Vector{Int64}}:
 0  2  3   4
 0  0  7   8
 0  0  0  12
 0  0  0   0

λY * Y

Apply the element λY onto Y.

Here λY is an element of a Lie group and Y is an element of a homogeneous space.

one(B::AbstractLieAlgHorMatrix)

The $N\times{}N$ identity, built with a KernelAbstractions kernel.

Base.one(::AbstractMatrix) writes the diagonal in a scalar-indexed loop, which is what a GPU array cannot serve; geodesic reaches this on every retraction. It existed for the Stiefel lift only, so the Grassmann retraction was taking the scalar-indexed path — the same hazard issue A19 recorded for GeometricOptimizers.𝔄, whose argument is a bare matrix and which reached Base.one until the $2n\times{}2n$ identities went through GeometricOptimizers.unit_matrix as well.

parent(B::AbstractLieAlgHorMatrix)

The tuple of blocks B's free parameters are stored in — (A, B) for a StiefelLieAlgHorMatrix, (B,) for a GrassmannLieAlgHorMatrix — and not the single array every other parent this package defines returns.

Every operation on a lift that is elementwise in the free parameters — as opposed to in the ambient $N\times{}N$ matrix, which has no setindex! and counts each off-diagonal block twice — is written once over this tuple rather than once per lift type. That is the four methods below, and l2norm, and the _difference! / _add! / _rac! / _div! / _square! family in named_tuple_wrapper.jl. They used to exist for the Stiefel lift alone, which is half of why a GrassmannManifold could not be driven through an Optimizer at all (issue A11).

The docstring is attached to the signature and not to the bare Base.parent: the package also defines parent for Manifold, SkewSymMatrix, SymmetricMatrix and AbstractTriangular, each of which returns the single array it wraps, and a signature-less docstring would be shown as the general meaning of parent for all of them.

rand(backend, manifold_type, N, n)

Draw random elements for a specific device.

Examples

Random elements of the manifold can be allocated on GPU. Call ...

rand(CUDABackend(), StiefelManifold{Float32}, N, n)

... for drawing elements on a CUDA device.

rand(manifold_type, N, n)

Draw random elements from the Stiefel and the Grassmann manifold.

Because both of these manifolds are compact spaces we can sample them uniformly [30].

Examples

When we call ...

using GeometricOptimizers
using GeometricOptimizers: _round # hide
import Random
Random.seed!(123)

N, n = 5, 3
Y = rand(StiefelManifold{Float32}, N, n)
_round(Y; digits = 5) # hide

# output

5×3 StiefelManifold{Float32, Matrix{Float32}}:
 -0.27575   0.32991   0.77275
 -0.62485  -0.33224  -0.0686
 -0.69333   0.36724  -0.18988
 -0.09295  -0.73145   0.46064
  0.2102    0.33301   0.38717

... the sampling is done by first allocating a random matrix of size $N\times{}n$ via Y = randn(Float32, N, n).

We then perform a QR decomposition Q, R = qr(Y) with the qr function from the LinearAlgebra package (this is using Householder reflections internally).

The final output are then the first n columns of the Q matrix.

vec(B::AbstractLieAlgHorMatrix)

The free parameters of B, laid out end to end and lazily — not the $N^2$ entries of the matrix B presents itself as.

Examples

using GeometricOptimizers

A = SkewSymMatrix([1, ], 2)
B = [2 3; ]
B̄ = StiefelLieAlgHorMatrix(A, B, 3, 2)
B̄ |> vec

# output

vcat(1-element Vector{Int64}, 2-element Vector{Int64}):
 1
 2
 3

Implementation

This is using Vcat from the package LazyArrays, so nothing is copied.

vec(A::AbstractTriangular)

Return the associated vector to $A$.

Examples

using GeometricOptimizers

M = [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
LowerTriangular(M) |> vec

# output

6-element Vector{Int64}:
  5
  9
 10
 13
 14
 15

vec(A)

Output the associated vector of A.

Examples

using GeometricOptimizers

M = [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
SkewSymMatrix(M) |> vec

# output

6-element Vector{Float64}:
 1.5
 3.0
 1.5
 4.5
 3.0
 1.5

update!(cache::NewtonOptimizerCache, x, g, hes)

Update an instance of NewtonOptimizerCache based on x.

This is used in update!(::OptimizerState, ::AbstractVector).

This sets:

\[\begin{aligned} % \bar{x}^\mathtt{cache} & \gets x, \\ x^\mathtt{cache} & \gets x, \\ g^\mathtt{cache} & \gets g, \\ \mathrm{rhs}^\mathtt{cache} & \gets -g, \\ H^\mathtt{cache} & \gets H(x), \\ \delta^\mathtt{cache} & \gets (H^\mathtt{cache})^{-1}\mathrm{rhs}^\mathtt{cache}, \end{aligned}\]

where we wrote $H$ for the Hessian (i.e. the input argument hes).

update!(state::NewtonOptimizerState, gradient, x)

Update an instance of NewtonOptimizerState based on x and gradient, where g is of type SimpleSolvers.Gradient.

Examples

If we only call update! once there are still NaNs for x̄, ḡ and f̄.

f(x) = sum(x.^2)
x = [1., 2.]
state = NewtonOptimizerState(x)
grad = GradientAutodiff{Float64}(f, length(x))
update!(state, grad, x)

# output

NewtonOptimizerState{Float64, Vector{Float64}, Vector{Float64}, GlobalSection{Float64, Vector{Float64}, Nothing}}(0, [1.0, 2.0], [NaN, NaN], [2.0, 4.0], [NaN, NaN], 5.0, NaN, GlobalSection{Float64, Vector{Float64}, Nothing}([1.0, 2.0], nothing))

update!(cache, x, g)

Update the BFGSCache based on x and g.

Extended help

The update rule used here can be found in [31] and [27]:

It does:

\[\begin{aligned} \delta & \gets x^{(k)} - x^{(k-1)}, \\ \gamma & \gets \nabla{}f^{(k)} - \nabla{}f^{(k-1)}, \\ T_1 & \gets \delta\gamma^TQ, \\ T_2 & \gets Q\gamma\delta^T, \\ T_3 & \gets (1 + \frac{\gamma^TQ\gamma}{\delta^T\gamma})\delta\delta^T,\\ Q & \gets Q - (T_1 + T_2 - T_3)/{\delta^T\gamma} \end{aligned}\]

update!(cache, x, g)

Update the DFPCache based on x and g.

Extended help

The update rule used here can be found in [31] and [27].

AdamOptimizerWithDecay(n_epochs, T; η₁, η₂, kwargs...)

Adam paired with a DecayingStatic line search whose learning rate decays geometrically from η₁ to η₂ over n_epochs, returned as a NamedTuple to splat into Optimizer:

method = AdamOptimizerWithDecay(1000)
opt    = Optimizer(x, problem; method...)
solve!(x, OptimizerState(method.algorithm, x), opt)

The state has to be built from method.algorithm, so the pairing is worth binding to a name rather than splatting it twice.

There is no new type here and no schedule of its own: this is exactly Adam(T) and DecayingStatic(T; η₁, η₂, n = n_epochs), under the one name that the two together used to have.

Why it exists

GeometricMachineLearning carried an AdamOptimizerWithDecay that bundled Adam's ρ₁, ρ₂, δ with a learning-rate schedule η₁, η₂, n_epochs, and computed the same $\gamma = \exp(\log(\eta_2/\eta_1)/n)$. Since 0.2.0 the two halves of that live in different places here — the direction in an OptimizerMethod, the step size in a SimpleSolvers.LinesearchMethod — which is the right split but left no single name to migrate that method to. This is it, and GML deleted its own copy in 0.5, where both packages exporting the name made using the two of them together an error; see GeometricMachineLearning#239.

Arguments

η₁, η₂ and n_epochs go to the line search; everything else is forwarded to Adam, so β₁, β₂ and δ — GML's ρ₁, ρ₂ and δ — keep Adam's own defaults rather than a second copy of them. T is the element type of the parameters and is positional, as it is for Adam and DecayingStatic.

Examples

AdamOptimizerWithDecay(1000).linesearch

# output

DecayingStatic from α = 0.01 to α = 1.0e-6 over 1000 iterations.

AdamW(args...)

Deliberately undefined; use AdamWithEuclideanDecay.

The decoupled weight decay of AdamW is the Euclidean $λx$, whose Riemannian gradient vanishes identically on every manifold of this package, so a method called AdamW here would be Adam under a second name for anyone optimizing on a manifold — and would say so nowhere. Rather than let that be discovered at run time, the name errors and points at the one that describes what it does. It is kept free for a Riemannian weight decay, should issue #28 conclude that one is called for.

_div!(C, A, B)

Divide A by B (elment-wise)

_dot(a, b)

The inner product of two gradients or directions, taken in the flattened coordinates.

Implementation

For an AbstractVecOrMat this is LinearAlgebra.dot. For a horizontal lift — or a NamedTuple of them — it is emphatically not: dot on an AbstractLieAlgHorMatrix is the ambient Frobenius product, which counts each of the off-diagonal blocks of the lift twice and so comes out exactly twice the product of the free parameters. The intrinsic coordinates are the ones every other quantity in this package is expressed in — Q is sized by the flattening, outer! flattens before it forms its outer product, and the α of a line search parameterizes a curve in them — so pairing a gradient with a direction has to happen there too.

Used by trial_slope for $\varphi'(\alpha)$, and by the quasi-Newton caches for $\delta^T\gamma$, whose value has to be consistent with the flattened T₁, T₂ and γ^TQγ it divides.

_is_decayable(x)

Whether Euclidean weight decay can move x at all.

An ordinary array is decayable and a NamedTuple is decayable as soon as one of its entries is, which is the case AdamWithEuclideanDecay exists for. The StiefelManifold and the GrassmannManifold are not: $\lambda{}x$ is the gradient of $\frac{\lambda}{2}||x||^2$, which is constant on them, so the decay vanishes there identically (see the Weight Decay on Manifolds page).

Implementation

This is the only place that geometric fact is recorded — _weight_decay! consults it rather than restating it, so the no-op and the warning in OptimizerCache cannot drift apart.

It is declared on the two concrete manifolds and not on Manifold on purpose. What makes the decay vanish is that these two are compact, with $||Y||_F^2 = \mathrm{tr}(Y^TY) = n$; it is not a consequence of being a manifold, and on a noncompact one Euclidean decay may well do something. A manifold added later therefore has to decide for itself, and gets an error saying so rather than inheriting a silent no-op that may be wrong for it.

_manifold_αmax(solution_blocks, direction_blocks, c)

The ceiling a step ceiling of c imposes on a solution made of several blocks: the smallest of the per-block step_αmax over the blocks that live on a Manifold, and Inf where none does.

One \alpha is applied to every block, so each manifold block needs $\|\alpha\delta_i\| \leq 2\pi{}c$ and the binding one is the largest $\|\delta_i\|$. A block that is an ordinary array contributes nothing: the $2\pi$ is the turn of a rotation and a vector space has no such scale, so it must neither impose a ceiling of its own nor inflate the norm that sets one for its neighbours.

That last part is what this function exists for. The ceiling used to be step_αmax(c, direction(cache)) over the whole direction, whose l2norm combines the blocks in quadrature — so a NamedTuple of ordinary arrays was bounded by a rotation that does not exist in its problem (measured: a Euclidean NamedTuple solve took 3 184 iterations against 1 for the same problem written as a vector), and in a mixed NamedTuple the Euclidean blocks tightened the manifold blocks' bound for no reason (measured on the mixed problem of test/named_tuple_parameters.jl: $\|\delta_Y\| = 2.5\times10^{-16}$ against a total of 3.9, bounding $\alpha$ at 1.6 where the geometry of the manifold block permits $2.6\times10^{16}$). That was catalogued as issue A15 and is what this closes.

Implementation

Recursive over the two tuples rather than a loop over zip, so that a heterogeneous NamedTuple — a StiefelManifold beside a Matrix beside a Vector — stays inferable and the returned parameters stay concrete.

_optimizer(x, problem, algorithm, linesearch, gradient, retraction, config)

Build the cache and the Hessian for algorithm and hand everything to Optimizer's inner constructor.

Takes every argument positionally on purpose; see the note on Julia 1.12 below Optimizer(x, F).

_rac!(B, A)

Compute the element-wise square-root of A.

_square!(B, A)

_squared_gradient_norm(method::ScalarMomentAdam, cache, gradient, x)

The ‖·‖² that ScalarMomentAdam's second moment accumulates, selected by method.ambient_norm.

false — the default — squares the horizontal lift already in the cache; true squares [26]'s own ambient Euclidean gradient, at one extra gradient evaluation per step. Gradient applied to a Manifold returns rgrad(Y, ∇L) and the ambient ∇L is not part of the optimizer protocol, so true recovers it from the flattened closure GradientAutodiff(F, ::Manifold) builds — which is what gradient(vec(x)) is. store_gradient!'s reuse does not help, because what it caches is the lift.

The two are not interchangeable up to a constant, and which to reach for when is on the Optimizer Methods page.

_weight_decay!(δ, x, λ)

Subtract the decoupled weight decay $\lambda{}x$ from the direction δ, in place.

This is called after the Adam direction has been formed, i.e. δ is $-m_1/(\sqrt{m_2} + \delta)$ on entry and $-m_1/(\sqrt{m_2} + \delta) - \lambda{}x$ on exit. It has to happen before the line search scales the direction, so that the decay is multiplied by the learning rate as [28] prescribes, and it must not touch the moments — that is the whole point of the decoupling, see AdamWithEuclideanDecay.

Implementation

The Manifold method is a genuine no-op and not an omission: δ is then an element of $\mathfrak{g}^\mathrm{hor}$ rather than something that could be added to x at all, and the quantity it would have to be corrected by — the Riemannian gradient of $\frac{\lambda}{2}||x||^2$ — is zero. It asserts _is_decayable rather than restating why, so that the no-op holds only for the manifolds on which that has been established; see issue #28 for what a decay that did something here would have to look like. test/adam_with_euclidean_decay.jl checks that rgrad(Y, λY) really does vanish rather than taking it on trust.

apply_section!(Y::AT, λY::GlobalSection{T, AT}, Y₂::AT) where {T, AT<:StiefelManifold{T}}

Apply λY to Y₂ and store the result in Y.

This is the inplace version of apply_section.

apply_section(λY::GlobalSection{T, AT}, Y₂::AT) where {T, AT <: StiefelManifold{T}}

Apply λY to Y₂.

Mathematically this is the group action of the element $\lambda{}Y\in{}G$ on the element $Y_2$ of the homogeneous space $\mathcal{M}$.

Internally it calls apply_section!.

cayley(B̄::GrassmannLieAlgHorMatrix)

Compute the Cayley retraction of B.

This is equivalent to the method of cayley for StiefelLieAlgHorMatrix.

See cayley(::StiefelLieAlgHorMatrix).

cayley(B̄::StiefelLieAlgHorMatrix)

Compute the Cayley retraction of B.

Implementation

Internally this is using

\[\mathrm{Cayley}(\bar{B}) = \mathbb{I} + \frac{1}{2} B' (\mathbb{I}_{2n} - \frac{1}{2} (B'')^T B')^{-1} (B'')^T (\mathbb{I} + \frac{1}{2} B),\]

with

\[\bar{B} = \begin{bmatrix} A & -B^T \\ B & \mathbb{O} \end{bmatrix} = \begin{bmatrix} \frac{1}{2}A & \mathbb{I} \\ B & \mathbb{O} \end{bmatrix} \begin{bmatrix} \mathbb{I} & \mathbb{O} \\ \frac{1}{2}A & -B^T \end{bmatrix} =: B'(B'')^T,\]

i.e. $\bar{B}$ is expressed as a product of two $N\times{}2n$ matrices.

cayley(Y::Manifold, Δ)

Take as input an element of a manifold Y and a tangent vector in Δ in the corresponding tangent space and compute the Cayley retraction.

In different notation: take as input an element $x$ of $\mathcal{M}$ and an element of $T_x\mathcal{M}$ and return $\mathrm{Cayley}(v_x).$

Examples

using GeometricOptimizers

Y = StiefelManifold([1. 0. 0.;]' |> Matrix)
Δ = [0. .5 0.;]' |> Matrix
Y₂ = GeometricOptimizers.cayley(Y, Δ)

Y₂' * Y₂ ≈ [1.;]

# output

true

See the example in geodesic(::Manifold, ::AbstractMatrix).

check(Y::Manifold)

Measure how far Y is from the manifold, as $\|Y^TY - \mathbb{I}\|$.

Both manifolds this package provides store a representative whose columns are orthonormal — for StiefelManifold that is the point itself, for GrassmannManifold it is the representative of the equivalence class — so the same expression measures both. A retraction maps onto the manifold by construction, so in exact arithmetic this is zero and what it actually returns is accumulated round-off.

This is the assertion the manifold tests rest on. It used to exist for StiefelManifold only, which is why the accuracy loss in GeometricOptimizers.𝔄 went unnoticed for so long: half the retraction paths had nothing that could have caught it.

Examples

using GeometricOptimizers
using GeometricOptimizers: check
import Random
Random.seed!(123)

check(rand(GrassmannManifold, 5, 3)) < 1e-14

# output

true

contains_nonfinite(a)

Whether a holds any value that is not finite.

This was contains_nan, and tested isnan only. NaN is the last thing a diverging solve produces: it reaches Inf first, and before that every finite magnitude on the way. On the SVD problem of test/optimizer_convergence/svd_optim.jl the diverging solve passed through f = 1.2e169 and check(Y) = 1.07e200 — both perfectly ordinary Float64s, neither of them NaN — and only went NaN on the iteration after that. By then it had been off the manifold for two iterations.

isfinite still does not catch 1e200, which is why it is not the only guard; see convergence_measures for the one that does.

convergence_measures(status, config)

Checks if the optimizer converged.

Here status is an OptimizerStatus object and config is an SimpleSolvers.Options object.

Extended help

curvature_is_usable(ΔxΔg, Δx, Δg)

Whether the secant pair $(\delta, \gamma)$ satisfies the curvature condition, i.e. whether the quasi-Newton update built from it keeps $Q$ positive definite.

Both the BFGS and the DFP update divide by $\delta^T\gamma$ and preserve positive definiteness of $Q$ only for $\delta^T\gamma > 0$ [27]. The guard this replaced was !iszero(ΔxΔg) && !isnan(ΔxΔg), which admits both signs and, more importantly, admits denominators that are zero to within round-off: on the SVD problem of test/optimizer_convergence/svd_optim.jl, $\delta^T\gamma$ took the values -12.8, -4.5e-16 and +1.5e-15 on consecutive iterations, all three of which !iszero accepts. Dividing a rank-two correction by 1.5e-15 is what drove $\lambda_\mathrm{max}(Q)$ from 3 to 442 there, and $\lambda_\mathrm{min}(Q)$ to -398 from another starting point.

The threshold is therefore relative – an absolute one cannot tell 1.5e-15 on a problem scaled to $10^0$ from a legitimate small pairing on a problem scaled to $10^{-15}$. Its exact value matters much less than its existence: 1e-8 and eps(T) behave identically on everything measured here, because the pairs being rejected are the ones that are non-positive rather than merely small.

default_gradient(problem, x)

Return the SimpleSolvers.Gradient that Optimizer uses if none is supplied.

Implementation

The NamedTuple method is not just a matter of the length: a Gradient built for a NamedTuple is called on the flattened parameters, so it has to be constructed from x itself (see GradientAutodiff(F, ::NamedTuple)), which composes problem.F with the unflatten that belongs to x. Sizing it with length(x) — the number of entries of the NamedTuple rather than the length of its flattening — used to make the first step fail with a DimensionMismatch.

default_linesearch(T, method)

Return the line search that Optimizer uses for method if none is supplied.

Everything except Adam defaults to SimpleSolvers.Backtracking(T; expand = true): the (quasi-)Newton methods because they build a direction with a scale of its own, and GradientMethod and MomentumMethod because a searching line search is what makes them converge rather than merely descend. All produce genuine descent directions, so a backtracking search always has an α to find.

expand = true is what lets that search lengthen a step as well as shorten it, and it is not the SimpleSolvers default — see the tip below for why it is the default here.

The AdamFamily methods are the exception and keep a fixed Static(DEFAULT_LEARNING_RATE). Adam's direction is $-m_1/(\sqrt{m_2} + \delta)$ and ScalarMomentAdam's is $-m_1/\sqrt{m_2 + \delta}$; either way it is a moving average that is deliberately not required to descend on any individual step, so a sufficient-decrease search has nothing to work with and would spend every such step reporting that it found no descent direction.

One caveat on sharing the value DEFAULT_LEARNING_RATE across the family: Adam normalizes componentwise and ScalarMomentAdam normalizes the whole lift at once, so the same $\eta$ is a step $\sqrt{\dim}$ shorter for the latter. See the warning in its docstring.

For AdamWithEuclideanDecay a fixed step is not merely the cheaper choice but the only one under which $\lambda$ means what it is documented to mean. Two reasons, and the first is the one that matters:

1. The merit the line search minimizes is the bare objective $f$. It is not $f + \frac{\lambda}{2}||x||^2$, because this package regularizes by adding $-\lambda{}x$ to the direction rather than by penalizing the objective — that is what decoupled means. So the penalty is invisible to the search, which will spend its $\alpha$ undoing the decay's contribution to $f$ as far as it can. The search is not approximating a line minimum of the function the method is actually descending.
2. The decay per step is $\alpha\lambda$, so a varying $\alpha$ makes the effective regularization strength whatever the search happened to pick on that step.

DecayingStatic is the exception to the exception: it varies $\alpha$ too, but on a schedule rather than in response to the merit, which is exactly AdamW's own learning-rate schedule and leaves $\lambda$ its meaning relative to $\eta$.

ensure_descent!(cache, method, config)

Replace the direction stored in cache by the steepest-descent direction if it does not descend.

A (quasi-)Newton step solves $H\delta = -\nabla{}f$ and is a descent direction only where $H$ is positive definite. Where it is not — an indefinite Hessian near a saddle, or a quasi-Newton $Q$ that has lost positive definiteness — $\delta$ points uphill. Since a line search only ever returns a non-negative step, the iteration then walks to the nearest stationary point, which may be a maximum.

F(x) = sum(sin.(x) .^ 2) started from x = ones(3) is the smallest example: $\sin^2$ has second derivative $2\cos(2x) = -0.83$ there, so the Newton direction ascends and the solve converges to $\pi/2$, where F is maximal. With this safeguard it descends to 0.

Descent is tested as $\nabla{}f\cdot\delta < 0$. rhs stores $-\nabla{}f$, so the test is $\mathrm{rhs}\cdot\delta > 0$; a direction that fails it — including a NaN one, for which every comparison is false — is replaced by rhs, which always descends. The substitution is reported at config.verbosity ≥ 2.

Examples

F(x) = sum(sin.(x) .^ 2)
x = ones(3)
algorithm = Newton()
state = OptimizerState(algorithm, x)
optimizer = Optimizer(x, F; algorithm = algorithm, linesearch = Bisection())

solve!(x, state, optimizer)
F(x) < 1e-20

# output

true

geodesic(B̄::AbstractLieAlgHorMatrix, algorithm = ScaledSquaring())

Compute the geodesic of a horizontal lift, i.e. $\exp(\bar{B})$.

Works for both StiefelLieAlgHorMatrix and GrassmannLieAlgHorMatrix — manifold_type supplies the manifold the result belongs to.

Implementation

Internally this is using

\[\exp(\bar{B}) = \mathbb{I} + B'\mathfrak{A}(B', B'')(B'')^T,\]

with $\bar{B} = B'(B'')^T$ the factorisation of lift_factors. The only matrix function this evaluates is $\mathfrak{A}$ on the $2n\times{}2n$ product $(B'')^TB'$, so the cost is set by $n$ and not by $N$, and so is the accuracy — see AbstractExponentialAlgorithm for the choice of algorithm and GeometricOptimizers.𝔄 for the implementations.

geodesic(Y::Manifold, Δ, algorithm = ScaledSquaring())

Take as input an element of a manifold Y and a tangent vector in Δ in the corresponding tangent space and compute the geodesic (exponential map).

In different notation: take as input an element $x$ of $\mathcal{M}$ and an element of $T_x\mathcal{M}$ and return $\mathtt{geodesic}(x, v_x) = \exp(v_x).$

algorithm selects how the exponential is evaluated and is passed straight through to geodesic(::AbstractLieAlgHorMatrix); see AbstractExponentialAlgorithm for the choice. It matters only for a large $\Delta$ — the default ScaledSquaring is accurate at every step size.

Examples

using GeometricOptimizers

Y = StiefelManifold([1. 0. 0.;]' |> Matrix)
Δ = [0. .5 0.;]' |> Matrix
Y₂ = GeometricOptimizers.geodesic(Y, Δ)

Y₂' * Y₂ ≈ [1.;]

# output

true

Implementation

Internally this geodesic method calls geodesic(::AbstractLieAlgHorMatrix).

global_rep(λY::GlobalSection{T, AT}, Δ::AbstractMatrix{T}) where {T, AT<:GrassmannManifold{T}}

Express Δ (an element of the tangent space of Y) as an instance of GrassmannLieAlgHorMatrix.

The method global_rep for GrassmannManifold is similar to that for StiefelManifold.

Examples

using GeometricOptimizers
using GeometricOptimizers: _round
import Random

Random.seed!(123)

Y = rand(GrassmannManifold, 6, 3)
Δ = rgrad(Y, randn(6, 3))
λY = GlobalSection(Y)

_round(global_rep(λY, Δ); digits = 3)

# output

6×6 GrassmannLieAlgHorMatrix{Float64, Matrix{Float64}}:
  0.0     0.0     0.0     0.981  -2.058   0.4
  0.0     0.0     0.0    -0.424   0.733  -0.919
  0.0     0.0     0.0    -1.815   1.409   1.085
 -0.981   0.424   1.815   0.0     0.0     0.0
  2.058  -0.733  -1.409   0.0     0.0     0.0
 -0.4     0.919  -1.085   0.0     0.0     0.0

global_rep(λY::GlobalSection{T, AT}, Δ::AbstractMatrix{T}) where {T, AT<:StiefelManifold{T}}

Express Δ (an the tangent space of Y) as an instance of StiefelLieAlgHorMatrix.

This maps an element from $T_Y\mathcal{M}$ to an element of $\mathfrak{g}^\mathrm{hor}$.

These two spaces are isomorphic where the isomorphism where the isomorphism is established through $\lambda(Y)\in{}G$ via:

\[T_Y\mathcal{M} \to \mathfrak{g}^{\mathrm{hor}}, \Delta \mapsto \lambda(Y)^{-1}\Omega(Y, \Delta)\lambda(Y).\]

Also see GeometricOptimizers.Ω.

Examples

using GeometricOptimizers
using GeometricOptimizers: _round
import Random

Random.seed!(123)

Y = rand(StiefelManifold, 6, 3)
Δ = rgrad(Y, randn(6, 3))
λY = GlobalSection(Y)

_round(global_rep(λY, Δ); digits = 3)

# output

6×6 StiefelLieAlgHorMatrix{Float64, SkewSymMatrix{Float64, Vector{Float64}}, Matrix{Float64}}:
  0.0     0.679   1.925   0.981  -2.058   0.4
 -0.679   0.0     0.298  -0.424   0.733  -0.919
 -1.925  -0.298   0.0    -1.815   1.409   1.085
 -0.981   0.424   1.815   0.0     0.0     0.0
  2.058  -0.733  -1.409   0.0     0.0     0.0
 -0.4     0.919  -1.085   0.0     0.0     0.0

Implementation

The function global_rep does in fact not perform the entire map $\lambda(Y)^{-1}\Omega(Y, \Delta)\lambda(Y)$ but only

\[\Delta \mapsto \mathrm{skew}(Y^T\Delta),\]

to get the small skew-symmetric matrix $A\in\mathcal{S}_\mathrm{skew}(n)$ and

\[\Delta \mapsto (\lambda(Y)_{[1:N, n:N]}^T \Delta)_{[1:(N-n), 1:n]},\]

to get the arbitrary matrix $B\in\mathbb{R}^{(N-n)\times{}n}$.

global_section(Y::GrassmannManifold)

Compute a matrix of size $N\times(N-n)$ whose columns are orthogonal to the columns in Y.

The method global_section for the Grassmann manifold is equivalent to that for the StiefelManifold (we represent the Grassmann manifold as an embedding in the Stiefel manifold).

See the documentation for global_section(Y::StiefelManifold{T}) where T.

global_section(Y::StiefelManifold)

Compute a matrix of size $N\times(N-n)$ whose columns are orthogonal to the columns in Y.

This matrix is also called $Y_\perp$ [6, 8, 9].

Examples

using GeometricOptimizers
using GeometricOptimizers: global_section
import Random

Random.seed!(123)

Y = StiefelManifold([1. 0.; 0. 1.; 0. 0.; 0. 0.])

round.(Matrix(global_section(Y)); digits = 3)

# output

4×2 Matrix{Float64}:
 0.0    -0.0
 0.0     0.0
 0.936  -0.353
 0.353   0.936

Further note that we convert the QRCompactWYQ object to a Matrix before we display it.

Implementation

The implementation is done with a QR decomposition (LinearAlgebra.qr!). Internally we do:

A = randn(N, N - n) # or the gpu equivalent
A = A - Y.A * (Y.A' * A)
qr!(A).Q

gradient(cache)

Return the stored gradient (array) of an instance of BFGSCache

gradient(cache)

Return the stored gradient (array) of an instance of DFPCache

gradient(::NewtonOptimizerCache)

Return the stored gradient (array) of an instance of NewtonOptimizerCache

gradient_difference!(cache, state)

Write the gradient difference ∇f(xᵏ) - ∇f(xᵏ⁻¹) into cache.Δg, for OptimizerStatus.

Implementation

Every cache in this package overrides this with the difference of the two gradients it holds itself, latest_gradient at $x_{k+1}$ and gradient at $x_k$, which is the successive difference the status prints and needs no state.ḡ. The default – cache.g against state.ḡ – was wrong for each of them in a different way:

• for the three first-order caches, update!(::MomentumState, ...) runs after the step and copies the cache's pre-step gradient into state.g, shifting the one before it into state.ḡ, so state.ḡ ends up two iterates behind cache.g rather than one: on $f(x) = \sum(x^2 + 0.1x^4)$ from [1.5, -0.8, 0.4] with MomentumMethod + Bisection, iteration three reported rgₐ = 4.976 where $\|\nabla{}f(x_k) - \nabla{}f(x_{k-1})\| = 0.295$. On the first iteration it differenced against the _similar memory MomentumState never writes.
• for BFGSCache and DFPCache it was the γ of the secant pair, $\nabla{}f(x_k) - \nabla{}f(x_{k-1})$, which those form inside update!(cache, ...) and which is one step behind the rg reported next to it.
• for NewtonOptimizerCache it was structurally zero: solver_step! advances state.ḡ at the same iterate the cache takes its gradient at, so the difference could only ever be 0.

In all three cases the two g rows of a status are now about one step rather than about two different ones; see convergence_measures for which iterate rg belongs to.

state.ḡ is still two iterates behind for the first-order methods, and Δf̃ above still reads it. That is recorded as open issue A10 in CHANGELOG.md.

invalidate_latest_gradient!(cache)

Declare that latest_gradient is no longer the gradient at solution(cache).

The default is a no-op, for the caches that never claim the pairing in the first place.

isaOptimizerState(alg)

Verify if an object implements the OptimizerState interface.

isconverged(status)

Whether any of the three convergence flags convergence_measures sets is set.

The flags are a disjunction on purpose: x_converged, f_converged and g_converged test different things and a solve is entitled to stop on any one of them. Note that solve! can also stop for reasons that are not convergence — the iteration cap, a non-finite iterate, an increase in f where one is not allowed — and none of those sets a flag here, so this is what tells the two apart; see meets_stopping_criteria.

x_converged is the one to be careful with: "the iterate stopped moving" is a statement about a ratio, and a diverging solve is where the denominator stops meaning anything. See the warning under convergence_measures for the two guards on it and for the one case they do not cover.

latest_gradient(cache)

The array holding the most recently evaluated gradient: a line search trial point while linesearch_problem's $\varphi'$ is being evaluated, and the accepted iterate once solver_step! has called refresh_latest_gradient! on it.

This is the array OptimizerStatus reports as rg, and it is deliberately not gradient: the latter is $\nabla{}f(x_k)$ at the iterate the step is built from, which the direction and the state updates need to keep reading unchanged.

Implementation

Every cache in this package carries a scratch array of its own for this and implements the three methods that go with it (refresh_latest_gradient!, latest_gradient_is_current, invalidate_latest_gradient!); see GradientCache for the field and for why it may not be shared with gradient.

The default here aliases the two together, which is what all six caches did before there was a name for the distinction, and is what a cache that does not refresh anything should keep doing: nothing sets the pairing store_gradient! relies on, latest_gradient_is_current defaults to false, and rg then means $\|\nabla{}f(x_k)\|$ at the iterate the step started from.

latest_gradient_is_current(cache, state, x)

Whether latest_gradient already holds $\mathrm{global\_rep}(\mathrm{section}(state), \nabla{}f(x))$, i.e. exactly what store_gradient! would otherwise evaluate.

The default is false: a cache that never refreshes latest_gradient has nothing to reuse.

lift_factors(B::AbstractLieAlgHorMatrix)

Factor a horizontal lift into the two $N\times{}2n$ matrices $B'$ and $B''$ with $\bar{B} = B'(B'')^T$:

\[\bar{B} = \begin{bmatrix} A & -B^T \\ B & \mathbb{O} \end{bmatrix} = \begin{bmatrix} \frac{1}{2}A & \mathbb{I} \\ B & \mathbb{O} \end{bmatrix} \begin{bmatrix} \mathbb{I} & \mathbb{O} \\ \frac{1}{2}A & -B^T \end{bmatrix} =: B'(B'')^T.\]

Every retraction in this package is built on this factorisation — it is what makes their cost scale with $n$ rather than with $N$, since the only matrix function either of them evaluates is on the $2n\times{}2n$ product $(B'')^TB'$. Both geodesic and cayley use it, for both manifolds.

For a GrassmannLieAlgHorMatrix this is the same expression with $A \equiv \mathbb{O}$.

lift_from_columns(B, V)

Rebuild a lift of the same type as B from the $N\times{}n$ block V, i.e. from the first $n$ columns of a skew matrix. Used by retraction_differential, which is the only place a lift has to be assembled from its action on $E$ rather than from its own blocks.

Not to be confused with the canonical horizontal lift GeometricOptimizers.Ω, which maps a tangent vector to an element of $\mathfrak{g}^\mathrm{hor}$. This one takes the columns of a matrix that is already in the Lie algebra.

linesearch_parameters(cache, x, state, c)

The parameters solver_step! hands the line search: the iterate x, the state whose section the trial steps retract from, and — where a manifold supplies a scale — the step ceiling αmax.

Implementation

Three methods, chosen on the type of the solution, as trial_iterate!'s two are.

For a Manifold, αmax is step_αmax(c, direction(cache)). This is the caller's half of the fix for issue A1b, and it has to be rebuilt at every solver step because $\|\delta\|$ changes at every solver step — including between the two searches of one solver_step!, since the second runs on a direction steepest_descent! has just replaced.

For a NamedTuple it is _manifold_αmax over the blocks, i.e. the same quantity derived per block and minimised over the manifold ones. A NamedTuple carrying no manifold block therefore gets Inf, which is what upstream reads as "the caller has no scale of its own" (SimpleSolvers.linesearch_αmax) and is exactly equivalent to the omission below.

For an AbstractVector the field is omitted, which is not the same as passing Inf only in spirit — upstream reads it through a hasproperty guard that constant-folds, so a caller that supplies nothing pays nothing. Omitting it is also the right answer and not merely the cheap one: a Euclidean $f(x + \alpha{}p)$ grows with $\alpha$, so the search's own sufficient-decrease test throws an over-long step out unaided, which is why this defect went unreported upstream for so long. There is no geometric scale to supply and the method's own SimpleSolvers.DEFAULT_LINESEARCH_αmax is the whole of the bound.

Every branch returns a concrete NamedTuple, so the merit closures stay type-stable. That is why the manifold-free NamedTuple passes Inf rather than joining the AbstractVector branch: deciding the shape of the parameters on the block types makes the return type a Union of two NamedTuples, which the merit closures then pay for on every evaluation.

linesearch_rejected(status)
linesearch_rejected(status, αmax)

Whether a line search reported that it could not decrease the merit along the direction it was given.

The second form takes the ceiling the caller passed as params.αmax and is what solver_step! uses; see A step the caller forbade is not a failed direction below.

SimpleSolvers.solve returns only the step length, so an outcome of LINESEARCH_FLOOR, LINESEARCH_EXHAUSTED or LINESEARCH_NO_DESCENT used to be indistinguishable from a successful search: the step came back and solver_step! took it. The three mean, in order, that the decrease achieved was no larger than the merit's own round-off resolution, that the budget ran out or the merit could not be bracketed, and that $\varphi'(0) \geq 0$. In none of them does the returned $\alpha$ carry a guarantee, and taking it anyway is how a solve can walk uphill.

The outcome comes from SimpleSolvers.solve_with_status, which solver_step! calls in place of solve for exactly this reason.

(Plain code and not an @extref: SimpleSolvers documents solve_with_status per method and not as a binding, so there is no binding-level entry in its inventory to link to. Documenter reports an unresolvable external link as an error, unlike the @ref-to-a-dead-signature case of issue D5, which it resolves silently to the wrong page.)

manifold_constructor(x::Manifold)

The one-argument constructor of x's manifold: GrassmannManifold for a GrassmannManifold{Float32, Matrix{Float32}}.

The type name and not typeof(x), because the array the result is applied to may have a different element type from x's: ParameterHandling.flatten's unflatten and the closure GradientAutodiff differentiates are both handed a vector of ForwardDiff.Duals. And x's manifold rather than a hardcoded StiefelManifold, which is what used to turn a GrassmannManifold into a StiefelManifold on a flatten round trip and to make a bare one a MethodError at Optimizer construction (issue A11).

Its three callers are ParameterHandling.flatten(::Type, ::Manifold), GeometricOptimizers._similar(::Manifold) and GradientAutodiff(F, ::Manifold).

manifold_type(B::AbstractLieAlgHorMatrix)

The manifold a retraction of B lands on.

$\mathfrak{g}^\mathrm{hor}$ is the horizontal component of the Lie algebra of a specific homogeneous space, so the lift already determines where the retraction maps to. This is what lets geodesic and cayley be written once for both manifolds rather than twice each.

meets_stopping_criteria(status, config, iterations)

Check if the optimizer has converged.

Implementation

meets_stopping_criteria checks if one of the following is true:

• converged (the output of SimpleSolvers.assess_convergence) is true and iterations $\geq$ config.min_iterations,
• if config.allow_f_increases is false: status.f_increased is true,
• iterations $\geq$ config.max_iterations,
• status.rfₐ $>$ config.f_abstol_break,
• any of status.x_nonfinite, status.f_nonfinite, status.g_nonfinite.

Extended help

metric(Y::GrassmannManifold, Δ₁::AbstractMatrix, Δ₂::AbstractMatrix)

Compute the metric for vectors Δ₁ and Δ₂ at Y.

The representation of the Grassmann manifold is realized as a quotient space of the Stiefel manifold.

The metric for the Grassmann manifold is:

\[g^{Gr}_Y(\Delta_1, \Delta_2) = g^{St}_Y(\Delta_1, \Delta_2) = \mathrm{Tr}(\Delta_1^T (\mathbb{I} - Y Y^T) \Delta_2) = \mathrm{Tr}(\Delta_1^T \Delta_2),\]

where we used that $Y^T\Delta_i$ for $i = 1, 2.$

metric(Y::StiefelManifold, Δ₁::AbstractMatrix, Δ₂::AbstractMatrix)

Compute the dot product for Δ₁ and Δ₂ at Y.

This uses the canonical Riemannian metric for the Stiefel manifold:

\[g_Y: (\Delta_1, \Delta_2) \mapsto \mathrm{Tr}(\Delta_1^T(\mathbb{I} - \frac{1}{2}YY^T)\Delta_2).\]

opnorm₁(X)

The induced 1-norm of X, i.e. its largest absolute column sum, as a reduction.

LinearAlgebra.opnorm(X, 1) is the natural spelling and is not used, because LinearAlgebra.opnorm1 is a double loop over X[i, j]. Scalar indexing is precisely what an array on a GPU backend cannot serve, and being free of it is the reason ScaledSquaring is the default algorithm — so the one norm that algorithm takes has to be expressible as sum and maximum, which every KernelAbstractions backend specializes.

The two agree to a few eps, not bitwise: opnorm1 accumulates each column sequentially in at least Float64, whereas sum is pairwise and accumulates in eltype(X). The value is only ever used to pick the number of halvings s = ⌈log₂(‖X‖₁/θ)⌉, so a difference of an ulp can at most shift s by one, and only for an argument that lands exactly on a power of two.

refresh_latest_gradient!(cache, gradient_instance)

Evaluate the gradient at the iterate cache currently holds and store it in latest_gradient.

solver_step! calls this once the accepted step has been taken, so that rg is a statement about the point the solve is about to report rather than about the one it started the step from.

All six caches in this package implement it; the default is a no-op so that a cache defined elsewhere keeps working, at the cost of an rg that means $\|\nabla{}f(x_k)\|$. See latest_gradient.

It is also what establishes the pairing store_gradient! relies on, so a cache that implements this has to implement latest_gradient_is_current with it.

restart!(state::BFGSState)

Reset the inverse-Hessian approximation to the identity, so that the next direction is the steepest-descent one.

This is the BFGSState method of restart!; DFPState is an alias for BFGSState, so it covers DFP as well. solver_step! calls it when a line search reports that it could not decrease the merit — see linesearch_rejected.

restart!(state)

Discard whatever curvature information state has accumulated and start again from the identity.

Only the quasi-Newton states carry any, so this is a no-op for everything else. See restart!(::BFGSState) and solver_step!.

retraction(R::AbstractRetraction, x)

Apply the retraction R to x, i.e. dispatch on the retraction type rather than calling geodesic or cayley by name.

This is what a caller who has been handed a retraction = … keyword uses: the two shipped types, Geodesic and Cayley, select the two functions, and a Geodesic also carries the AbstractExponentialAlgorithm its geodesic is evaluated with. R(x) is the same thing written as a call.

x is an AbstractLieAlgHorMatrix on a manifold — the retractions this package is about map $\mathfrak{g}^\mathrm{hor}\to{}G$ — a NamedTuple of parameters, or an ordinary array, on which every retraction is the identity because the extended retraction on a vector space is addition and update_section! does the adding.

An R and an x that do not go together is an error that says so, rather than a nothing that fails further downstream.

retraction_differential(retraction, B, α)

The horizontal lift $D(\alpha)$ that generates the velocity of the curve a line search walks along.

trial_iterate! builds its trial point as $x(\alpha) = \Lambda_0\mathrm{retract}(\alpha\bar{B})E$ with $\Lambda_0$ the section the step starts from, so with $W(\alpha) = \Lambda_0\mathrm{retract}(\alpha\bar{B})$ the frame the cache holds afterwards, this returns the $D(\alpha)$ for which

\[\frac{dx}{d\alpha} = W(\alpha)D(\alpha)E .\]

trial_slope pairs the gradient against it, which is what makes $\varphi'(\alpha)$ the derivative of $\varphi$.

Implementation

For Geodesic this is $\bar{B}$ itself. The exponential is a one-parameter subgroup, so $\frac{d}{d\alpha}\exp(\alpha\bar{B}) = \exp(\alpha\bar{B})\bar{B}$ and the frame factors out whole. The same holds for an ordinary array under either retraction, where the retraction is addition and $\varphi$ is affine in $\alpha$.

Cayley is not a one-parameter subgroup, and this is what that costs. With $M = (\mathbb{I} - \frac{\alpha}{2}\bar{B})^{-1}$,

\[\frac{d}{d\alpha}\mathrm{Cayley}(\alpha\bar{B}) = M\bar{B}M , \qquad D(\alpha) = \mathrm{Cayley}(\alpha\bar{B})^{-1}M\bar{B}M = M^T\bar{B}M ,\]

using $M^T = (\mathbb{I} + \frac{\alpha}{2}\bar{B})^{-1}$ for skew $\bar{B}$. $D$ is skew but not horizontal; only its first $n$ columns ever enter — $D(\alpha)E$ — and the horizontal projection has the same ones, so what lift_from_columns is handed is that projection.

For a GrassmannLieAlgHorMatrix the top $n\times{}n$ block is dropped rather than kept, and that is not an approximation. $M$, $M^T$ and $\bar{B}$ are all rational functions of $\bar{B}$ and so commute, which gives

\[D(\alpha) = \bar{B}\left(\mathbb{I} - \tfrac{\alpha^2}{4}\bar{B}^2\right)^{-1} ,\]

and for a Grassmann lift — where $A \equiv \mathbb{O}$, so $\bar{B}^2$ is block-diagonal — the top-left block of that product is identically zero at every $\alpha$. The drop is exact whatever $D$ is paired against, not merely exact against a gradient that happens to be orthogonal to $Y$.

Both inverses are taken through lift_factors and the Woodbury identity, exactly as cayley does, so the cost is $O(Nn^2 + n^3)$ and no $N\times{}N$ matrix is formed.

α = 0 returns B unchanged, and that is the case that matters for cost: SimpleSolvers.Backtracking — the default of default_linesearch — evaluates $\varphi'$ only at $\alpha = 0$, so it never reaches the general branch.

rgrad(Y::GrassmannManifold, ∇L::AbstractMatrix)

Compute the Riemannian gradient for the Grassmann manifold at Y based on ∇L.

Here $Y$ is a representation of $\mathrm{span}(Y)\in{}Gr(n, N)$ and $\nabla{}L\in\mathbb{R}^{N\times{}n}$ is the Euclidean gradient.

This gradient has the property that it is orthogonal to the space spanned by $Y$.

The precise form of the mapping is:

\[\mathtt{rgrad}(Y, \nabla{}L) \mapsto \nabla{}L - YY^T\nabla{}L.\]

Note the property $Y^T\mathrm{rgrad}(Y, \nabla{}L) = \mathbb{O}.$

Also see rgrad(::StiefelManifold, ::AbstractMatrix).

Examples

using GeometricOptimizers

Y = GrassmannManifold([1 0 ; 0 1 ; 0 0; 0 0])
Δ = [1 2; 3 4; 5 6; 7 8]
rgrad(Y, Δ)

# output

4×2 Matrix{Int64}:
 0  0
 0  0
 5  6
 7  8

rgrad(Y::StiefelManifold, ∇L::AbstractMatrix)

Compute the Riemannian gradient for the Stiefel manifold at Y based on ∇L.

Here $Y\in{}St(N,n)$ and $\nabla{}L\in\mathbb{R}^{N\times{}n}$ is the Euclidean gradient.

The function computes the Riemannian gradient with respect to the canonical metric: metric(::StiefelManifold, ::AbstractMatrix, ::AbstractMatrix).

The precise form of the mapping is:

\[\mathtt{rgrad}(Y, \nabla{}L) \mapsto \nabla{}L - Y(\nabla{}L)^TY\]

Note the property $Y^T\mathtt{rgrad}(Y, \nabla{}L)\in\mathcal{S}_\mathrm{skew}(n).$

Examples

using GeometricOptimizers

Y = StiefelManifold([1 0 ; 0 1 ; 0 0; 0 0])
Δ = [1 2; 3 4; 5 6; 7 8]
rgrad(Y, Δ)

# output

4×2 Matrix{Int64}:
 0  -1
 1   0
 5   6
 7   8

rhs(cache)

Return the right hand side of an instance of BFGSCache

rhs(cache)

Return the right hand side of an instance of DFPCache

rhs(cache)

Return the right hand side of an instance of NewtonOptimizerCache

solution_scale(x)

The norm the relative change in the iterate is measured against, i.e. the denominator of $\|x - x'\|/\|x'\|$.

For an ordinary array this is the l2norm of the iterate, which is the only scale a Euclidean problem supplies. For a Manifold it is not: a point of $St(N, n)$ or $Gr(n, N)$ satisfies $Y^TY = \mathbb{I}_n$, so $\|Y\|_F = \sqrt{\mathrm{tr}(Y^TY)} = \sqrt{n}$ exactly, and the scale is a constant the geometry supplies rather than something to be measured off the iterate.

That distinction is what keeps x_converged honest on a solve that has diverged; see convergence_measures. While the iterate is on the manifold the two agree to round-off, so this changes no measurement of any converging solve — it changes what is reported once the iterate is somewhere $\sqrt{n}$ no longer describes.

A NamedTuple combines its blocks in quadrature, exactly as l2norm does, so a NamedTuple holding a manifold next to an ordinary array uses the nominal scale for the one and the measured scale for the other.

solver_step!(x, state, opt)

Compute a full iterate for an Optimizer.

Examples

julia> f(x) = sum(x .^ 2 + x .^ 3 / 3);

julia> x = [1f0, 2f0]
2-element Vector{Float32}:
 1.0
 2.0

julia> opt = Optimizer(x, f; algorithm = Newton());

julia> state = NewtonOptimizerState(x);

julia> update!(state, gradient(opt), x);

julia> solver_step!(x, state, opt)
2-element Vector{Float32}:
 0.25
 0.6666666

Extended help

steepest_descent!(cache)

Replace the direction stored in cache by $-\nabla{}f$, the one direction that always descends.

Implementation

rhs is $-\nabla{}f$ on the (quasi-)Newton caches, so the default copies it across — which is what solver_step! did inline before this existed.

It is not on the three first-order caches, where rhs is an alias for direction itself, so copying would be a silent no-op. GradientCache's direction already is $-\nabla{}f$ and the no-op is correct there; MomentumCache and AdamCache hold a moving average and have to be given the gradient explicitly.

step_size(method, t)

The step size a DecayingStatic takes in iteration t, where t is the iteration number the solve reports — so t = 1 for the first step, and t = 0 evaluates the schedule at a point no solve asks for (see the remark on $\alpha(0)$ in DecayingStatic).

step_αmax(c, δ)

The ceiling on the line search's $\alpha$ that a step ceiling of c imposes on the direction δ: $c\,2\pi/\|\delta\|$.

c is in multiples of $2\pi$ because that is the scale the geometry supplies. Retracting a horizontal lift is a rotation, so $\alpha\mapsto\mathrm{retract}(\alpha\bar{B})$ has nothing left to reach once $\|\alpha\bar{B}\|$ passes $2\pi$ — Cayley has converged to a fixed rotation by then and Geodesic is periodic. Everything past it is round-off, and on the SVD problem enough of it to leave the manifold; see DEFAULT_STEP_CEILING and issue A1b.

δ is one direction and not a NamedTuple of them: the $2\pi$ is a property of a rotation, so a solution built of several blocks needs one ceiling per block and the smallest of them, which is _manifold_αmax's job. See linesearch_parameters.

Implementation

A direction whose norm is zero or not finite yields Inf, which SimpleSolvers.linesearch_αmax reads as "the caller has no scale of its own" and which leaves the method's own ceiling standing. That is a guard and not a formality: upstream raises an ArgumentError on a NaN or non-positive params.αmax before it evaluates the merit — correctly, since silently ignoring one would hand back exactly the unbounded step the caller was ruling out — so a vanishing direction has to be special-cased here rather than there.

store_gradient!(cache, state, gradient_instance, x)

Put $\mathrm{global\_rep}(\mathrm{section}(state), \nabla{}f(x))$ into gradient_array(cache), which is what the three first-order update!(cache, ...) methods build their direction from.

Implementation

This reuses latest_gradient when latest_gradient_is_current says it already holds that value, which in a solve! loop is every iteration but the first: solver_step! refreshes it at the accepted iterate, and the next update! is asked for the gradient at that same iterate in the same frame. The two are not merely close, they are the same computation – update_section!(Λᵗ, Λ⁽ᵗ⁻¹⁾, B, retraction) has the body the two-argument form update!(::MomentumState, ...) uses, so section(cache) after solver_step! is bit-for-bit section(state) after update!(state, opt, x).

Without the reuse the refresh doubles the gradient evaluations of a first-order step: on the SVD problem of test/optimizer_convergence/svd_optim.jl, Adam + Static over 2 000 iterations costs 124 ms without the refresh, 167 ms with it and 128 ms with it reused — one trajectory throughout. With the reuse, the refresh is the step's gradient evaluation.

trace(result)

The per-iteration record solve! kept, or an empty vector if it was not asked to keep one.

Options(store_trace = true) is what asks. Before this existed the option was accepted and silently ignored — by this package and by SimpleSolvers 0.11, where it is a field of SimpleSolvers.Options that nothing reads — so code that set it got no trace and no error either.

Examples

There is one entry per iteration, and the last one agrees with the status the solve reports:

julia> f(x) = sum(x .^ 2);

julia> x = [1.0, 2.0];

julia> state = OptimizerState(Newton(), x);

julia> result = solve!(x, state, Optimizer(x, f; algorithm = Newton(), store_trace = true));

julia> length(trace(result)) == iteration_number(state)
true

julia> last(trace(result)).rg == status(result).rg
true

julia> first(trace(result)).iteration
1

Without the option there is no trace, and no error either:

julia> x = [1.0, 2.0];

julia> result = solve!(x, OptimizerState(Newton(), x), Optimizer(x, f; algorithm = Newton()));

julia> isempty(trace(result))
true

Implementation

The entries are taken from the OptimizerStatus that solve! already computes on every iteration, so nothing extra is evaluated and the cost when store_trace is unset is one Bool test per iteration.

What it is for: a statistic that has to be phase-independent. Adam at a fixed learning rate does not converge to the minimizer, it orbits it at a distance of order $\alpha$, so the error at any one iteration is a sample of an arbitrary phase on that orbit and moves with the last bits of the floating-point arithmetic. Averaging over a stretch of the orbit measures its radius instead, which is a property of $\alpha$ and the problem — across Julia 1.10, 1.12 and 1.13 the final-iterate error on the problem in test/optimizer_convergence/svd_optim.jl spans a factor of 3.0 and the mean over the last five hundred iterations spans 1.06.

trial_iterate!(cache, params, α, retraction)

Write the iterate a step of length α along the current direction would produce into solution(cache), and return it. This is what the merit function of the line search evaluates.

Implementation

There are two methods, chosen on the type of the solution.

For an AbstractVector the iterate is $x_k + \alpha{}p_k$, i.e. SimpleSolvers.compute_new_iterate!, which is what this used to do unconditionally.

For a Manifold — or a NamedTuple that contains one — that is both undefined and wrong: adding $\alpha{}p_k$ to a point of $St(N, n)$ leaves the manifold, and the direction is an AbstractLieAlgHorMatrix of a different shape than the point to begin with. The step has to go through the retraction, exactly as solver_step! does once the line search has picked its α, so this reproduces that: scale the direction, push the state's section through the retraction into the cache's section, and read the point back out.

trial_slope(gradient_instance, cache, retraction, α)

Return $\varphi'(\alpha) = \langle\nabla{}f(x(\alpha)), D(\alpha)\rangle$ at the iterate currently held in solution(cache), for the derivative of the line search's merit. $D(\alpha)$ is retraction_differential(retraction, direction(cache), α).

Implementation

On a Manifold the gradient that gradient_instance produces lives in the tangent space of the embedding while direction(cache) is a horizontal lift, so the two are not even the same shape and have to be brought together by global_rep first — which is the same map update!(::GradientCache, …) applies to the gradient before storing it.

The pairing is then _dot and not dot: the α of the line search parameterizes a curve in the intrinsic coordinates of the lift, and dot on a lift is the ambient Frobenius product, which is exactly twice that. With dot this returned 2\varphi'(\alpha) — invisible to SimpleSolvers.Bisection, which only looks for a sign change, and swamped by the Armijo slack in SimpleSolvers.Backtracking, but wrong, and wrong in a way a curvature condition or a quadratic fit would act on.

The direction the gradient is paired with is retraction_differential(retraction, B, α) and not B itself. $\varphi'(\alpha) = \langle\nabla{}f(x(\alpha)), B\rangle$ holds only where $\alpha \mapsto \mathrm{retract}(\alpha{}B)$ is a one-parameter subgroup, which Geodesic is and Cayley is not; under Cayley the generator of the curve's velocity turns with the step, and the differential is what supplies it. Against a central difference of the merit this is now exact for both retractions at every $\alpha$.

unit_matrix(backend, T, n)
unit_matrix(A::AbstractMatrix)

The $n\times{}n$ identity of element type T on backend, with the diagonal written by write_ones_kernel! above.

Base.one(::AbstractMatrix) is the natural spelling and is not used: Base._one allocates and then writes the diagonal in a scalar-indexed loop, which is precisely what an array on a KernelAbstractions backend cannot serve. It is the same hazard GeometricOptimizers.opnorm₁ exists to avoid one level up, and every identity this package builds goes through here rather than through Base.one — Base.one for SymmetricMatrix, SkewSymMatrix, the two AbstractTriangulars and AbstractLieAlgHorMatrix, and the $2n\times{}2n$ identities that GeometricOptimizers.𝔄 and NativePade need.

LinearAlgebra.I covers some of those uses and is not enough either. GPUArrays supplies a kernel-based +(::AbstractGPUMatrix, ::UniformScaling), so X + I is portable on the array types that package covers — but KernelAbstractions is the interface this package is written against, and a backend of its own is under no obligation to be one of them.

The matrix form takes the backend and the element type from A and its size from LinearAlgebra.checksquare, so it throws on a non-square argument exactly as Base.one does.

update_section!(Λᵗ, Λ⁽ᵗ⁻¹⁾, B⁽ᵗ⁻¹⁾, retraction)
update_section!(Λ, B, retraction)

Transport the GlobalSection Λ⁽ᵗ⁻¹⁾ along the step B⁽ᵗ⁻¹⁾ and write the result into Λᵗ:

\[\Lambda^{(t)} \leftarrow \Lambda^{(t-1)}\mathrm{Retraction}(B^{(t-1)}).\]

This is the fourth of the five steps an optimizer step consists of on a homogeneous space — see Optimization on Homogeneous Spaces — and the one that makes the cache independent of the iterate: the new point is read back out of the section afterwards with apply_section!, rather than the section being recomputed at the new point.

B⁽ᵗ⁻¹⁾ is the final velocity the optimizer method produced, already scaled by the step length. It is an AbstractLieAlgHorMatrix when the parameter is on a manifold, and of the parameter's own type when it is not — a vector-space parameter carries no $\lambda$, and the extended retraction on a vector space is addition, so retraction is then ignored.

The three-argument form is the two-argument section written in place, Λᵗ === Λ⁽ᵗ⁻¹⁾. A NamedTuple of parameters is walked leaf by leaf.

value(obj::AbstractOptimizerProblem, x)

Evaluates the value at x (i.e. computes obj.F(x)).

Ω(Y::GrassmannManifold{T}, Δ::AbstractMatrix{T}) where T

Perform the canonical horizontal lift for the Grassmann manifold:

\[ \Delta \mapsto \Omega^{St}(\Delta),\]

where $\Omega^{St}$ is the canonical horizontal lift for the Stiefel manifold.

using GeometricOptimizers
using GeometricOptimizers: StiefelProjection

E = GrassmannManifold(StiefelProjection(5, 2))
Δ = [0. 0.; 0. 0.; 2. 3.; 4. 5.; 6. 7.]
GeometricOptimizers.Ω(E, Δ)

# output

5×5 SkewSymMatrix{Float64, Vector{Float64}}:
 0.0  -0.0  -2.0  -4.0  -6.0
 0.0   0.0  -3.0  -5.0  -7.0
 2.0   3.0   0.0  -0.0  -0.0
 4.0   5.0   0.0   0.0  -0.0
 6.0   7.0   0.0   0.0   0.0

Ω(Y::StiefelManifold{T}, Δ::AbstractMatrix{T}) where T

Perform canonical horizontal lift for the Stiefel manifold:

\[ \Delta \mapsto (\mathbb{I} - \frac{1}{2}YY^T)\Delta{}Y^T - Y\Delta^T(\mathbb{I} - \frac{1}{2}YY^T).\]

Internally this performs

SkewSymMatrix(2 * (I(n) - .5 * Y * Y') * Δ * Y')

It uses SkewSymMatrix to save memory.

Examples

using GeometricOptimizers
using GeometricOptimizers: StiefelProjection

E = StiefelManifold(StiefelProjection(5, 2))
Δ = [0. -1.; 1. 0.; 2. 3.; 4. 5.; 6. 7.]
GeometricOptimizers.Ω(E, Δ)

# output

5×5 SkewSymMatrix{Float64, Vector{Float64}}:
 0.0  -1.0  -2.0  -4.0  -6.0
 1.0   0.0  -3.0  -5.0  -7.0
 2.0   3.0   0.0  -0.0  -0.0
 4.0   5.0   0.0   0.0  -0.0
 6.0   7.0   0.0   0.0   0.0

Note that the output of Ω is a skew-symmetric matrix, i.e. an element of $\mathfrak{g}$.

𝔄(B̂, B̄)
𝔄(B̂, B̄, algorithm)

Compute $\mathfrak{A}(B', B'') := \sum_{n=1}^\infty \frac{1}{n!} ((B'')^TB')^{n-1}.$

This expression has the property $\mathbb{I} + B'\mathfrak{A}(B', B'')(B'')^T = \exp(B'(B'')^T).$

Note that the argument $(B'')^TB'$ is only $2n\times{}2n$, so this is where the cost of a retraction is set by $n$ rather than by $N$.

Examples

using GeometricOptimizers
using GeometricOptimizers: 𝔄
import Random
Random.seed!(123)

B = rand(StiefelLieAlgHorMatrix, 10, 2)
B̂ = hcat(vcat(.5 * B.A, B.B), vcat(one(B.A), zero(B.B)))
B̄ = hcat(vcat(one(B.A), zero(B.B)), vcat(-.5 * B.A, -B.B))

one(B̂ * B̄') + B̂ * 𝔄(B̂, B̄) * B̄' ≈ exp(Matrix(B))

# output

true

𝔄(X, algorithm)

Compute $\mathfrak{A}(X)$ with the requested AbstractExponentialAlgorithm.

All algorithms compute the same function and differ only in accuracy at a large X, in cost, and in which backends they run on. See AbstractExponentialAlgorithm for the comparison and ScaledSquaring, which is the default.

Examples

The five agree wherever the unscaled series is still accurate, and only four of them agree beyond that:

using GeometricOptimizers
using GeometricOptimizers: 𝔄, ScaledSquaring, NativePade, AugmentedPade, TaylorSeries
import Random
Random.seed!(123)

X = randn(6, 6)

isapprox(𝔄(X, ScaledSquaring()), 𝔄(X, NativePade()); rtol = 1e-12) &&
    isapprox(𝔄(X, NativePade()), 𝔄(X, AugmentedPade()); rtol = 1e-12) &&
    isapprox(𝔄(X, ScaledSquaring()), 𝔄(X, TaylorSeries()); rtol = 1e-12)

# output

true

𝔄(A)

Compute $\mathfrak{A}(A) := \sum_{n=1}^\infty \frac{1}{n!} (A)^{n-1}.$

Implementation

This uses a Taylor expansion that iteratively adds terms with

while norm(Aⁿ) > ε
mul!(A_temp, Aⁿ, A)
Aⁿ .= A_temp
rmul!(Aⁿ, T(inv(n)))

𝔄A += Aⁿ
n += 1
end

until the norm of Aⁿ becomes smaller than machine precision. The counter n in the above algorithm is initialized as 2 The matrices Aⁿ and 𝔄 are initialized as the identity matrix.

𝔄exp(B̂, B̄, algorithm = ScaledSquaring())

Compute $\exp(B'(B'')^T)$ as $\mathbb{I} + B'\mathfrak{A}(B', B'')(B'')^T$, i.e. the identity 𝔄 exists for, packaged as the exponential it computes.

This is what makes a geodesic retraction cheap: the argument handed to 𝔄 is $(B'')^TB'$, which is $2n\times{}2n$, so the cost is set by $n$ and not by $N$ even though the result is $N\times{}N$. geodesic computes this same product inline — it needs to wrap the result in manifold_type(B) and to take the lift factors apart itself — so this is for callers that want the matrix exponential of a low-rank product on its own.

algorithm is forwarded to 𝔄, which supplies TaylorSeries, ScaledSquaring, NativePade and AugmentedPade. ProjectedSkew is not among them: it is a geodesic-level algorithm with its own branch there and no 𝔄 method, so 𝔄exp(B̂, B̄, ProjectedSkew()) fails inside 𝔄 exactly as 𝔄(B̂, B̄, ProjectedSkew()) does.

Implementation

The default is ScaledSquaring and not the unscaled series, for the reason given under geodesic: the series cancels catastrophically once $\|\bar{B}\| \gtrsim 50$, which is not a regime a function that presents itself as an exponential may quietly get wrong. Relative error against exp(Matrix(B)) for B = scale * rand(StiefelLieAlgHorMatrix, 10, 2), as test/retractions/exponential_accuracy.jl draws it:

Note that this differs from 𝔄(B̂, B̄), which has no algorithm argument at all and is the unscaled series. That one is a kernel, and 𝔄 carries the warning; this one is a result.

Examples

using GeometricOptimizers
using GeometricOptimizers: 𝔄exp
import Random
Random.seed!(123)

B = rand(StiefelLieAlgHorMatrix, 10, 2)
B̂ = hcat(vcat(.5 * B.A, B.B), vcat(one(B.A), zero(B.B)))
B̄ = hcat(vcat(one(B.A), zero(B.B)), vcat(-.5 * B.A, -B.B))

𝔄exp(B̂, B̄) ≈ exp(Matrix(B))

# output

true

direction(cache)

Return the direction of the gradient step (i.e. Δx) of an instance of BFGSCache.

direction(cache)

Return the direction of the gradient step (i.e. Δx) of an instance of DFPCache.

direction(cache)

Return the direction of the gradient step (i.e. Δx) of an instance of NewtonOptimizerCache.

linesearch_problem(problem, gradient, cache, retraction)

Create a SimpleSolvers.LinesearchProblem for the linesearch algorithm.

The variable on which this problem depends is $\alpha$. The trial iterate is built by trial_iterate!, which is what makes this work for manifold parameters and not only for AbstractVectors.

Example

julia> x = [1, 0., 0.]
3-element Vector{Float64}:
 1.0
 0.0
 0.0

julia> f = x -> sum(x .^ 3 / 6 + x .^ 2 / 2);

julia> obj = OptimizerProblem(f, x);

julia> grad = GradientAutodiff{Float64}(obj.F, length(x));

julia> hess = HessianAutodiff{Float64}(obj.F, length(x));

julia> cache = NewtonOptimizerCache(x);

julia> state = NewtonOptimizerState(x); update!(state, grad, x);

julia> params = (x = state.x, state = state);

julia> update!(cache, state, grad, hess, x);

julia> ls_obj = linesearch_problem(obj, grad, cache, Cayley());

julia> ls_obj.F(0., params)
0.6666666666666666

julia> ls_obj.D(0., params)
-1.125

outer!(m, g₁, g₂)

The outer product of two horizontal lifts, written into m.

Like the ArrayNamedTuple method above, this exists because the quasi-Newton Q is sized by the intrinsic dimension of the parameters — the length of their flattening — while the direction and the gradient are handed around in the ambient horizontal-lift representation. For a bare StiefelManifold of size (3, 1) those are 2 and 3 × 3 respectively, so SimpleSolvers.outer!, which indexes its arguments linearly against axes(m), would assert on the mismatch. Flattening first is what the NamedTuple case has always done; without the same method here BFGS and DFP cannot run on a bare Manifold at all.

solve!(x, state, opt)

Solve the optimization problem described by opt::Optimizer and store the result in x.

Examples

julia> f(x) = sum(x .^ 2 + x .^ 3 / 3);

julia> x = [1f0, 2f0]
2-element Vector{Float32}:
 1.0
 2.0

julia> opt = Optimizer(x, f; algorithm = Newton());

julia> state = NewtonOptimizerState(x);

julia> solve!(x, state, opt)
GeometricOptimizers.OptimizerResult{Float32, Float32, Vector{Float32}, GeometricOptimizers.OptimizerStatus{Float32, Float32}}( * Convergence measures

    |x - x'|               = 7.82e-03
    |x - x'|/|x'|          = 2.56e+02
    |f(x) - f(x')|         = 6.18e-05
    |f(x) - f(x')|/|f(x')| = 6.63e+04
    |g(x) - g(x')|         = 1.57e-02
    |g(x)|                 = 6.10e-05
, Float32[4.6478817f-8, 3.0517578f-5], 9.313341f-10, GeometricOptimizers.OptimizerTraceEntry{Float32, Float32}[])

julia> x
2-element Vector{Float32}:
 4.6478817f-8
 3.0517578f-5

julia> iteration_number(state)
4

The trailing empty vector is the per-iteration trace, which is only filled if Options(store_trace = true) asked for it.

Two rows of that status moved in 0.2.0 and both are the same fix (issue A8). |g(x)| is now $\|\nabla{}f\|$ at the iterate the solve returns, so g_converged fires when the residual has actually reached f_reltol — this solve stops on iteration 4 where it used to run a fifth and report the residual of the fourth. And |g(x) - g(x')| used to be 0.00e+00 structurally for Newton: solver_step! advanced state.ḡ at the same iterate the cache took its gradient at, so the difference it printed could not be anything else. See gradient_difference!.

Also see solver_step!.