Gradient descent

Manopt.gradient_descentFunction
gradient_descent(M, f, grad_f, p=rand(M); kwargs...)
gradient_descent(M, gradient_objective, p=rand(M); kwargs...)
gradient_descent!(M, f, grad_f, p; kwargs...)
gradient_descent!(M, gradient_objective, p; kwargs...)

perform the gradient descent algorithm

\[p_{k+1} = \operatorname{retr}_{p_k}\bigl( s_k\operatorname{grad}f(p_k) \bigr), \qquad k=0,1,…\]

where $s_k > 0$ denotes a step size.

The algorithm can be performed in-place of p.

Input

  • M::AbstractManifold: a Riemannian manifold $\mathcal{M}$
  • f: a cost function $f: \mathcal{M}→ ℝ$ implemented as (M, p) -> v
  • grad_f: the (Riemannian) gradient $\operatorname{grad}f: \mathcal{M} → T_{p}\mathcal{M}$ of f as a function (M, p) -> X or a function (M, X, p) -> X computing X in-place
  • p::P: a point on the manifold $\mathcal{M}$

Alternatively to f and grad_f you can provide the corresponding AbstractManifoldFirstOrderObjective gradient_objective directly.

Keyword arguments

All other keyword arguments are passed to decorate_state! for state decorators or decorate_objective! for objective decorators, respectively.

If you provide the ManifoldFirstOrderObjective directly, the evaluation= keyword is ignored. The decorations are still applied to the objective.

If you activate tutorial mode (cf. is_tutorial_mode), this solver provides additional debug warnings.

Output

The obtained approximate minimizer $p^*$. To obtain the whole final state of the solver, see get_solver_return for details, especially the return_state= keyword.

source
Manopt.gradient_descent!Function
gradient_descent(M, f, grad_f, p=rand(M); kwargs...)
gradient_descent(M, gradient_objective, p=rand(M); kwargs...)
gradient_descent!(M, f, grad_f, p; kwargs...)
gradient_descent!(M, gradient_objective, p; kwargs...)

perform the gradient descent algorithm

\[p_{k+1} = \operatorname{retr}_{p_k}\bigl( s_k\operatorname{grad}f(p_k) \bigr), \qquad k=0,1,…\]

where $s_k > 0$ denotes a step size.

The algorithm can be performed in-place of p.

Input

  • M::AbstractManifold: a Riemannian manifold $\mathcal{M}$
  • f: a cost function $f: \mathcal{M}→ ℝ$ implemented as (M, p) -> v
  • grad_f: the (Riemannian) gradient $\operatorname{grad}f: \mathcal{M} → T_{p}\mathcal{M}$ of f as a function (M, p) -> X or a function (M, X, p) -> X computing X in-place
  • p::P: a point on the manifold $\mathcal{M}$

Alternatively to f and grad_f you can provide the corresponding AbstractManifoldFirstOrderObjective gradient_objective directly.

Keyword arguments

All other keyword arguments are passed to decorate_state! for state decorators or decorate_objective! for objective decorators, respectively.

If you provide the ManifoldFirstOrderObjective directly, the evaluation= keyword is ignored. The decorations are still applied to the objective.

If you activate tutorial mode (cf. is_tutorial_mode), this solver provides additional debug warnings.

Output

The obtained approximate minimizer $p^*$. To obtain the whole final state of the solver, see get_solver_return for details, especially the return_state= keyword.

source

State

Manopt.GradientDescentStateType
GradientDescentState{P,T} <: AbstractGradientSolverState

Describes the state of a gradient based descent algorithm.

Fields

  • callbacks::D: provided callback functions given as a dictionary with symbols as keys
  • p::P: a point on the manifold $\mathcal{M}$ storing the current iterate
  • X::T: a tangent vector at the point $p$ on the manifold $\mathcal{M}$ storing the gradient at the current iterate
  • stop::StoppingCriterion: a functor indicating that the stopping criterion is fulfilled
  • stepsize::Stepsize: a functor inheriting from Stepsize to determine a step size
  • direction::DirectionUpdateRule : a processor to handle the obtained gradient and compute a direction to “walk into”.
  • retraction_method::AbstractRetractionMethod: a retraction $\operatorname{retr}$ to use, see the section on retractions

Constructor

