Optimization

Filtering and Smoothing

class pystopt.filtering.FilterType(value, names=None, *, module=None, qualname=None, type=None, start=1, boundary=None)

Bases: Enum

Full = 1

Normal = 2

Unfiltered = 3

pystopt.filtering.apply_filter_spectral(actx: ArrayContext, discr: SpectralDiscretization, x: ArrayOrContainerT, *, method: str = 'ideal', **kwargs: Any) → ArrayOrContainerT

Apply a spectral filter to an array x.

For curves in two-dimensional space, the following Fourier-based filters are available

• "ideal" uses apply_filter_fourier_ideal().

• "trapezoidal" uses apply_filter_fourier_trapezoidal().

• "butterworth" uses apply_filter_fourier_butterworth().

For surfaces in three-dimensional space, the following Spherical Harmonic-based filters are available

• "ideal" uses apply_filter_spharm_ideal().

• "tikhonov" uses apply_filter_spharm_tikhonov().

• "tikhonov_opt" uses apply_filter_spharm_tikhonov_opt().

• "tikhonov_ideal" uses apply_filter_spharm_tikhonov_ideal().

Parameters:

• x – array container to filter.

• method – filtering method to use.

pystopt.filtering.apply_normal_gradient_filter_spectral(actx: ArrayContext, places: pytential.GeometryCollection, g: ndarray, *, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None, only_normal: bool = False, method: str = 'ideal', **kwargs: Any) → ndarray

Wrapper around apply_filter_spectral() that can be used to only filter the normal component of the vector g.

This is mainly meant for shape optimization, where the gradient is of the norm \(g \mathbf{n}\). Filtering the whole vector g would result in a gradient that does not point purely in the normal direction and can unnaturally distort the mesh.

Parameters:

only_normal – if True, the tangential component is discarded completely after filtering.

pystopt.filtering.apply_filter_by_type(actx: ArrayContext, places: pytential.GeometryCollection, ary: ArrayOrContainerT, *, filter_type: FilterType = FilterType.Unfiltered, filter_arguments: Dict[str, Any] | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None, force_spectral: bool = False) → ArrayOrContainerT

pystopt.filtering.modal_decay_estimate(actx: ArrayContext, discr: SpectralDiscretization, x: ArrayOrContainerT) → float

pystopt.filtering.interp_error_estimate(actx: ArrayContext, discr: SpectralDiscretization, x: ArrayOrContainerT) → ArrayOrContainerT

Yukawa Kernel Smoothing