GradientDescentState(M::AbstractManifold; kwargs...)

Initialize the gradient descent solver state, where

Input

Keyword arguments

See also

gradient_descent

source

Direction update rules

A field of the options is the direction, a DirectionUpdateRule, which by default IdentityUpdateRule just evaluates the gradient. The following ones are furthermore available

Manopt.AverageGradientMethod
AverageGradient(; kwargs...)
AverageGradient(M::AbstractManifold; kwargs...)

Add an average of gradients to a gradient processor. A set of previous directions (from the inner processor) and the last iterate are stored, average is taken after vector transporting them to the current iterates tangent space.

Input

Keyword arguments

Info

This function generates a ManifoldDefaultsFactory for AverageGradientRule. For default values, that depend on the manifold, this factory postpones the construction until the manifold from for example a corresponding AbstractManoptSolverState is available.

source
Manopt.MomentumGradientMethod
MomentumGradient(args...; kwargs...)

Append a momentum to a gradient processor.

The last direction and last iterate are stored and the new is composed as $η_i = m*η_{i-1}' - s d_i$, where $sd_i$ is the current (inner) direction and $η_{i-1}'$ is the vector transported last direction multiplied by momentum $m$.

Input

  • M (optional)

Keyword arguments

Info

This function generates a ManifoldDefaultsFactory for MomentumGradientRule. For default values, that depend on the manifold, this factory postpones the construction until the manifold from for example a corresponding AbstractManoptSolverState is available.

source
Manopt.NesterovMethod
Nesterov(; kwargs...)
Nesterov(M::AbstractManifold; kwargs...)

Assume $f$ is $L$-Lipschitz and $μ$-strongly convex. Given

  • a step size $h_k<\frac{1}{L}$ (from the GradientDescentState
  • a shrinkage parameter $β_k$
  • and a current iterate $p_k$
  • as well as the interim values $γ_k$ and $v_k$ from the previous iterate.

This compute a Nesterov type update using the following steps, see [ZS18]

  1. Compute the positive root $α_k∈(0,1)$ of $α^2 = h_k\bigl((1-α_k)γ_k+α_k μ\bigr)$.
  2. Set $\barγ_k+1 = (1-α_k)γ_k + α_kμ$
  3. $y_k = \operatorname{retr}_{p_k}\Bigl(\frac{α_kγ_k}{γ_k + α_kμ}\operatorname{retr}^{-1}_{p_k}v_k \Bigr)$
  4. $x_{k+1} = \operatorname{retr}_{y_k}(-h_k \operatorname{grad}f(y_k))$
  5. $v_{k+1} = \operatorname{retr}_{y_k}\Bigl(\frac{(1-α_k)γ_k}{\barγ_k}\operatorname{retr}_{y_k}^{-1}(v_k) - \frac{α_k}{\barγ_{k+1}}\operatorname{grad}f(y_k) \Bigr)$
  6. $γ_{k+1} = \frac{1}{1+β_k}\barγ_{k+1}$

Then the direction from $p_k$ to $p_k+1$ by $d = \operatorname{retr}^{-1}_{p_k}p_{k+1}$ is returned.

Input

Keyword arguments

Info

This function generates a ManifoldDefaultsFactory for NesterovRule. For default values, that depend on the manifold, this factory postpones the construction until the manifold from for example a corresponding AbstractManoptSolverState is available.

source
Manopt.PreconditionedDirectionMethod
PreconditionedDirection(preconditioner; kwargs...)
PreconditionedDirection(M::AbstractManifold, preconditioner; kwargs...)

Add a preconditioner to a gradient processor following the motivation for optimization, as a linear invertible map $P: T_{p}\mathcal{M} → T_{p}\mathcal{M}$ that usually should be

  • symmetric: $⟨X, P(Y)⟩ = ⟨P(X), Y⟩$
  • positive definite $⟨X, P(X)⟩ > 0$ for $X$ not the zero-vector

The gradient is then preconditioned as $P(X)$, where $X$ is either the gradient of the objective or the result of a previous (internally stored) gradient processor.

For example if you provide as the preconditioner the inverse of the Hessian $\operatorname{Hess}^{-1} f$, you turn a gradient descent into a Newton method.

Arguments

  • M::AbstractManifold: a Riemannian manifold $\mathcal{M}$ (optional)
  • preconditioner: preconditioner function, either as a (M, p, X) -> Y allocating or (M, Y, p, X) -> Y mutating function

Keyword arguments

Info

This function generates a ManifoldDefaultsFactory for PreconditionedDirectionRule. For default values, that depend on the manifold, this factory postpones the construction until the manifold from for example a corresponding AbstractManoptSolverState is available.

source

Internal functions

Manopt.AverageGradientRuleType
AverageGradientRule <: DirectionUpdateRule

Add an average of gradients to a gradient processor.

A set of previous directions (from the inner processor) and the last iterate are stored. The average is taken after vector transporting them to the current iterates tangent space.

Fields

Constructors

AverageGradientRule(    M::AbstractManifold;    p::P=rand(M);    n::Int=10    direction::Union{<:DirectionUpdateRule,ManifoldDefaultsFactory}=IdentityUpdateRule(),    gradients = fill(zero_vector(p.M, o.x),n),    last_iterate = deepcopy(x0),    vector_transport_method = default_vector_transport_method(M, typeof(p)))AverageGradientRule(M::AbstractManifold, p; kwargs...)

Add average to a gradient problem, where

source
Manopt.IdentityUpdateRuleType
IdentityUpdateRule <: DirectionUpdateRule

The default gradient direction update is the identity, usually it just evaluates the gradient.

You can also use Gradient() to create the corresponding factory, though this only delays this parameter-free instantiation to later.

source
Manopt.MomentumGradientRuleType
MomentumGradientRule <: DirectionUpdateRule

Store the necessary information to compute the MomentumGradient direction update.

Fields

Constructors

MomentumGradientRule(M::AbstractManifold; kwargs...)MomentumGradientRule(M::AbstractManifold, p; kwargs...)

Initialize a momentum gradient rule to s, where p and X are memory for interim values.

Keyword arguments

See also

MomentumGradient

source
Manopt.NesterovRuleType
NesterovRule <: DirectionUpdateRule

Compute a Nesterov inspired direction update rule. See Nesterov for details

Fields

Constructor

NesterovRule(M::AbstractManifold; kwargs...)NesterovRule(M::AbstractManifold, p; kwargs...)

Keyword arguments

See also

Nesterov

source
Manopt.PreconditionedDirectionRuleType
PreconditionedDirectionRule <: DirectionUpdateRule

Add a preconditioning as gradient processor, see PreconditionedDirection for more mathematical background.

Fields

  • direction: internal DirectionUpdateRule to determine directions to apply the preconditioning to
  • preconditioner: the preconditioner function

Constructors

PreconditionedDirectionRule(    M::AbstractManifold,    preconditioner;    direction::Union{<:DirectionUpdateRule,ManifoldDefaultsFactory}=IdentityUpdateRule(),    evaluation::AbstractEvaluationType=AllocatingEvaluation())

Add preconditioning to a gradient problem.

Input

  • M::AbstractManifold: a Riemannian manifold $\mathcal{M}$
  • preconditioner: preconditioner function, either as a (M, Y, p, X) -> Y mutating function

Keyword arguments

source

Technical details

The gradient_descent solver requires the following functions of a manifold to be available

  • A retract!(M, q, p, X); it is recommended to set the default_retraction_method to a favourite retraction. If this default is set, a retraction_method= does not have to be specified.
  • By default gradient descent uses ArmijoLinesearch which requires max_stepsize(M) to be set and an implementation of inner(M, p, X).
  • By default the stopping criterion uses the norm as well, to stop when the norm of the gradient is small, but if you implemented inner, the norm is provided already.
  • By default the tangent vector storing the gradient is initialized calling zero_vector(M,p).

Literature

[ABG06]
P.-A. Absil, C. Baker and K. Gallivan. Trust-Region Methods on Riemannian Manifolds. Foundations of Computational Mathematics 7, 303–330 (2006).
[AMS08]
P.-A. Absil, R. Mahony and R. Sepulchre. Optimization Algorithms on Matrix Manifolds (Princeton University Press, 2008), available online at press.princeton.edu/chapters/absil/.
[AOT22]
S. Adachi, T. Okuno and A. Takeda. Riemannian Levenberg-Marquardt Method with Global and Local Convergence Properties. ArXiv Preprint (2022).
[ABBC20]
N. Agarwal, N. Boumal, B. Bullins and C. Cartis. Adaptive regularization with cubics on manifolds. Mathematical Programming (2020).
[AdCOdO20]
Y. T. Almeida, J. X. da Cruz Neto, P. R. Oliveira and J. C. de Oliveira Souza. A modified proximal point method for DC functions on Hadamard manifolds. Computational Optimization and Applications 76, 649–673 (2020).
[BB26]
[BBP26]
M. Baran, R. Bergmann and P. Przybysz. A Riemannian quasi-Newton algorithm for optimization with Euclidean bounds (2026), arXiv:2605.10573 [math.OC].
[Bac14]
M. Bačák. Computing medians and means in Hadamard spaces. SIAM Journal on Optimization 24, 1542–1566 (2014), arXiv:1210.2145.
[Bea72]
E. M. Beale. A derivation of conjugate gradients. In: Numerical methods for nonlinear optimization, edited by F. A. Lootsma (Academic Press, London, London, 1972); pp. 39–43.
[BFNZ25]
R. Bergmann, O. P. Ferreira, S. Z. Németh and J. Zhu. On projection mappings and the gradient projection method on hyperbolic space forms. Preprint, in preparation (2025).
[BFSS24]
R. Bergmann, O. P. Ferreira, E. M. Santos and J. C. Souza. The difference of convex algorithm on Hadamard manifolds. Journal of Optimization Theory and Applications (2024).
[BH19]
R. Bergmann and R. Herzog. Intrinsic formulation of KKT conditions and constraint qualifications on smooth manifolds. SIAM Journal on Optimization 29, 2423–2444 (2019), arXiv:1804.06214.
[BHJ24]
R. Bergmann, R. Herzog and H. Jasa. The Riemannian Convex Bundle Method, preprint (2024), arXiv:2402.13670.
[BHS+21]
R. Bergmann, R. Herzog, M. Silva Louzeiro, D. Tenbrinck and J. Vidal-Núñez. Fenchel duality theory and a primal-dual algorithm on Riemannian manifolds. Foundations of Computational Mathematics 21, 1465–1504 (2021), arXiv:1908.02022.
[BJJP25a]
R. Bergmann, H. Jasa, P. John and M. Pfeffer. The Intrinsic Riemannian Proximal Gradient Method for Convex Optimization, preprint (2025), arXiv:2507.16055.
[BJJP25b]
R. Bergmann, H. Jasa, P. John and M. Pfeffer. The Intrinsic Riemannian Proximal Gradient Method for Nonconvex Optimization, preprint (2025), arXiv:2506.09775.
[BPS16]
R. Bergmann, J. Persch and G. Steidl. A parallel Douglas Rachford algorithm for minimizing ROF-like functionals on images with values in symmetric Hadamard manifolds. SIAM Journal on Imaging Sciences 9, 901–937 (2016), arXiv:1512.02814.
[Ber15]
D. P. Bertsekas. Convex Optimization Algorithms (Athena Scientific, 2015); p. 576.
[BIA10]
P. B. Borckmans, M. Ishteva and P.-A. Absil. A Modified Particle Swarm Optimization Algorithm for the Best Low Multilinear Rank Approximation of Higher-Order Tensors. In: 7th International Conference on Swarm INtelligence (Springer Berlin Heidelberg, 2010); pp. 13–23.
[Bou23]
[BLNZ95]
R. H. Byrd, P. Lu, J. Nocedal and C. Zhu. A Limited Memory Algorithm for Bound Constrained Optimization. SIAM Journal on Scientific Computing 16, 1190–1208 (1995). Publisher: Society for Industrial and Applied Mathematics.
[BNS94]
R. H. Byrd, J. Nocedal and R. B. Schnabel. Representations of quasi-Newton matrices and their use in limited memory methods. Mathematical Programming 63, 129–156 (1994).
[Car92]
M. P. do Carmo. Riemannian Geometry. Mathematics: Theory & Applications (Birkhäuser Boston, Inc., Boston, MA, 1992); p. xiv+300.
[CP11]
A. Chambolle and T. Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision 40, 120–145 (2011).
[CFFS10]
[CGT00]
A. R. Conn, N. I. Gould and P. L. Toint. Trust Region Methods (Society for Industrial and Applied Mathematics, 2000).
[DY99]
Y. H. Dai and Y. Yuan. A Nonlinear Conjugate Gradient Method with a Strong Global Convergence Property. SIAM Journal on Optimization 10, 177–182 (1999).
[DL21]
W. Diepeveen and J. Lellmann. An Inexact Semismooth Newton Method on Riemannian Manifolds with Application to Duality-Based Total Variation Denoising. SIAM Journal on Imaging Sciences 14, 1565–1600 (2021), arXiv:2102.10309.
[DSN24]
D. Dodd, L. Sharrock and C. Nemeth. Learning-rate-free stochastic optimization over Riemannian manifolds, arXiv preprint arXiv:2406.02296 (2024).
[Dre07]
D. W. Dreisigmeyer. Direct Search Algorithms over Riemannian Manifolds (Optimization Online, 2007).
[ETTZ96]
A. S. El-Bakry, R. A. Tapia, T. Tsuchiya and Y. Zhang. On the formulation and theory of the Newton interior-point method for nonlinear programming. Journal of Optimization Theory and Applications 89, 507–541 (1996).
[Fan06]
J. Fan. Convergence Rate of The Trust Region Method for Nonlinear Equations Under Local Error Bound Condition. Computational Optimization and Applications 34, 215–227 (2006).
[FHS+21]
S. Feng, W. Huang, L. Song, S. Ying and T. Zeng. Proximal gradient method for nonconvex and nonsmooth optimization on Hadamard manifolds. Optimization Letters 16, 2277–2297 (2021).
[FO98]
O. Ferreira and P. R. Oliveira. Subgradient algorithm on Riemannian manifolds. Journal of Optimization Theory and Applications 97, 93–104 (1998).
[FO02]
O. Ferreira and P. R. Oliveira. Proximal point algorithm on Riemannian manifolds. Optimization. A Journal of Mathematical Programming and Operations Research 51, 257–270 (2002).
[Fle87]
R. Fletcher. Practical Methods of Optimization. 2 Edition, A Wiley-Interscience Publication (John Wiley & Sons Ltd., 1987).
[FR64]
R. Fletcher and C. M. Reeves. Function minimization by conjugate gradients. The Computer Journal 7, 149–154 (1964).
[GILO26]
J. Glaubitz, A. Iske, J. Lampert and P. Öffner. Why summation by parts is not enough (02 2026), arXiv:2602.10786.
[GS23]
[HZ06a]
W. W. Hager and H. Zhang. A survey of nonlinear conjugate gradient methods. Pacific Journal of Optimization 2, 35–58 (2006).
[HZ06b]
W. W. Hager and H. Zhang. Algorithm 851: CG_DESCENT, a conjugate gradient method with guaranteed descent. ACM Transactions on Mathematical Software 32, 113–137 (2006).
[HZ05]
W. W. Hager and H. Zhang. A New Conjugate Gradient Method with Guaranteed Descent and an Efficient Line Search. SIAM Journal on Optimization 16, 170–192 (2005).
[Hag89]
W. W. Hager. A derivative-based bracketing scheme for univariate minimization and the conjugate gradient method. Computers & Mathematics with Applications 18, 779–795 (1989).
[Han23]
N. Hansen. The CMA Evolution Strategy: A Tutorial. ArXiv Preprint (2023).
[HS52]
M. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards 49, 409 (1952).
[HNP23]
N. Hoseini Monjezi, S. Nobakhtian and M. R. Pouryayevali. A proximal bundle algorithm for nonsmooth optimization on Riemannian manifolds. IMA Journal of Numerical Analysis 43, 293–325 (2023).
[HU17]
S. Hosseini and A. Uschmajew. A Riemannian Gradient Sampling Algorithm for Nonsmooth Optimization on Manifolds. SIAM Journal on Optimization 27, 173–189 (2017).
[Hua14]
W. Huang. Optimization algorithms on Riemannian manifolds with applications. Ph.D. Thesis, Florida State University (2014).
[HAG18]
W. Huang, P.-A. Absil and K. A. Gallivan. A Riemannian BFGS method without differentiated retraction for nonconvex optimization problems. SIAM Journal on Optimization 28, 470–495 (2018).
[HGA15]
W. Huang, K. A. Gallivan and P.-A. Absil. A Broyden class of quasi-Newton methods for Riemannian optimization. SIAM Journal on Optimization 25, 1660–1685 (2015).
[IP17]
B. Iannazzo and M. Porcelli. The Riemannian Barzilai–Borwein method with nonmonotone line search and the matrix geometric mean computation. IMA Journal of Numerical Analysis 38, 495–517 (2017).
[Kar77]
H. Karcher. Riemannian center of mass and mollifier smoothing. Communications on Pure and Applied Mathematics 30, 509–541 (1977).
[LY24]
Z. Lai and A. Yoshise. Riemannian Interior Point Methods for Constrained Optimization on Manifolds. Journal of Optimization Theory and Applications 201, 433–469 (2024), arXiv:2203.09762.
[LB19]
C. Liu and N. Boumal. Simple algorithms for optimization on Riemannian manifolds with constraints. Applied Mathematics & Optimization (2019), arXiv:1091.10000.
[LS91]
Y. Liu and C. Storey. Efficient generalized conjugate gradient algorithms, part 1: Theory. Journal of Optimization Theory and Applications 69, 129–137 (1991).
[Lue72]
D. G. Luenberger. The gradient projection method along geodesics. Management Science 18, 620–631 (1972).
[Ngu23]
D. Nguyen. Operator-Valued Formulas for Riemannian Gradient and Hessian and Families of Tractable Metrics in Riemannian Optimization. Journal of Optimization Theory and Applications 198, 135–164 (2023), arXiv:2009.10159.
[NW06]
J. Nocedal and S. J. Wright. Numerical Optimization. 2 Edition (Springer, New York, 2006).
[Pee93]
R. Peeters. On a Riemannian version of the Levenberg-Marquardt algorithm. Serie Research Memoranda 0011 (VU University Amsterdam, Faculty of Economics, Business Administration and Econometrics, 1993).
[PR69]
E. Polak and G. Ribière. Note sur la convergence de méthodes de directions conjuguées. Revue française d’informatique et de recherche opérationnelle 3, 35–43 (1969).
[Pow77]
M. J. Powell. Restart procedures for the conjugate gradient method. Mathematical Programming 12, 241–254 (1977).
[SI20]
H. Sakai and H. Iiduka. Hybrid Riemannian conjugate gradient methods with global convergence properties. Computational Optimization and Applications 77, 811–830 (2020).
[SI21]
H. Sakai and H. Iiduka. Sufficient Descent Riemannian Conjugate Gradient Methods. Journal of Optimization Theory and Applications 190, 130–150 (2021).
[SO15]
J. C. Souza and P. R. Oliveira. A proximal point algorithm for DC functions on Hadamard manifolds. Journal of Global Optimization 63, 797–810 (2015).
[TMHF00]
B. Triggs, P. F. McLauchlan, R. I. Hartley and A. W. Fitzgibbon. Bundle Adjustment — A Modern Synthesis. In: Vision Algorithms: Theory and Practice, Vol. 1883, edited by G. Goos, J. Hartmanis, J. Van Leeuwen, B. Triggs, A. Zisserman and R. Szeliski (Springer Berlin Heidelberg, Berlin, Heidelberg, 2000); pp. 298–372.
[WS22]
M. Weber and S. Sra. Riemannian Optimization via Frank-Wolfe Methods. Mathematical Programming 199, 525–556 (2022).
[WBS25]
L. Weigl, R. Bergmann and A. Schiela. Newton's method for nonlinear mappings into vector bundles Part II: Application to variational problems, preprint (2025), arXiv:2507.13836.
[WS24]
L. Weigl and A. Schiela. Newton's method for nonlinear mappings into vector bundles, preprint (2024), arXiv:2404.04073.
[Yua15]
Y.-x. Yuan. Recent advances in trust region algorithms. Mathematical Programming 151, 249–281 (2015).
[ZS18]
H. Zhang and S. Sra. Towards Riemannian accelerated gradient methods, preprint (2018), arXiv:1806.02812.
[ZBLN97]
C. Zhu, R. H. Byrd, P. Lu and J. Nocedal. Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization. ACM Trans. Math. Softw. 23, 550–560 (1997).