pystopt.filtering.apply_yukawa_smoothing(actx: ArrayContext, places: pytential.GeometryCollection, x: ArrayOrContainerT, *, method: str = 'beltrami', strength: float | None = None, gmres_arguments: Dict[str, Any] | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ArrayOrContainerT

Filter a given quantity on a surface.

The filtering is done by solving the Yukawa-Beltrami equation

\[-\Delta_\Sigma y + \frac{1}{k^2} y = x.\]

where the input strength is \(k^2\) on a surface \(\Sigma\). A smaller strength means less smoothing.

Parameters:

• x – array container with DOFArray leaves.

• strength – strength of the smoothing in \((0, \infty]\). If it is \(\infty\), then the equation reduces to a Laplace-Beltrami equation.

pystopt.filtering.apply_yukawa_filter(actx: ArrayContext, places: pytential.GeometryCollection, x: ArrayOrContainerT, *, strength: float | None = None, gmres_arguments: Dict[str, Any] | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ArrayOrContainerT

Filter a given quantity on a surface.

Parameters:

• x – array container with DOFArray leaves.

• strength – strength of the smoothing in \((0, \infty]\). If it is \(\infty\), then the equation reduces to a Laplace-Beltrami equation.

Fourier Spectral Filtering

pystopt.filtering.apply_filter_fourier_ideal(actx: ArrayContext, discr: SpectralDiscretization, xk: arraycontext.with_container_arithmetic, *, kmax: int | Sequence[int] | None = None) → arraycontext.with_container_arithmetic

Ideal filter in Fourier space based on

\[\begin{split}\hat{f}_k = \begin{cases} f_k, & \quad |k| < k_{max}, \\ 0, & \quad \text{otherwise}. \end{cases}\end{split}\]

Parameters:

• x – array to filter.

• kmax – cutoff wavenumber above which all coefficients are set to zero.

pystopt.filtering.apply_filter_fourier_trapezoidal(actx: ArrayContext, discr: SpectralDiscretization, xk: arraycontext.with_container_arithmetic, *, ka: int | Sequence[int] | None = None, kb: int | Sequence[int] | None = None) → arraycontext.with_container_arithmetic

Trapezoidal filter in Fourier space based on

\[\begin{split}\hat{f}_k = \begin{cases} f_k, & \quad |k| <= k_a, \\ \frac{k_b - k}{k_b - k_a} f_k, & \quad k_a < |k| < k_b, \\ 0, & \quad |k| >= k_b. \end{cases}\end{split}\]

Parameters:

• x – array to filter.

• ka – wavenumber at which the ramp should start.

• kb – wavenumber at which the ramp should end, above which all coefficients are set to zero.

pystopt.filtering.apply_filter_fourier_butterworth(actx: ArrayContext, discr: SpectralDiscretization, xk: arraycontext.with_container_arithmetic, *, kmax: int | Sequence[int] | None = None, p: float | Sequence[float] | None = None) → arraycontext.with_container_arithmetic

Butterworth filter in Fourier space based on

\[\hat{f}_k = \left[1 + \left(\frac{k}{k_{max}}\right)^p\right] f_k.\]

where a smaller \(k_{max}\) has the effect of localizing the wavenumbers around \(0\) and a larger \(p\) has the effect of making the transition sharper.

Parameters:

• x – array to filter.

• kmax – cutoff wavenumber above which all coefficients are set to zero.

Spherical Harmonic Spectral Filtering

pystopt.filtering.apply_filter_spharm_ideal(actx: ArrayContext, discr: SpectralDiscretization, xk: arraycontext.with_container_arithmetic, *, kmax: int | Sequence[int] | None = None, p: float | str | None = None) → arraycontext.with_container_arithmetic

Ideal low-pass filtering for Spherical Harmonic coefficients.

All coefficients for which \((n^p + m^p)^{1/p} > k_{max}\) are set to zero, where \(n\) is the order and \(m\) is the degree of the corresponding spherical harmonic.

Parameters:

• kmax – cut-off for spherical harmonic orders.

• p – order of the norm used for the spherical harmonic indices.

pystopt.filtering.apply_filter_spharm_tikhonov(actx: ArrayContext, discr: SpectralDiscretization, xk: arraycontext.with_container_arithmetic, *, alpha: float | None = None, p: int | None = None) → arraycontext.with_container_arithmetic

Tikhonov regularization based filtering, obtained by optimizing

\[\frac{1}{2} \|f - I[f]\|_2 + \frac{\alpha}{2} \|\Delta^p_S I[f]\|_2,\]

where \(\Delta_S\) is the Laplacian on the unit sphere and \(I[f]\) is the projection of \(f\) into the spherical harmonic basis. The solution filters as

\[\hat{f}^m_n = \frac{f^m_n}{1 + \alpha n^{2p} (n + 1)^{2 p}}.\]

Parameters:

• alpha – weight used in the regularization; a larger alpha implies a smaller cutoff region in the filter.

• p – order of the Laplacian used in the regularization; a larger p implies a steeper cutoff in the filter.

pystopt.filtering.apply_filter_spharm_tikhonov_opt(actx: ArrayContext, discr: SpectralDiscretization, xk: arraycontext.with_container_arithmetic, *, p: int | None = None, alpha: float | None = None, maxiter: int = 32, optimize: str | None = None, method: str | None = None) → arraycontext.with_container_arithmetic

Tikhonov regularization based filtering, where the weight is obtained by Generalized Cross Validation.

Parameters:

• maxiter – maximum number of iterations to perform in the optimization of the GCV functional.

• optimize – optimization method to use, can be one of "minimize_scalar" or "minimize" from scipy.

• method – optimization method to use when optimize is "minimize".

pystopt.filtering.apply_filter_spharm_tikhonov_ideal(actx: ArrayContext, discr: SpectralDiscretization, xk: arraycontext.with_container_arithmetic, *, kmax: int | Sequence[int] | None = None, alpha: float | None = None, p: int | None = None) → arraycontext.with_container_arithmetic

Applies a combination of an ideal filter and a Tikhonov regularization.

The ideal filter is only applied in \(m\) and the Tikhonov regularization, by definition, is applied in \(n\). This allows handling the two sources of error coming from the DG interpolation.

See apply_filter_spharm_ideal() and apply_filter_spharm_tikhonov().

Manifolds and Shape Spaces

This module implements some helpers to handle Riemannian optimization in shape spaces. For an introduction to these concepts see [RingWirth2012].

[RingWirth2012] (1,2,3,4)

W. Ring, B. Wirth, Optimization Methods on Riemannian Manifolds and Their Application to Shape Space, SIAM Journal on Optimization, Vol. 22, pp. 596–627, 2012, DOI.

The abstract classes used to represent the manifold are

class pystopt.optimize.RiemannianManifold

Bases: object

Generic Riemannian manifold.

ambient_dim: int

Dimension of the ambient embedding space.

dim: int

Dimension of the manifold.

class pystopt.optimize.ShapeSpace(places: pytential.GeometryCollection, dofdesc: pytential.symbolic.primitives.DOFDescriptor)

Bases: RiemannianManifold

places: pytential.GeometryCollection

A geometry collection containing the surfaces.

dofdesc: pytential.symbolic.primitives.DOFDescriptor

A descriptor for the surface considered in the shape space.

They implement the following methods through a functools.singledispatch() mechanism.

pystopt.optimize.norm(m: RiemannianManifold, x: ArrayOrContainerT, u: ArrayOrContainerT) → Array

pystopt.optimize.norm(m: ShapeSpace, x: ArrayOrContainerT, u: ArrayOrContainerT) → Array

pystopt.optimize.norm(m: CartesianEucledianSpace, x: ArrayOrContainerT, u: ArrayOrContainerT) → Array

Parameters:

• x – a point on the manifold m;

• u – a vector in the tangent space of m at x.

Returns:

the norm of the vector u using the metric() on m.

pystopt.optimize.metric(m: RiemannianManifold, x: ArrayOrContainerT, u: ArrayOrContainerT, v: ArrayOrContainerT) → Array

pystopt.optimize.metric(m: CartesianEucledianSpace, x: ArrayOrContainerT, u: ArrayOrContainerT, v: ArrayOrContainerT) → Array

pystopt.optimize.metric(m: L2ShapeSpace, x: ArrayOrContainerT, u: ArrayOrContainerT, v: ArrayOrContainerT) → Array

pystopt.optimize.metric(m: WeightedL2ShapeSpace, x: ArrayOrContainerT, u: ArrayOrContainerT, v: ArrayOrContainerT) → ArrayOrContainerT

Parameters:

• x – a point on the manifold m;

• u – a vector in the tangent space of m at x.

• v – a vector in the tangent space of m at x.

Returns:

the value of \(g_x(u, v)\), where \(g_x\) is the metric at a point \(x\). This essentially reduces to an inner product in most cases.

pystopt.optimize.representation(m: RiemannianManifold, x: ArrayOrContainerT, u: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.representation(m: CartesianEucledianSpace, x: ArrayOrContainerT, u: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.representation(m: WeightedL2ShapeSpace, x: ArrayOrContainerT, u: ArrayOrContainerT) → ArrayOrContainerT

Parameters:

• x – a point on the manifold m;

• u – a vector in the tangent space of m at x.

Returns:

a representation of the vector u in the manifold. This essentially corresponds to “raising the index” in tensor calculus.

pystopt.optimize.retract(m: RiemannianManifold, x: ArrayOrContainerT, u: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.retract(m: ShapeSpace, x: ArrayOrContainerT, u: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.retract(m: CartesianEucledianSpace, x: ArrayOrContainerT, u: ArrayOrContainerT) → ArrayOrContainerT

Computes the retraction of the vector u from the tangent space at x.

A retraction \(R_x\) is a map from the tangent space \(T_x M\) back to the manifold \(M\). This map must satisfy

\[\begin{split}\begin{cases} R_x(0) = x, \\ D R_x(0)[u] = u, \end{cases}\end{split}\]

so that the zero vector is just mapped to the point itself and the directional derivative in the direction of the vector u is just the vector itself (i.e. the directional derivative at \(0\) is the identity map).

In the context of optimization, the retraction is used to generalize the update formula

\[x^{(k + 1)} = x^{(k)} + \alpha d^{(k)}\]

to

\[x^{(k + 1)} = R_{x^{(k)}}(\alpha d^{(k)}).\]

Parameters:

• x – a point on the manifold m;

• u – a vector in the tangent space of m at x.

Returns:

the retraction of the vector u from the tangent space at x.

pystopt.optimize.transport(m: RiemannianManifold, x: ArrayOrContainerT, y: ArrayOrContainerT, v: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.transport(m: ShapeSpace, x: ArrayOrContainerT, y: ArrayOrContainerT, v: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.transport(m: CartesianEucledianSpace, x: ArrayOrContainerT, y: ArrayOrContainerT, v: ArrayOrContainerT) → ArrayOrContainerT

Transports the vector v from the tangent space at x to y.

A transport \(T_{x, y}\) is a map from the tangent space \(T_x M\) to the tangent space \(T_y M\). In general, this should be consistent with the chosen retraction in retract().

In the context of optimization, the transport map is used to more advanced conjugate gradient algorithms, where we have a gradient and a descent direction. For example, the standard expression

\[d^{(k + 1)} = -g^{(k + 1)} + \beta d^{(k)}\]

to

\[d^{(k + 1)} = -g^{(k + 1)} + T_{x^{(k)}, x^{(k + 1)}}(d^{(k)}),\]

since \(g^{(k + 1)}\) and \(d^{(k)}\) are in two different tangent spaces and cannot be meaningfully added for curved manifolds.

Parameters:

• x – a point on the manifold m;

• y – a point on the manifold m;

• v – a vector in the tangent space of m at x.

Returns:

the vector v transported to the tangent space at y.

Explicit implementations of the RiemannianManifold class are given below.

class pystopt.optimize.CartesianEucledianSpace(vdot: Callable[[ArrayOrContainerT, ArrayOrContainerT], Array])

Bases: RiemannianManifold

A flat Eucledian manifold in a Cartesian coordinate system.

class pystopt.optimize.L2ShapeSpace(places: pytential.GeometryCollection, dofdesc: pytential.symbolic.primitives.DOFDescriptor, transported: bool)

Bases: ShapeSpace

A flat shape space using the \(L^2\) inner product.

This space only works on GeometryStates. The inner product on this space is given by

\[g(u, v) = \int_{\mathbb{S}^d} u \cdot v \,\mathrm{d}S,\]

where \(\mathbb{S}^d\) is the \(d\)-dimensional sphere. Note that this choice is not sufficient to describe a manifold (see [RingWirth2012]).

Besides that, it implements a non-trivial transport. Following [RingWirth2012], each point in the manifold is a closed surface and the tangent vectors are isomorphic to scalings of the normal vector field. Then, we set

\[T_{x, y}[\mathbb{v}] = (\mathbb{v} \cdot \mathbf{n}_x) \mathbf{n}_y,\]

where \(\mathbf{n}\) are the normal vectors on the surfaces represented by \(x\) and \(y\).

transported: bool

If False, the tangent vector transport to the normal direction on the given shape is not performed.

class pystopt.optimize.WeightedL2ShapeSpace(places: pytential.GeometryCollection, dofdesc: pytential.symbolic.primitives.DOFDescriptor, transported: bool, alpha: float)

Bases: L2ShapeSpace

A shape space using a weighted \(L^2\) inner product.

Unlike L2ShapeSpace, this manifold is no longer flat (see [RingWirth2012]). The inner product is given by

\[g(u, v) = \int_{\mathbb{S}^d} (1 + \alpha \kappa) u \cdot v \,\mathrm{d}S,\]

where \(\alpha\) is a real coefficient and \(\kappa\) is the sum of the principal curvatures. Setting \(\alpha = 0\) recovers the degenerate case of L2ShapeSpace.

alpha: float

Weight used for the curvature-based term.

Line Search

All line searches can be accessed programtically using

pystopt.optimize.get_line_search_from_name(name: str, **kwargs: Any) → LineSearch

Parameters:

• name – a string identifier for the line search method. The identifiers follow the naming convention "word-word-...", where each word is a lowercase part of the class name without "LineSearch". For example, for StabilizedBarzilaiBorweinLineSearch, the identifier is "stabilized-barzilai-borwein".

• kwargs – keyword arguments passed directly to the class constructor.

Line searches must be subclasses of LineSearch and implemented using functools.singledispatch() on line_search().

class pystopt.optimize.LineSearch(fun: Callable[[ArrayOrContainerT], float], jac: Callable[[ArrayOrContainerT], ArrayOrContainerT] | None = None)

Generic line search algorithm for steepest descent algorithms.

fun: Callable[[ArrayOrContainerT], float]

Functional that is being optimized.

jac: Callable[[ArrayOrContainerT], ArrayOrContainerT] | None = None

Gradient of fun.

pystopt.optimize.line_search(ls: LineSearch, m: RiemannianManifold, x: ArrayOrContainerT, f: float, g: ArrayOrContainerT, d: ArrayOrContainerT) → Tuple[float, ArrayOrContainerT]

pystopt.optimize.line_search(ls: FixedLineSearch, m: RiemannianManifold, x: ArrayOrContainerT, f: float, g: ArrayOrContainerT, d: ArrayOrContainerT) → float

pystopt.optimize.line_search(ls: BacktrackingLineSearch, m: RiemannianManifold, x: ArrayOrContainerT, f: float, g: ArrayOrContainerT, d: ArrayOrContainerT) → float

pystopt.optimize.line_search(ls: StabilizedBarzilaiBorweinLineSearch, m: RiemannianManifold, x: ArrayOrContainerT, f: float, g: ArrayOrContainerT, d: ArrayOrContainerT) → float

pystopt.optimize.line_search(ls: ModifiedBarzilaiBorweinLineSearch, m: RiemannianManifold, x: ArrayOrContainerT, f: float, g: ArrayOrContainerT, d: ArrayOrContainerT) → float

Approximation of the steepest descent step size.

The update is expected to be

\[x^{(k + 1)} = x^{(k)} + \alpha^{(k)} d^{(k)},\]

where \(\alpha^{(k)}\) is the step size approximated by this function. The descent direction d is expected to be a descent direction, i.e. \(\operatorname{dot}(d, g) < 0\). The manifold m is used to provide the specific meaning to all the operations above.

Parameters:

• x – previous point at which to evaluate functions.

• f – function value evaluated at x.

• g – gradient value evaluated at x.

• d – descent direction evaluated at x.

Returns:

an approximation of the step size and, optionally, the next iterate.

Concrete implementations include

class pystopt.optimize.FixedLineSearch(fun: Callable[[ArrayOrContainerT], float], jac: Callable[[ArrayOrContainerT], ArrayOrContainerT] | None = None, alpha: float | None = None)

Bases: LineSearch

Basic fixed step size without an actual line search.

If not provided, the step size is approximated by the following

\[\begin{split}\alpha = \begin{cases} \frac{\|x^{(0)}\|}{\|g^{(0)}\|}, & \quad \|x^{(0)}\| > 0, \\ 2 \frac{|f^{(0)}|}{\|g^{(0)}\|}, & \quad \text{otherwise}. \end{cases}\end{split}\]

alpha: float | None = None

A fixed given step size to return during optimization.

class pystopt.optimize.BacktrackingLineSearch(fun: Callable[[ArrayOrContainerT], float], jac: Callable[[ArrayOrContainerT], ArrayOrContainerT] | None = None, maxit: int = 64, max_alpha: float = 1.0, min_alpha: float = 1e-08, beta: float = 0.0001, factor: float = 0.85, nhistory: int = 5)

Bases: LineSearch

Implements the classic backtracking line search based on the Armijo-Goldstein condition.

The algorithm decreases the step size by factor while

\[f(\mathbf{x} + \alpha \mathbf{d}) - f(\mathbf{x}) > \alpha \beta (\mathbf{d} \cdot \mathbf{g}),\]

where \(\mathbf{g}\) is is the gradient and \(\mathbf{d}\) is the descent direction at \(\mathbf{x}\). If an \(\alpha\) is not found after maxit iterations, it just returns the last value.

The initial guess for the upper bound is given by

\[\tilde{\alpha}_{max} = \frac{1}{N_{history}} \sum_{k} \alpha^{(k)},\]

i.e. an average over the last nhistory terms. As it is an average, it is guaranteed to never exceed max_alpha.

nhistory: int = 5

Number of step sizes to keep as history for the average.

maxit: int = 64

Maximum number of iterations of line search to be performed.

max_alpha: float = 1.0

Upper bound for the step size.

min_alpha: float = 1e-08

Lower bound for the step size.

beta: float = 0.0001

Parameter used in the line search Armijo condition.

factor: float = 0.85

Decrease in step size if the Armijo condition is not satisfied.

class pystopt.optimize.GrowingBacktrackingLineSearch(fun: Callable[[ArrayOrContainerT], float], jac: Callable[[ArrayOrContainerT], ArrayOrContainerT] | None = None, maxit: int = 64, max_alpha: float = 1.0, min_alpha: float = 1e-08, beta: float = 0.0001, factor: float = 0.85, nhistory: int = 5, factor_max_alpha: float = 5.0)

Bases: BacktrackingLineSearch

A growing adaptive version of BacktrackingLineSearch.

The adaptation is in terms of the initial upper bound used to start the line search. it uses an upper bound given by

\[\tilde{\alpha}_{max} = \min \left( c \alpha_{history}, \frac{\alpha_{max}}{\min(1, \|d\|)} \right)\]

where \(\alpha_{history}\) is the upper bound from BacktrackingLineSearch.

factor_max_alpha: float = 5.0

A factor multiplying max_alpha to allow larger steps in some scenarios, but still control the magnitute.

class pystopt.optimize.StabilizedBarzilaiBorweinLineSearch(fun: ~typing.Callable[[~pystopt.optimize.linesearch.ArrayOrContainerT], float], jac: ~typing.Callable[[~pystopt.optimize.linesearch.ArrayOrContainerT], ~pystopt.optimize.linesearch.ArrayOrContainerT] | None = None, method: int = 2, first_alpha: float | None = None, max_alpha: float = 1.0, delta: float | None = None, delta_factor: float = 0.5, nhistory: int = 5, history: ~typing.List[float] = <factory>, require_monotonicity: bool = False)

Bases: LineSearch

Implements the step size estimate from [Burdakov2019].

[Burdakov2019] (1,2,3)

O. Burdakov, Y.-H. Dai, N. Huang, Stabilized Barzilai-Borwein Method, ARXIV.

According to the discussion in Section 4 of [Burdakov2019], the stabilization is more important when using method BB1. We use method BB2 with stabilization by default to provide best results at no additional cost.

method: int = 2

Version of the Barzilai-Borwein estimate, see [Burdakov2019] for explicit formulae. This is an integer value in {±1, ±2}, where negative values use the selected method but disable stabilization.

first_alpha: float | None = None

Guess for the initial step size. If not provided, it is computed from the available data.

max_alpha: float = 1.0

Maximum value of the step size, if the predicted step size grows out of control. This can happen as the optimum is reached.

delta: float | None = None

Parameter in stabilization formula. If not provided, it is estimated from the value of the first nhistory terms.

delta_factor: float = 0.5

Factor used in the estimate of delta based on the history (named \(c\) in [Burdakov2019]).

nhistory: int = 5

Number of items in the history of iterates that can be used to estimate delta.

history: List[float]

A list of the initial nhistory step sizes used to estimate delta.

require_monotonicity: bool = False

A flag to control if the method requires monotonicity in the function values. Since the step sizes are approximated, this does not always guarantee a decrease in the function values, so when True, we also perform a simple backtracking line search.

class pystopt.optimize.ModifiedBarzilaiBorweinLineSearch(fun: Callable[[ArrayOrContainerT], float], jac: Callable[[ArrayOrContainerT], ArrayOrContainerT] | None = None, max_alpha: float = inf, first_alpha: float | None = None, factor: float = 1.0)

Bases: LineSearch

Modified stabilized version of Barzilai-Borwein from [Wang2021].

[Wang2021]

TODO

This method is very similar to StabilizedBarzilaiBorweinLineSearch. The main difference is a simplification that removes the \(Delta\) parameter and simply uses the previous step for an estimate.

In general, this method is expected to work equally well.

max_alpha: float = inf

Maximum value of the step size, if the predicted step size grows out of control. This can happen as the optimum is reached.

first_alpha: float | None = None

Step size at the first step. If not given, it is approximated from the available data.

factor: float = 1.0

A factor that can is multiplied into the estimated step size. Can be used to control the magnitude of the step.

In the context of shape optimization, the following methods can be used to get an estimate of the maximum allowable step size.

pystopt.optimize.estimate_line_search_step_dogan(ambient_dim: int, g: Expression | ndarray, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

An estimate for the CG step size based on Section 4.1 in [Dogan2007].

[Dogan2007]

G. Doǧan, P. Morin, R. H. Nochetto, M. Verani, Discrete Gradient Flows for Shape Optimization and Applications, Computer Methods in Applied Mechanics and Engineering, Vol. 196, pp. 3898–3914, 2007, DOI.

Parameters:

g – shape gradient or only its normal component.

Returns:

a symbolic expression for the estimate of the step size.

pystopt.optimize.estimate_line_search_step_feppon(ambient_dim: int, g: Expression | ndarray, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

An estimate for the CG step size based on Section 5.3.2 in [Feppon2020].

[Feppon2020]

F. Feppon, G. Allaire, C. Dapogny, Null Space Gradient Flows for Constrained Optimization With Applications to Shape Optimization, ESAIM: Control, Optimisation and Calculus of Variations, Vol. 26, pp. 90, 2020, DOI.

Parameters:

g – shape gradient or only its normal component.

Returns:

a symbolic expression for the estimate of the step size.

Descent Direction

Note

Some descent directions should exclusively be used with a specific implementation of pystopt.optimize.LineSearch.

All descent directions can be obtained programmatically using

pystopt.optimize.get_descent_direction_from_name(name: str, **kwargs) → DescentDirection

Parameters:

• name – a string identifier for the descent direction method. The identifiers follow the naming convention "word-word-..."`, where each word is the lowercase part of the class name. For example, FletcherReevesDirection is identified by "fletcher-reeves".

• kwargs – keyword arguments passed directly to the class constructor.

Descent directions must be subclasses of

class pystopt.optimize.DescentDirection

The descent direction is computed as

\[\mathbf{d}^{(k + 1)} = -\mathbf{g}^{(k + 1)} + \beta^{(k)} \mathbf{d}^{(k)}\]

where the CG update parameter \(\beta\) depends on

\[\beta^{(k)} \equiv \beta^{(k)}( \mathbf{g}^{(k + 1)}, \mathbf{g}^{(k)}, \mathbf{d}^{(k)} )\]

and implement the following function through functools.singledispatch()

pystopt.optimize.descent_direction(ds: DescentDirection, m: RiemannianManifold, x: ArrayOrContainerT, g: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.descent_direction(ds: SteepestDescentDirection, m: RiemannianManifold, x: ArrayOrContainerT, g: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.descent_direction(ds: HestenesStiefelDirection, m: RiemannianManifold, x: ArrayOrContainerT, g: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.descent_direction(ds: FletcherReevesDirection, m: RiemannianManifold, x: ArrayOrContainerT, g: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.descent_direction(ds: PolakRibiereDirection, m: RiemannianManifold, x: ArrayOrContainerT, g: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.descent_direction(ds: FletcherDirection, m: RiemannianManifold, x: ArrayOrContainerT, g: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.descent_direction(ds: LiuStoreyDirection, m: RiemannianManifold, x: ArrayOrContainerT, g: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.descent_direction(ds: DaiYuanDirection, m: RiemannianManifold, x: ArrayOrContainerT, g: ArrayOrContainerT) → ArrayOrContainerT

pystopt.optimize.descent_direction(ds: DaiYuanDirection, m: RiemannianManifold, x: ArrayOrContainerT, g: ArrayOrContainerT) → ArrayOrContainerT

Constructs \(\mathbf{d}^{(k + 1)}\) from \(\mathbf{g}^{(k + 1)}\). and the datum at the previous step in DescentDirection.

Concrete implementations include

class pystopt.optimize.SteepestDescentDirection

Bases: DescentDirection

class pystopt.optimize.HestenesStiefelDirection

Bases: DescentDirection

class pystopt.optimize.FletcherReevesDirection

Bases: DescentDirection

class pystopt.optimize.PolakRibiereDirection

Bases: DescentDirection

class pystopt.optimize.FletcherDirection

Bases: DescentDirection

class pystopt.optimize.LiuStoreyDirection

Bases: DescentDirection

class pystopt.optimize.DaiYuanDirection

Bases: DescentDirection

class pystopt.optimize.HagerZhangDirection

Bases: DescentDirection

Helpers

pystopt.optimize.flatten_to_numpy(actx: ArrayContext | None, ary: ArrayOrContainer) → ndarray

Flatten and transfer an array container to numpy.

This is a composition of to_numpy() and flatten().

pystopt.optimize.unflatten_from_numpy_like(actx: ArrayContext | None, ary: ndarray, prototype: ArrayOrContainerT) → ArrayOrContainerT

A reverse operator to flatten_to_numpy().

pystopt.optimize.get_vdot_for_container(x: ArrayOrContainerT, actx: ArrayContext | None = None) → Callable[[ArrayOrContainerT, ArrayOrContainerT], float]

Uses heuristics to get a vdot() callable for the type of x.

pystopt.optimize.get_norm_for_container(x: ArrayOrContainerT, actx: ArrayContext | None = None) → Callable[[ArrayOrContainerT], float]

Uses heuristics to get a norm() for the type of x.

class pystopt.optimize.CGResult(x: ndarray, success: bool, status: int, message: str, fun: float, jac: float, nfev: int, njev: int, nit: int)

Bases: object

Based on scipy.optimize.OptimizeResult.

x: ndarray

Solution of the optimization.

success: bool

Flag to denote a successful exit.

status: int

Termination status of the optimize.

message: str

Description of the termination status in status.

fun: float

Function value at the end of the optimization.

jac: float

Norm of the gradient (Jacobian) at the end of the optimization.

nfev: int

Number of function evaluations.

njev: int

Number of gradient (Jacobian) evaluation.

nit: int

Number of iterations performed by the optimizer.

pretty() → str

A fancier stringify of the results.

class pystopt.optimize.CGCallbackInfo(it: int, alpha: float, x: ndarray, f: float, g: ndarray, d: ndarray)

Bases: object

it: int

Current iteration.

alpha: float

Current step size, such that \(x = x_p + \alpha d\), where \(x_p\) is the previous iterate.

x: ndarray

Solution at current iteration.

f: float

Function value at current iteration.

g: ndarray

Gradient at the current iteration.

d: ndarray

Descent direction at the current iteration. This will usually not be the same as the gradient and can be used for debugging.

class pystopt.optimize.CGLogCallback(norm: Callable[[Any], float] | None, log: Callable[[...], Any] | None)

Bases: LogCallback

__call__(*args: Any, **kwargs: Any) → int

Parameters:

info – an instance of CGCallbackInfo.

class pystopt.optimize.CGHistoryCallback(norm: Callable[[Any], float] | None)

Bases: HistoryCallback

Store values from CGCallbackInfo from all iterations.

The values stored as the function values, the gradient norm and the descent direction norm.

__call__(info, *args, **kwargs)

Parameters:

info – an instance of CGCallbackInfo.

f: List[float]

History of function values.

g: List[float]

History of \(\ell^2\) norms of the gradient.

d: List[float]

History of \(\ell^2\) norms of the descent direction.

alpha: List[float]

History of step sizes.