From f34d81eed3abab910e76c1b09b318238cce6c571 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Mon, 2 Nov 2026 07:24:36 -0800 Subject: [PATCH] Add multivariate Newton and Cuyt's multivariate Halley root finders One LU of J serves both solves; F''[a, a] comes from a user closure. Line search on |F|^2. Closes #1616 Co-Authored-By: Claude Opus 5.5 --- doc/roots/multivariate_roots.qbk | 164 ++++++++ doc/roots/roots_overview.qbk | 4 +- example/Jamfile.v2 | 1 + example/multivariate_roots_ellipse.cpp | 103 +++++ .../boost/math/tools/multivariate_roots.hpp | 357 ++++++++++++++++ module/math.cppm | 1 + nanobind/src/roots.cpp | 66 ++- nanobind/tests/test_roots.py | 28 ++ .../multivariate_roots_performance.cpp | 113 +++++ test/Jamfile.v2 | 2 + .../test_tools_multivariate_roots_incl.cpp | 23 ++ test/test_multivariate_roots.cpp | 390 ++++++++++++++++++ 12 files changed, 1248 insertions(+), 4 deletions(-) create mode 100644 doc/roots/multivariate_roots.qbk create mode 100644 example/multivariate_roots_ellipse.cpp create mode 100644 include/boost/math/tools/multivariate_roots.hpp create mode 100644 reporting/performance/multivariate_roots_performance.cpp create mode 100644 test/compile_test/test_tools_multivariate_roots_incl.cpp create mode 100644 test/test_multivariate_roots.cpp diff --git a/doc/roots/multivariate_roots.qbk b/doc/roots/multivariate_roots.qbk new file mode 100644 index 0000000000..7a258b6242 --- /dev/null +++ b/doc/roots/multivariate_roots.qbk @@ -0,0 +1,164 @@ +[section:multivariate_roots Multivariate Newton and Halley] + +[import ../../example/multivariate_roots_ellipse.cpp] + +[h4 Synopsis] + +`` +#include +`` + + namespace boost { namespace math { namespace tools { + + template + Container multivariate_newton(F f, Container x, int digits); + + template + Container multivariate_newton(F f, Container x, int digits, std::uintmax_t& max_iter); + + template + Container multivariate_halley(F f, Container x, int digits); + + template + Container multivariate_halley(F f, Container x, int digits, std::uintmax_t& max_iter); + + }}} // namespaces boost::math::tools + +[h4 Description] + +These functions find a root of a system of /n/ equations in /n/ unknowns, [*F](/x/) = 0 with /x/ in ℝ[super /n/], +starting from the guess `x`. +They are the multivariate analogues of `newton_raphson_iterate` and `halley_iterate` of __root_finding_with_derivatives. + +* `multivariate_newton` solves /J/(/x/)/a/ = −[*F](/x/) for the Newton step /a/, and converges quadratically. +* `multivariate_halley` is the multivariate Halley method of Cuyt, which converges cubically. +With the Newton step /a/ and /b/ = /J/[super −1][*F]″(/x/)\[/a/, /a/\], it updates each component as + +[expression /x/[sub /i/] ← /x/[sub /i/] + /a/[sub /i/][super 2] \/ (/a/[sub /i/] + /b/[sub /i/]\/2).] + +This is the abstract Padé approximant of order (1, 1) to the inverse of [*F], with the vector operations taken componentwise: +when /n/ = 1 it is Halley's method. +One LU factorization of /J/ serves both linear solves, and the method needs [*F]″ only along the direction /a/, +so the /n/[super 3] entries of the second derivative tensor are never formed. + +[variablelist Parameters +[[F f] [A callable that takes a `Container` /x/ and returns a pair or tuple. +For Newton it holds [*F](/x/) and the Jacobian /J/(/x/), /J/[sub /ij/] = ∂/F/[sub /i/]\/∂/x/[sub /j/]. +For Halley it holds a third element, a callable `d2` with `d2(v)` = [*F]″(/x/)\[/v/, /v/\], +the second derivative of /t/ ↦ [*F](/x/ + /t v/) at /t/ = 0; Newton ignores a third element, so one callable serves both. +[*F](/x/) and `d2(v)` may be any containers with `size()` and `operator[]`. +/J/ may be a container of rows (`std::vector>`, `std::array, N>`, ...), +or a row-major array of /n/[super 2] entries (`std::vector`, `std::array`, ...).]] +[[Container x] [The initial guess, a `std::vector`, `std::array`, or similar random access container +of a floating point type, which may be a multiprecision type.]] +[[int digits] [The desired number of binary digits of precision.]] +[[uintmax_t& max_iter] [An optional maximum number of iterations; on exit, the number of iterations performed. +It is the number of factorizations of /J/, and does not count the line search's evaluations of `f`.]] +] + +The iteration stops when the Newton step /a/ satisfies ‖/a/‖[sub ∞] ≤ 2[super 1−digits] ‖/x/‖[sub ∞] +(the step is then taken), when [*F](/x/) = 0 exactly, or when `max_iter` iterations have been performed. +It also stops when steps below √ε‖/x/‖[sub ∞] no longer shrink at all, since they are then rounding noise; +at a multiple root, where the steps shrink by a constant factor, the iteration goes on. +As for the scalar functions, the tolerance is relative to /x/, so a root at the origin is approached until +[*F] underflows to zero. + +[h4 Globalization] + +Far from a root a full Newton step can increase ‖[*F]‖, so both functions minimize ½‖[*F]‖[super 2] along their steps. +Newton backtracks along /a/ until the Armijo condition ½‖[*F](/x/ + /t a/)‖[super 2] ≤ (1 − 2×10[super −4] /t/)½‖[*F](/x/)‖[super 2] holds, +choosing each new /t/ from the quadratic through ½‖[*F]‖[super 2] and its slope at /t/ = 0, kept in \[/t/\/10, /t/\/2\]. +Halley tries its full step first, and if that does not decrease ‖[*F]‖ it takes the backtracking Newton step instead. +A trial point where [*F] is not finite is treated as an increase, so `f` may return infinities or NaNs outside its domain. + +The componentwise division can blow up where /a/[sub /i/] + /b/[sub /i/]\/2 is near zero. +Cuyt and Rall take 0\/0 as 0 (the Halley component is zero when the Newton component is). +In addition, where /q/ = 1 + /b/[sub /i/]\/(2/a/[sub /i/]) < 1\/2, so that the Halley component points against the Newton +component or is more than twice as long, the Newton component /a/[sub /i/] is used in its place. +Near a root /b/ = /O/(‖/a/‖[super 2]) and /q/ → 1, so this does not affect the cubic convergence. + +[h4 Errors] + +An __evaluation_error is raised through the default policy (and so is thrown by default) when + +* the Jacobian is singular: a pivot of its LU factorization with partial pivoting is not above /n/ε times the largest +entry of its column, +* no step along the Newton direction decreases ‖[*F]‖, as happens near a local minimum of ‖[*F]‖ that is not a root +(where /J/ is singular) or when /J/ is wrong, +* [*F] is not finite at the initial guess, or [*F], /J/ or `d2(v)` has the wrong size. + +A domain error is raised when `x` is empty. + +When no step decreases ‖[*F]‖ but the Newton step is below ε[super 1\/3]‖/x/‖[sub ∞], `x` is returned without an error: +at a local minimum of ‖[*F]‖ that is not a root, /J/ is singular and the Newton step large or undefined, whereas a step +this small is the rounding noise of an accurate /x/, about cond(/J/)ε‖/x/‖, for cond(/J/) up to ε[super −2\/3] +(10[super 10] in double precision). + +With a policy that does not throw, every component of the result is the value returned by the error handler, NaN by default. + +[h4 Supplying the derivatives] + +[*F]″(/x/)\[/v/, /v/\] is the second derivative of the univariate function /t/ ↦ [*F](/x/ + /t v/) at /t/ = 0, +so a single second order forward-mode pass of [link math_toolkit.autodiff autodiff] in the one variable /t/ computes it, for any /n/. +The example [@../../example/multivariate_roots_ellipse.cpp multivariate_roots_ellipse.cpp] finds the point of an ellipse +closest to a point /q/: at the closest point /p/, /p/ is on the ellipse and /p/ − /q/ is along the normal: + +[multivariate_roots_ellipse_system] + +The Jacobian is built one column at a time by first order passes along the coordinate directions, and `d2` is a +second order pass along `v`, computed only when the solver asks for it: + +[multivariate_roots_autodiff] + +[multivariate_roots_ellipse_solve] + +The closure `d2` may capture whatever it needs from the evaluation at /x/; for the Rosenbrock gradient system +∇((1 − /x/)[super 2] + 100(/y/ − /x/[super 2])[super 2]) = 0, written out by hand, + + auto f = [](std::vector const& x) { + double u = x[0], w = x[1]; + std::vector F{-2 * (1 - u) - 400 * u * (w - u * u), 200 * (w - u * u)}; + std::vector> J{{2 - 400 * (w - u * u) + 800 * u * u, -400 * u}, {-400 * u, 200.0}}; + auto d2 = [u](std::vector const& v) { + return std::vector{-800 * v[0] * v[1] + 2400 * u * v[0] * v[0], -400 * v[0] * v[0]}; + }; + return std::make_tuple(F, J, d2); + }; + std::uintmax_t iterations = 1000; + auto x = boost::math::tools::multivariate_halley(f, std::vector{-1.2, 1.0}, 53, iterations); + +[h4 Performance and caveats] + +A Halley iteration costs one more triangular solve and one evaluation of `d2` than a Newton iteration, and saves +about a third of the iterations near the root; which is faster overall depends on the cost of `f` and `d2` +relative to the /O/(/n/[super 3]) factorization. +Far from the root Halley's steps follow curved valleys better: from the classic start (−1.2, 1), Newton with the +line search needs 253 iterations on the Rosenbrock gradient system, as its straight steps creep along the valley, +and Halley needs 68. + +Newton's method is invariant under affine changes of the variables /x/ ↦ /Ax/; the multivariate Halley method is not, +because its update divides componentwise. It is invariant under scaling of the individual variables, and, like Newton, +under any invertible linear change of the equations [*F] ↦ /B/[*F]. + +The LU factorization is dense; for large sparse systems, solve with a sparse factorization in a loop of your own. +The line search guarantees that ‖[*F]‖ decreases, not that the root found is the one nearest the guess. + +[h4 References] + +* A. Cuyt and P. Van der Cruyssen, ['Abstract Padé approximants for the solution of a system of nonlinear equations], +Computers & Mathematics with Applications 9 (1983) 617–624. +* A. Cuyt, ['Numerical stability of the Halley-iteration for the solution of a system of nonlinear equations], +Mathematics of Computation 38 (1982) 171–179. +* A. A. M. Cuyt and L. B. Rall, ['Computational implementation of the multivariate Halley method for solving +nonlinear systems of equations], ACM Transactions on Mathematical Software 11 (1985) 20–36. +* J. E. Dennis and R. B. Schnabel, ['Numerical Methods for Unconstrained Optimization and Nonlinear Equations], +SIAM (1996), section 6.3. + +[endsect] [/section:multivariate_roots Multivariate Newton and Halley] + +[/ + Copyright 2026 Nick Thompson. + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or copy at + http://www.boost.org/LICENSE_1_0.txt). +] diff --git a/doc/roots/roots_overview.qbk b/doc/roots/roots_overview.qbk index c6749516e0..75acb17abf 100644 --- a/doc/roots/roots_overview.qbk +++ b/doc/roots/roots_overview.qbk @@ -5,7 +5,8 @@ Some __root_finding_without_derivatives methods are __bisection, __bracket_solve, including use of __root_finding_TOMS748. For __root_finding_with_derivatives the methods of -__newton, __halley, and __schroder are implemented. +__newton, __halley, and __schroder are implemented, +and for systems of equations, [link math_toolkit.multivariate_roots multivariate Newton and Halley]. For locating minima of a function, a __brent_minima_example is provided. @@ -17,6 +18,7 @@ There are several fully-worked __root_finding_examples, including: [include roots_without_derivatives.qbk] [include roots.qbk] +[include multivariate_roots.qbk] [include cubic_roots.qbk] [include quartic_roots.qbk] [include polynomial_roots.qbk] diff --git a/example/Jamfile.v2 b/example/Jamfile.v2 index ea8e32a9a8..a87ed70348 100644 --- a/example/Jamfile.v2 +++ b/example/Jamfile.v2 @@ -151,6 +151,7 @@ test-suite examples : [ run reciprocal_fibonacci_constant.cpp : : : [ check-target-builds ../config//has_mpfr "MPFR Support" : "-lmpfr -lgmp" : no ] ] [ run reverse_mode_linear_regression_example.cpp ] [ run autodiff_reverse_black_scholes.cpp ] + [ run multivariate_roots_ellipse.cpp : : : [ requires cxx14_generic_lambdas cxx14_return_type_deduction ] ] [ compile leja_parallel_assembly.cpp : [ requires cxx17_if_constexpr cxx11_hdr_thread ] [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] ; diff --git a/example/multivariate_roots_ellipse.cpp b/example/multivariate_roots_ellipse.cpp new file mode 100644 index 0000000000..69490b03ea --- /dev/null +++ b/example/multivariate_roots_ellipse.cpp @@ -0,0 +1,103 @@ +// Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +// The point of the ellipse x^2/a^2 + y^2/b^2 = 1 closest to q, by Newton's method and Cuyt's multivariate Halley +// method, with the Jacobian and the second directional derivative from forward-mode autodiff. + +#include +#include +#include +#include +#include +#include +#include +#include + +//[multivariate_roots_ellipse_system +// At the closest point p, p is on the ellipse and p - q is along the normal (x/a^2, y/b^2). +struct ellipse_closest_point +{ + double a; + double b; + std::array q; + + template + std::array operator()(std::array const& p) const + { + return {(p[0] - q[0]) * p[1] / (b * b) - (p[1] - q[1]) * p[0] / (a * a), + p[0] * p[0] / (a * a) + p[1] * p[1] / (b * b) - 1}; + } +}; +//] + +//[multivariate_roots_autodiff +using boost::math::differentiation::make_fvar; + +// F(x + t v) with its first Order derivatives in t at t = 0: F, J v and F''[v, v]. One forward pass in the single +// variable t, however many variables F has. +template +auto along(System const& F, std::array const& x, std::array const& v) +{ + auto const t = make_fvar(0); + std::array, N> p; + for (std::size_t k = 0; k < N; ++k) + p[k] = x[k] + v[k] * t; + return F(p); +} + +// (F, J, d2) for multivariate_halley: J a column at a time, d2(v) = F''[v, v] when the solver asks for it. +template +auto with_derivatives(System const& F, std::array const& x) +{ + std::array value{}; + std::array, N> J{}; + for (std::size_t j = 0; j < N; ++j) + { + std::array e{}; + e[j] = 1; + auto const Fj = along<1>(F, x, e); + for (std::size_t i = 0; i < N; ++i) + { + value[i] = Fj[i].derivative(0); + J[i][j] = Fj[i].derivative(1); + } + } + auto d2 = [F, x](std::array const& v) { + auto const Fv = along<2>(F, x, v); + std::array w; + for (std::size_t i = 0; i < N; ++i) + w[i] = Fv[i].derivative(2); + return w; + }; + return std::make_tuple(value, J, d2); +} +//] + +int main() +{ + using boost::math::tools::multivariate_halley; + using boost::math::tools::multivariate_newton; + double const a = 3; + double const b = 1; + std::cout << std::setprecision(17); + for (std::array q : {std::array{4, 2}, {1, 0.25}, {-0.5, 3}, {2.9, -0.1}}) + { + //[multivariate_roots_ellipse_solve + ellipse_closest_point F{a, b, q}; + auto f = [&F](std::array const& p) { return with_derivatives(F, p); }; + // Start from the point of the ellipse at the angle of q, stretched to the ellipse's shape: + double const theta = std::atan2(a * q[1], b * q[0]); + std::array const p0{a * std::cos(theta), b * std::sin(theta)}; + std::uintmax_t halley_iterations = 50; + auto p = multivariate_halley(f, p0, std::numeric_limits::digits, halley_iterations); + //] + std::uintmax_t newton_iterations = 50; + auto pn = multivariate_newton(f, p0, std::numeric_limits::digits, newton_iterations); + std::cout << "q = (" << q[0] << ", " << q[1] << "): closest point (" << p[0] << ", " << p[1] + << "), distance " << std::hypot(p[0] - q[0], p[1] - q[1]) << "\n" + << " Halley iterations: " << halley_iterations << ", Newton iterations: " << newton_iterations + << ", |Halley - Newton| = " << std::hypot(p[0] - pn[0], p[1] - pn[1]) << "\n"; + } +} diff --git a/include/boost/math/tools/multivariate_roots.hpp b/include/boost/math/tools/multivariate_roots.hpp new file mode 100644 index 0000000000..8b89c9faf2 --- /dev/null +++ b/include/boost/math/tools/multivariate_roots.hpp @@ -0,0 +1,357 @@ +// (C) Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +// Newton's method and Cuyt's multivariate Halley method for a system F(x) = 0, x in R^n, globalized by a +// backtracking line search on |F|^2/2. +// +// Halley: with the Newton step a = -J^-1 F and b = J^-1 F''[a, a], x_i += a_i^2 / (a_i + b_i/2). +// A. Cuyt and P. Van der Cruyssen, Abstract Pade approximants for the solution of a system of nonlinear +// equations, Comput. Math. Appl. 9 (1983) 617-624. +// A. Cuyt and L. B. Rall, Computational implementation of the multivariate Halley method for solving +// nonlinear systems of equations, ACM TOMS 11 (1985) 20-36. + +#ifndef BOOST_MATH_TOOLS_MULTIVARIATE_ROOTS_HPP +#define BOOST_MATH_TOOLS_MULTIVARIATE_ROOTS_HPP + +#include +#include +#include +#include +#include +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +BOOST_MATH_NAMESPACE_BEGIN +namespace tools { +namespace detail { + +// J[i][j] for a container of rows, J[i*n + j] for a row-major array. +template +struct is_nested_matrix : std::false_type {}; + +template +struct is_nested_matrix()[0][0]))> : std::true_type {}; + +template +bool copy_jacobian(M const& J, std::vector& lu, std::size_t n, std::true_type) +{ + if (J.size() != n) + return false; + for (std::size_t i = 0; i < n; ++i) + { + if (J[i].size() != n) + return false; + for (std::size_t j = 0; j < n; ++j) + lu[i * n + j] = J[i][j]; + } + return true; +} + +template +bool copy_jacobian(M const& J, std::vector& lu, std::size_t n, std::false_type) +{ + if (J.size() != n * n) + return false; + for (std::size_t k = 0; k < n * n; ++k) + lu[k] = J[k]; + return true; +} + +// LU with partial pivoting of the row-major n x n matrix a, in place. False when a pivot is not above n eps times +// the largest entry of its column, so the test does not depend on the scaling of the variables. +template +bool lu_factor(std::vector& a, std::vector& piv, std::vector& colmax, std::size_t n) +{ + using std::abs; + for (std::size_t j = 0; j < n; ++j) + { + colmax[j] = 0; + for (std::size_t i = 0; i < n; ++i) + colmax[j] = (std::max)(colmax[j], Real(abs(a[i * n + j]))); + } + Real const tiny = n * tools::epsilon(); + for (std::size_t k = 0; k < n; ++k) + { + std::size_t p = k; + for (std::size_t i = k + 1; i < n; ++i) + if (abs(a[i * n + k]) > abs(a[p * n + k])) + p = i; + piv[k] = p; + if (p != k) + for (std::size_t j = 0; j < n; ++j) + std::swap(a[k * n + j], a[p * n + j]); + Real const d = a[k * n + k]; + if (!(abs(d) > tiny * colmax[k]) || !(BOOST_MATH_NAMESPACE::isfinite)(d)) + return false; + for (std::size_t i = k + 1; i < n; ++i) + { + Real const l = a[i * n + k] /= d; + for (std::size_t j = k + 1; j < n; ++j) + a[i * n + j] -= l * a[k * n + j]; + } + } + return true; +} + +template +void lu_solve(std::vector const& a, std::vector const& piv, std::size_t n, V& b) +{ + for (std::size_t k = 0; k < n; ++k) + std::swap(b[k], b[piv[k]]); + for (std::size_t i = 1; i < n; ++i) + for (std::size_t j = 0; j < i; ++j) + b[i] -= a[i * n + j] * b[j]; + for (std::size_t i = n; i-- > 0;) + { + for (std::size_t j = i + 1; j < n; ++j) + b[i] -= a[i * n + j] * b[j]; + b[i] /= a[i * n + i]; + } +} + +template +auto inf_norm(V const& v) -> typename std::decay::type +{ + using std::abs; + typename std::decay::type r = 0; + for (std::size_t i = 0; i < v.size(); ++i) + r = (std::max)(r, decltype(r)(abs(v[i]))); + return r; +} + +// The Halley component a^2 / (a + b/2) = a / q, q = 1 + b/(2a). Cuyt and Rall (1985), eq. (7.2), take 0/0 as 0. +// Where q < 1/2 the Pade step points against the Newton step or is more than twice as long, and the iteration is +// far from the cubic regime, where q -> 1: take the Newton component a instead. +template +Real halley_component(Real const& a, Real const& b) +{ + if (a == 0) + return 0; + Real const q = 1 + b / (2 * a); + return q >= Real(0.5) ? Real(a / q) : a; +} + +// b = F''[a, a], from the closure f returned as its third element; nothing for Newton. False if it has the wrong size. +template +bool second_derivative(std::true_type, R& r, V const& a, V& b) +{ + auto const& w = std::get<2>(r)(a); + if (w.size() != a.size()) + return false; + for (std::size_t i = 0; i < a.size(); ++i) + b[i] = w[i]; + return true; +} + +template +bool second_derivative(std::false_type, R&, V const&, V&) +{ + return true; +} + +// x with every component the value of the evaluation error, NaN unless the policy's handler says otherwise. +template +Container fail(Container x, char const* function, char const* message) +{ + using Real = typename std::decay::type; + Real const v = policies::raise_evaluation_error(function, message, std::numeric_limits::quiet_NaN(), + policies::policy<>()); + for (std::size_t i = 0; i < x.size(); ++i) + x[i] = v; + return x; +} + +template +Container multivariate_root_finder(F& f, Container x, int digits, std::uintmax_t& max_iter, char const* function) +{ + using std::sqrt; + using std::ldexp; + using Real = typename std::decay::type; + using R = typename std::decay::type; + using nested = is_nested_matrix(std::declval()))>::type>; + std::size_t const n = x.size(); + std::uintmax_t const limit = max_iter; + max_iter = 0; + if (n == 0) + { + policies::raise_domain_error(function, "x must not be empty", Real(0), policies::policy<>()); + return x; + } + if (digits > tools::digits() || digits < 1) + digits = tools::digits(); + Real const tol = ldexp(Real(1), 1 - digits); + Real const eps = tools::epsilon(); + using std::cbrt; + // A Newton step this small that cannot decrease |F| is rounding noise; see the line search below. + Real const noise = cbrt(eps); + Real const armijo = Real(1) / 10000; + + std::vector lu(n * n), colmax(n); + std::vector piv(n); + Container fx(x), a(x), b(x), xt(x); + // F into fx, J into lu, |F|^2/2 into phi, or infinity if F has a non-finite entry, and whether F = 0 exactly into + // zero (phi can underflow). False if F or J has the wrong shape. + // Every trial point is unpacked, so fx and lu hold the accepted point's values only once it is accepted; within an + // iteration only the steps a and b, already solved for, are used. + bool zero = false; + auto unpack = [&](R const& r, Real& phi) { + auto const& Fv = std::get<0>(r); + if (Fv.size() != n || !copy_jacobian(std::get<1>(r), lu, n, nested())) + return false; + phi = 0; + zero = true; + for (std::size_t i = 0; i < n; ++i) + { + fx[i] = Fv[i]; + phi += fx[i] * fx[i] / 2; + zero = zero && fx[i] == 0; + } + if (!(BOOST_MATH_NAMESPACE::isfinite)(phi)) + phi = std::numeric_limits::infinity(); + return true; + }; + + // The tuple holds a closure for F''[v, v] for Halley, which can be copied but not assigned. + std::unique_ptr r(new R(f(x))); + Real phi; + Real last = std::numeric_limits::infinity(); + if (!unpack(*r, phi)) + return fail(x, function, "F(x) must have x.size() entries and J(x) x.size() rows of x.size()"); + if (phi == std::numeric_limits::infinity()) + return fail(x, function, "F is not finite at the initial guess"); + while (max_iter < limit && !zero) + { + if (!lu_factor(lu, piv, colmax, n)) + return fail(x, function, "The Jacobian is singular"); + for (std::size_t i = 0; i < n; ++i) + a[i] = -fx[i]; + lu_solve(lu, piv, n, a); + ++max_iter; + Real const xnorm = inf_norm(x); + Real const anorm = inf_norm(a); + if (anorm <= tol * xnorm) + { + for (std::size_t i = 0; i < n; ++i) + x[i] += a[i]; + break; + } + // Steps at the rounding level of x that have stopped shrinking are noise: x is as good as it gets. (At a + // multiple root the steps shrink linearly, so the iteration goes on.) + if (anorm <= sqrt(eps) * xnorm && anorm >= last) + break; + last = anorm; + Real phit = 0; + bool newton = true; + BOOST_MATH_IF_CONSTEXPR (Halley) + { + if (!second_derivative(std::integral_constant(), *r, a, b)) + return fail(x, function, "d2(v) must have x.size() entries"); + lu_solve(lu, piv, n, b); + for (std::size_t i = 0; i < n; ++i) + { + Real const c = halley_component(a[i], b[i]); + newton = newton && c == a[i]; + xt[i] = x[i] + c; + } + if (!newton) + { + std::unique_ptr rt(new R(f(xt))); + if (!unpack(*rt, phit)) + return fail(x, function, "F(x) must have x.size() entries and J(x) x.size() rows of x.size()"); + if (phit <= (1 - 2 * armijo) * phi) + { + x = xt; + r = std::move(rt); + phi = phit; + continue; + } + } + } + // Backtrack along the Newton step until |F|^2/2 decreases enough. On the line x + t a, |F|^2/2 has + // slope -2 phi at t = 0; minimize the quadratic through that and the last trial, kept in [t/10, t/2]. + Real t = 1; + for (;;) + { + for (std::size_t i = 0; i < n; ++i) + xt[i] = x[i] + t * a[i]; + std::unique_ptr rt(new R(f(xt))); + if (!unpack(*rt, phit)) + return fail(x, function, "F(x) must have x.size() entries and J(x) x.size() rows of x.size()"); + if (phit <= (1 - 2 * armijo * t) * phi) + { + x = xt; + r = std::move(rt); + phi = phit; + break; + } + if (t * anorm <= tol * xnorm || t < eps) + { + // At a local minimum of |F| that is not a root J is singular, so the Newton step is large or undefined. + // A step below eps^(1/3) |x| that cannot decrease |F| is instead the rounding noise of an accurate x, + // about cond(J) eps |x|, for cond(J) up to eps^(-2/3), 1e10 in double. + if (anorm <= noise * xnorm) + return x; + return fail(x, function, "No step along the Newton direction decreases |F|: x is near a local minimum of |F| that is not a root, or J is wrong"); + } + Real s = t / 10; + if (phit != std::numeric_limits::infinity()) + { + Real const curvature = phit - phi + 2 * phi * t; + if (curvature > 0) + s = phi * t * t / curvature; + } + t = (std::min)((std::max)(s, t / 10), t / 2); + } + } + return x; +} + +} // namespace detail + +// A root of F near x. f(x) returns a pair or tuple (F(x), J(x)), with J a container of rows or a row-major array +// of n*n entries. On exit max_iter is the number of iterations taken. +BOOST_MATH_EXPORT template +Container multivariate_newton(F f, Container x, int digits, std::uintmax_t& max_iter) +{ + return detail::multivariate_root_finder(f, std::move(x), digits, max_iter, + "boost::math::tools::multivariate_newton<%1%>"); +} + +BOOST_MATH_EXPORT template +Container multivariate_newton(F f, Container x, int digits) +{ + std::uintmax_t max_iter = (std::numeric_limits::max)(); + return multivariate_newton(f, std::move(x), digits, max_iter); +} + +// As multivariate_newton, but f(x) returns (F(x), J(x), d2), where d2(v) is F''(x)[v, v], the second derivative +// of t -> F(x + t v) at t = 0. +BOOST_MATH_EXPORT template +Container multivariate_halley(F f, Container x, int digits, std::uintmax_t& max_iter) +{ + return detail::multivariate_root_finder(f, std::move(x), digits, max_iter, + "boost::math::tools::multivariate_halley<%1%>"); +} + +BOOST_MATH_EXPORT template +Container multivariate_halley(F f, Container x, int digits) +{ + std::uintmax_t max_iter = (std::numeric_limits::max)(); + return multivariate_halley(f, std::move(x), digits, max_iter); +} + +} // namespace tools +BOOST_MATH_NAMESPACE_END + +#endif // BOOST_MATH_TOOLS_MULTIVARIATE_ROOTS_HPP diff --git a/module/math.cppm b/module/math.cppm index 2fafcbb2ae..482bcb13f2 100644 --- a/module/math.cppm +++ b/module/math.cppm @@ -152,6 +152,7 @@ extern "C++" { #include #include #include +#include #include #include #include diff --git a/nanobind/src/roots.cpp b/nanobind/src/roots.cpp index c2c8353472..5f5c1eb843 100644 --- a/nanobind/src/roots.cpp +++ b/nanobind/src/roots.cpp @@ -2,14 +2,15 @@ // Distributed under the Boost Software License, Version 1.0. // https://www.boost.org/LICENSE_1_0.txt // -// The iterative root finders of tools/roots.hpp and tools/toms748_solve.hpp on a Python callable, and the -// closed-form quadratic, cubic and quartic solvers. The callable is called with the GIL held, and an -// exception it raises propagates. Boost.Math errors (no sign change, roots out of order, ...) follow +// The iterative root finders of tools/roots.hpp, tools/multivariate_roots.hpp and tools/toms748_solve.hpp on a +// Python callable, and the closed-form quadratic, cubic and quartic solvers. The callable is called with the GIL +// held, and an exception it raises propagates. Boost.Math errors (no sign change, roots out of order, ...) follow // boost.math.errstate as for a scalar. Tolerances are a number of binary digits, as in Boost. #include "bindings.hpp" #include +#include #include #include #include @@ -23,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -59,6 +61,44 @@ bt::eps_tolerance tolerance(int digits) return bt::eps_tolerance(static_cast(digits)); } +// f(x) -> (F, J) or (F, J, d2) for the multivariate solvers: F and d2(v) = F''(x)[v, v] of shape (n,), J of shape +// (n, n), passed on row-major. +struct py_system +{ + nb::callable f; + std::size_t n; + + auto operator()(std::vector const& x) const + { + recording suspended(nullptr); + nb::object r = f(vector_to_array(x)); + std::size_t const n_ = n; + // a as a float64 array of shape (n,), or (n, n) when square, flattened row-major. + auto shaped = [n_](nb::handle a, bool square, char const* what) { + nb::object arr = numpy().attr("asarray")(a, nb::arg("dtype") = numpy().attr("float64")); + nb::object shape = arr.attr("shape"); + bool ok = nb::len(shape) == (square ? 2u : 1u); + for (std::size_t i = 0; ok && i < nb::len(shape); ++i) + ok = nb::cast(shape[nb::int_(i)]) == n_; + if (!ok) + { + std::string got = nb::str(shape).c_str(); + throw nb::value_error((std::string(what) + " has shape " + got + ", expected (" + std::to_string(n_) + + (square ? ", " + std::to_string(n_) : std::string(",")) + ")").c_str()); + } + return array_to_vector(arr); + }; + nb::object d2 = nb::len(r) > 2 ? nb::object(r[2]) : nb::none(); + auto second = [d2, shaped](std::vector const& v) { + recording inner(nullptr); + if (d2.is_none()) + throw nb::value_error("multivariate_halley: f(x) must return (F, J, d2)"); + return shaped(d2(vector_to_array(v)), false, "d2(v)"); + }; + return std::make_tuple(shaped(r[0], false, "F(x)"), shaped(r[1], true, "J(x)"), second); + } +}; + // The real roots, ascending; Boost pads its fixed-size result with NaN. template auto real_roots(std::array const& roots) @@ -139,6 +179,26 @@ void bind_roots(nb::module_& m) "A root of the analytic f near guess by Newton's method; f(z) returns (f(z), f'(z)); nan+nanj if it fails." BMPY_URL("roots_deriv")); + auto def_multivariate = [&](char const* name, auto iterate, char const* sig, char const* doc) { + m.def(name, [newton_result, name, iterate](nb::callable f, nb::handle x0, int digits, std::uintmax_t max_iterations) { + std::vector x = array_to_vector(x0); + if (x.empty()) + throw nb::value_error((std::string(name) + ": x0 must not be empty").c_str()); + py_system system{f, x.size()}; + auto root = solve(name, [&] { return iterate(system, x, digits, max_iterations); }); + return newton_result(vector_to_array(std::move(root)), max_iterations); + }, arg(), arg(), nb::kw_only(), arg("digits") = default_digits, arg("max_iterations") = default_max_iterations, + nb::sig(sig), doc); + }; + def_multivariate("multivariate_newton", [](auto&&... a) { return bt::multivariate_newton(a...); }, + "def multivariate_newton(f, x0, /, *, digits=53, max_iterations=1000) -> NewtonResult", + "A root of the system F(x) = 0 near x0 by Newton's method with a line search; f(x) returns (F(x), J(x)), " + "arrays of shape (n,) and (n, n), J[i, j] = dF_i/dx_j." BMPY_URL("multivariate_roots")); + def_multivariate("multivariate_halley", [](auto&&... a) { return bt::multivariate_halley(a...); }, + "def multivariate_halley(f, x0, /, *, digits=53, max_iterations=1000) -> NewtonResult", + "A root of the system F(x) = 0 near x0 by Cuyt's multivariate Halley method; f(x) returns (F(x), J(x), d2), " + "where d2(v) is F''(x)[v, v], the second derivative of F(x + t v) in t at 0." BMPY_URL("multivariate_roots")); + m.def("quadratic_roots", [](double a, double b, double c) { auto [x0, x1] = bt::quadratic_roots(a, b, c); return real_roots(std::array{x0, x1}); diff --git a/nanobind/tests/test_roots.py b/nanobind/tests/test_roots.py index fb128bcda1..25607631a9 100644 --- a/nanobind/tests/test_roots.py +++ b/nanobind/tests/test_roots.py @@ -56,3 +56,31 @@ def test_polynomials(): # x^2 - 2, coefficients from the constant term up. for solve in (roots.polynomial_roots, roots.polynomial_roots_aberth, roots.polynomial_roots_companion_matrix): assert np.allclose(sorted(solve([-2.0, 0.0, 1.0]).real), [-SQRT2, SQRT2]) + + +def test_multivariate(): + # x^2 + y^2 = 4 and x = 2y; the root from (1, 1) is (4, 2)/sqrt(5). + def f(x): + F = np.array([x[0] ** 2 + x[1] ** 2 - 4, x[0] - 2 * x[1]]) + J = np.array([[2 * x[0], 2 * x[1]], [1.0, -2.0]]) + return F, J, lambda v: np.array([2 * (v[0] ** 2 + v[1] ** 2), 0.0]) + + expected = np.array([4, 2]) / math.sqrt(5) + for solve in (roots.multivariate_newton, roots.multivariate_halley): + r = solve(f, [1.0, 1.0]) + assert type(r) is roots.NewtonResult and np.allclose(r.root, expected, rtol=1e-15, atol=0) + assert 0 < r.iterations < 10 + # Newton ignores d2, and accepts a pair: + assert np.allclose(roots.multivariate_newton(lambda x: f(x)[:2], [1.0, 1.0]).root, expected, rtol=1e-15) + assert roots.multivariate_halley(f, [3.0, 1.0]).iterations <= roots.multivariate_newton(f, [3.0, 1.0]).iterations + # J = 0 at the origin: + with pytest.raises(RuntimeError, match="singular"): + roots.multivariate_newton(lambda x: (x * x - 1, np.diag(2 * x)), [0.0, 0.0]) + with pytest.raises(ValueError, match="J"): + roots.multivariate_newton(lambda x: (x * x - 1, np.ones(3)), [1.0, 2.0]) + with pytest.raises(ValueError, match="J"): + roots.multivariate_newton(lambda x: (x * x - 1, np.ones(4)), [1.0, 2.0]) + with pytest.raises(ValueError, match="F"): + roots.multivariate_newton(lambda x: (np.diag(x * x - 1), np.diag(2 * x)), [1.0, 2.0]) + with pytest.raises(ValueError, match="d2"): + roots.multivariate_halley(lambda x: f(x)[:2], [1.0, 1.0]) diff --git a/reporting/performance/multivariate_roots_performance.cpp b/reporting/performance/multivariate_roots_performance.cpp new file mode 100644 index 0000000000..c6d606db8a --- /dev/null +++ b/reporting/performance/multivariate_roots_performance.cpp @@ -0,0 +1,113 @@ +// Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +/* + * Google Benchmark suite: time to solution of multivariate_newton and multivariate_halley, with the iterations + * and evaluations of f as counters. + * + * Example: + * c++ -O2 -DNDEBUG -std=c++14 -Iinclude \ + * reporting/performance/multivariate_roots_performance.cpp \ + * -lbenchmark -lpthread -o multivariate_roots_bench + */ + +#include +#include + +#include +#include +#include + +namespace { + +using boost::math::tools::multivariate_halley; +using boost::math::tools::multivariate_newton; +using vector = std::vector; + +// Broyden's tridiagonal function, J row-major. +struct broyden_tridiagonal +{ + std::size_t* calls; + + auto operator()(vector const& x) const + { + ++*calls; + std::size_t const n = x.size(); + vector F(n); + vector J(n * n, 0.0); + for (std::size_t i = 0; i < n; ++i) + { + double const left = i > 0 ? x[i - 1] : 0.0; + double const right = i + 1 < n ? x[i + 1] : 0.0; + F[i] = (3 - 2 * x[i]) * x[i] - left - 2 * right + 1; + J[i * n + i] = 3 - 4 * x[i]; + if (i > 0) + J[i * n + i - 1] = -1; + if (i + 1 < n) + J[i * n + i + 1] = -2; + } + auto d2 = [](vector const& v) { + vector w(v.size()); + for (std::size_t i = 0; i < v.size(); ++i) + w[i] = -4 * v[i] * v[i]; + return w; + }; + return std::make_tuple(F, J, d2); + } +}; + +// The gradient of the Rosenbrock function, from the classic start (-1.2, 1). +struct rosenbrock_gradient +{ + std::size_t* calls; + + auto operator()(vector const& x) const + { + ++*calls; + double const u = x[0]; + double const w = x[1]; + vector F{-2 * (1 - u) - 400 * u * (w - u * u), 200 * (w - u * u)}; + vector J{2 - 400 * (w - u * u) + 800 * u * u, -400 * u, -400 * u, 200.0}; + auto d2 = [u](vector const& v) { return vector{-800 * v[0] * v[1] + 2400 * u * v[0] * v[0], -400 * v[0] * v[0]}; }; + return std::make_tuple(F, J, d2); + } +}; + +template +void run(benchmark::State& state, vector const& x0) +{ + std::size_t calls = 0; + std::uintmax_t iterations = 0; + for (auto _ : state) + { + calls = 0; + iterations = 10000; + System f{&calls}; + auto x = Halley ? multivariate_halley(f, x0, 53, iterations) : multivariate_newton(f, x0, 53, iterations); + benchmark::DoNotOptimize(x.data()); + } + state.counters["iterations"] = static_cast(iterations); + state.counters["f_calls"] = static_cast(calls); +} + +template +void BM_broyden(benchmark::State& state) +{ + run(state, vector(static_cast(state.range(0)), -1.0)); +} +BENCHMARK_TEMPLATE(BM_broyden, false)->Arg(10)->Arg(100)->Unit(benchmark::kMicrosecond); +BENCHMARK_TEMPLATE(BM_broyden, true)->Arg(10)->Arg(100)->Unit(benchmark::kMicrosecond); + +template +void BM_rosenbrock_gradient(benchmark::State& state) +{ + run(state, vector{-1.2, 1.0}); +} +BENCHMARK_TEMPLATE(BM_rosenbrock_gradient, false)->Unit(benchmark::kMicrosecond); +BENCHMARK_TEMPLATE(BM_rosenbrock_gradient, true)->Unit(benchmark::kMicrosecond); + +} // namespace + +BENCHMARK_MAIN(); diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index c6fedef6c5..42076f8665 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -1103,6 +1103,7 @@ test-suite misc : [ run test_root_iterations.cpp pch /boost/test//boost_unit_test_framework : : : [ requires cxx11_hdr_tuple ] ] [ run test_root_finding_concepts.cpp /boost/test//boost_unit_test_framework ] [ run test_toms748_solve.cpp pch /boost/test//boost_unit_test_framework ] + [ run test_multivariate_roots.cpp : : : [ requires cxx14_return_type_deduction cxx14_generic_lambdas ] ] [ run compile_test/test_interpolators_cubic_spline_incl.cpp compile_test_main : : : [ requires cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations ] [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/test_interpolators_barycentric_rational_incl.cpp compile_test_main : : : [ requires cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations cxx11_unified_initialization_syntax ] [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run test_octonion.cpp @@ -1830,6 +1831,7 @@ test-suite concepts : [ compile compile_test/test_tools_real_cast_inc.cpp : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ compile compile_test/test_tools_remez_inc.cpp : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ compile compile_test/test_tools_roots_inc.cpp : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] + [ compile compile_test/test_tools_multivariate_roots_incl.cpp : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ compile compile_test/test_tools_series_inc.cpp : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ compile compile_test/test_tools_solve_inc.cpp : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ compile compile_test/test_tools_stats_inc.cpp : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] diff --git a/test/compile_test/test_tools_multivariate_roots_incl.cpp b/test/compile_test/test_tools_multivariate_roots_incl.cpp new file mode 100644 index 0000000000..5685282db2 --- /dev/null +++ b/test/compile_test/test_tools_multivariate_roots_incl.cpp @@ -0,0 +1,23 @@ +// Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +// Basic sanity check that header +// #includes all the files that it needs to. +// +#include +// +// Note this header includes no other headers, this is +// important if this test is to be meaningful: +// +#include "test_compile_result.hpp" + +void compile_and_link_test() +{ + using V = std::vector; + auto f = [](V const& x) { return std::make_pair(V{x[0] - 1}, V{1.0}); }; + std::uintmax_t max_iter = 10; + check_result(boost::math::tools::multivariate_newton(f, V{0.0}, 53, max_iter)); + check_result(boost::math::tools::multivariate_newton(f, V{0.0}, 53)); +} diff --git a/test/test_multivariate_roots.cpp b/test/test_multivariate_roots.cpp new file mode 100644 index 0000000000..5152c9d8b4 --- /dev/null +++ b/test/test_multivariate_roots.cpp @@ -0,0 +1,390 @@ +// Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include "math_unit_test.hpp" +#include +#include +#include +#include +#include +#include +#include +#include + +using boost::math::tools::multivariate_halley; +using boost::math::tools::multivariate_newton; + +// Each system returns (F, J, F''[v, v]); Newton ignores the third element. + +// x^2 + y^2 = 4 and x = 2y, with roots +-(4, 2)/sqrt(5). J is a std::array of rows. +template +auto circle_line(std::array const& x) +{ + std::array F{x[0] * x[0] + x[1] * x[1] - 4, x[0] - 2 * x[1]}; + std::array, 2> J{{{2 * x[0], 2 * x[1]}, {1, -2}}}; + auto d2 = [](std::array const& v) { return std::array{2 * (v[0] * v[0] + v[1] * v[1]), 0}; }; + return std::make_tuple(F, J, d2); +} + +// The gradient of the Rosenbrock function (1 - x)^2 + 100(y - x^2)^2, zero only at (1, 1). +template +auto rosenbrock_gradient(std::vector const& x) +{ + Real const u = x[0]; + Real const w = x[1]; + std::vector F{-2 * (1 - u) - 400 * u * (w - u * u), 200 * (w - u * u)}; + std::vector> J{{2 - 400 * (w - u * u) + 800 * u * u, -400 * u}, {-400 * u, Real(200)}}; + auto d2 = [u](std::vector const& v) { + return std::vector{-800 * v[0] * v[1] + 2400 * u * v[0] * v[0], -400 * v[0] * v[0]}; + }; + return std::make_tuple(F, J, d2); +} + +// Broyden's tridiagonal function (More, Garbow and Hillstrom 1981, #30); J is row-major in a std::vector. +template +auto broyden_tridiagonal(std::vector const& x) +{ + std::size_t const n = x.size(); + std::vector F(n); + std::vector J(n * n, Real(0)); + for (std::size_t i = 0; i < n; ++i) + { + Real const left = i > 0 ? x[i - 1] : Real(0); + Real const right = i + 1 < n ? x[i + 1] : Real(0); + F[i] = (3 - 2 * x[i]) * x[i] - left - 2 * right + 1; + J[i * n + i] = 3 - 4 * x[i]; + if (i > 0) + J[i * n + i - 1] = -1; + if (i + 1 < n) + J[i * n + i + 1] = -2; + } + auto d2 = [](std::vector const& v) { + std::vector w(v.size()); + for (std::size_t i = 0; i < v.size(); ++i) + w[i] = -4 * v[i] * v[i]; + return w; + }; + return std::make_tuple(F, J, d2); +} + +// Powell's badly scaled function (More, Garbow and Hillstrom 1981, #3), root near (1.098e-5, 9.106). +template +auto powell_badly_scaled(std::vector const& x) +{ + using std::exp; + Real const ex = exp(-x[0]); + Real const ey = exp(-x[1]); + std::vector F{10000 * x[0] * x[1] - 1, ex + ey - Real(1.0001)}; + std::vector> J{{10000 * x[1], 10000 * x[0]}, {-ex, -ey}}; + auto d2 = [ex, ey](std::vector const& v) { + return std::vector{20000 * v[0] * v[1], ex * v[0] * v[0] + ey * v[1] * v[1]}; + }; + return std::make_tuple(F, J, d2); +} + +// x^2 = 1, y^2 = 1: J is singular wherever x = 0 or y = 0. +template +auto squares(std::vector const& x) +{ + std::vector F{x[0] * x[0] - 1, x[1] * x[1] - 1}; + std::vector J{2 * x[0], 0, 0, 2 * x[1]}; + auto d2 = [](std::vector const& v) { return std::vector{2 * v[0] * v[0], 2 * v[1] * v[1]}; }; + return std::make_tuple(F, J, d2); +} + +template +Real residual(std::vector const& F) +{ + using std::abs; + Real r = 0; + for (Real f : F) + r = (std::max)(r, Real(abs(f))); + return r; +} + +template +void test_circle_line() +{ + int const digits = std::numeric_limits::digits; + Real const s = 1 / std::sqrt(Real(5)); + for (bool halley : {false, true}) + { + std::uintmax_t iters = 100; + std::array x0{1, 1}; + auto f = [](std::array const& x) { return circle_line(x); }; + auto x = halley ? multivariate_halley(f, x0, digits, iters) : multivariate_newton(f, x0, digits, iters); + CHECK_ULP_CLOSE(4 * s, x[0], 2); + CHECK_ULP_CLOSE(2 * s, x[1], 2); + CHECK_LE(iters, std::uintmax_t(8)); + // From the other side, the other root: + x = halley ? multivariate_halley(f, std::array{-3, 1}, digits) : multivariate_newton(f, std::array{-3, 1}, digits); + CHECK_ULP_CLOSE(-4 * s, x[0], 2); + CHECK_ULP_CLOSE(-2 * s, x[1], 2); + } +} + +template +void test_rosenbrock() +{ + int const digits = std::numeric_limits::digits; + auto f = [](std::vector const& x) { return rosenbrock_gradient(x); }; + std::uintmax_t newton_iters = 1000; + std::uintmax_t halley_iters = 1000; + auto xn = multivariate_newton(f, std::vector{-1.2, 1}, digits, newton_iters); + auto xh = multivariate_halley(f, std::vector{-1.2, 1}, digits, halley_iters); + for (auto const& x : {xn, xh}) + { + CHECK_ULP_CLOSE(Real(1), x[0], 4); + CHECK_ULP_CLOSE(Real(1), x[1], 8); + } + // The line search makes Newton creep along the curved valley; Halley's steps bend with it. + CHECK_LE(newton_iters, std::uintmax_t(300)); + CHECK_LE(halley_iters, std::uintmax_t(100)); +} + +template +void test_broyden_tridiagonal() +{ + int const digits = std::numeric_limits::digits; + auto f = [](std::vector const& x) { return broyden_tridiagonal(x); }; + std::uintmax_t newton_iters = 100; + std::uintmax_t halley_iters = 100; + auto xn = multivariate_newton(f, std::vector(10, -1), digits, newton_iters); + auto xh = multivariate_halley(f, std::vector(10, -1), digits, halley_iters); + Real const tol = 8 * std::numeric_limits::epsilon(); + CHECK_LE(residual(std::get<0>(broyden_tridiagonal(xn))), tol); + CHECK_LE(residual(std::get<0>(broyden_tridiagonal(xh))), tol); + for (std::size_t i = 0; i < 10; ++i) + CHECK_ULP_CLOSE(xn[i], xh[i], 8); + CHECK_LE(newton_iters, std::uintmax_t(8)); + // In float the last iterations are at the rounding level, and their number is noise. + CHECK_LE(halley_iters, newton_iters + (std::is_same::value ? 2 : 0)); +} + +template +void test_powell_badly_scaled() +{ + int const digits = std::numeric_limits::digits; + auto f = [](std::vector const& x) { return powell_badly_scaled(x); }; + for (bool halley : {false, true}) + { + std::uintmax_t iters = 200; + std::vector x0{0, 1}; + auto x = halley ? multivariate_halley(f, x0, digits, iters) : multivariate_newton(f, x0, digits, iters); + CHECK_LE(iters, std::uintmax_t(200 - 1)); + // 1.0001 is not representable, so compare with the root of the system as rounded: + CHECK_LE(residual(std::get<0>(powell_badly_scaled(x))), 8 * std::numeric_limits::epsilon()); + if (std::is_same::value) + { + CHECK_RELATIVE_ERROR(Real(1.0981593296997149e-05L), x[0], Real(1e-12)); + CHECK_RELATIVE_ERROR(Real(9.1061467398665240L), x[1], Real(1e-12)); + } + } +} + +template +void test_singular() +{ + int const digits = std::numeric_limits::digits; + auto f = [](std::vector const& x) { return squares(x); }; + // J = 0 at the start: + CHECK_THROW(multivariate_newton(f, std::vector{0, 0}, digits), boost::math::evaluation_error); + CHECK_THROW(multivariate_halley(f, std::vector{0, 0.5}, digits), boost::math::evaluation_error); + // Nearly singular: the Newton step is huge, and the line search pulls it back. + for (bool halley : {false, true}) + { + std::vector x0{Real(1) / 1024, Real(0.5)}; + auto x = halley ? multivariate_halley(f, x0, digits) : multivariate_newton(f, x0, digits); + CHECK_ULP_CLOSE(Real(1), x[0], 1); + CHECK_ULP_CLOSE(Real(1), x[1], 1); + } + // x^2 + 1 = 0 has no real root; |F| is least at 0, where J is singular. + auto g = [](std::array const& x) { + return std::make_pair(std::array{x[0] * x[0] + 1}, std::array{2 * x[0]}); + }; + CHECK_THROW(multivariate_newton(g, std::array{1}, digits), boost::math::evaluation_error); +} + +// The max-norm errors of the iterates from (2.5, 0.5) on circle_line, the 0th to the 7th. +template +std::vector circle_line_errors(bool halley) +{ + using std::abs; + using std::sqrt; + auto f = [](std::array const& x) { return circle_line(x); }; + Real const s = 1 / sqrt(Real(5)); + std::vector e; + for (std::uintmax_t k = 0; k < 8; ++k) + { + std::uintmax_t iters = k; + std::array x0{Real(2.5), Real(0.5)}; + auto x = halley ? multivariate_halley(f, x0, std::numeric_limits::digits, iters) + : multivariate_newton(f, x0, std::numeric_limits::digits, iters); + e.push_back((std::max)(abs(x[0] - 4 * s), abs(x[1] - 2 * s))); + } + return e; +} + +// Near the root e_{k+1} <= e_k^2 for Newton and e_k^3 for Halley (the constants are about 0.3 and 0.5), as long +// as the error is above rounding. +template +void test_convergence_order() +{ + for (bool halley : {false, true}) + { + auto e = circle_line_errors(halley); + int const p = halley ? 3 : 2; + int checked = 0; + for (std::size_t k = 0; k + 1 < e.size(); ++k) + { + if (e[k] < Real(0.1) && e[k + 1] > 10 * std::numeric_limits::epsilon()) + { + CHECK_LE(e[k + 1], Real(std::pow(e[k], p))); + ++checked; + } + } + CHECK_GE(checked, 1); + } +} + +// The order itself, log(e_{k+2}/e_{k+1}) / log(e_{k+1}/e_k), in 50 digits, from the last three errors well above +// rounding. +void test_convergence_order_multiprecision() +{ + using Real = boost::multiprecision::cpp_bin_float_50; + for (bool halley : {false, true}) + { + auto e = circle_line_errors(halley); + Real const floor = 1e-40; + std::size_t k = 0; + while (k + 3 < e.size() && e[k + 3] > floor) + ++k; + Real const order = log(e[k + 2] / e[k + 1]) / log(e[k + 1] / e[k]); + CHECK_LE(abs(order - (halley ? 3 : 2)), Real(0.1)); + } +} + +// At a double or triple root Newton and Halley converge only linearly, with steps shrinking by a constant factor; they +// must not be mistaken for rounding noise. Here F = (x - 1)^m is computed exactly near the root, so the iteration +// can go all the way (with rounding errors in F, only about digits/m digits would be attainable). +template +void test_multiple_roots() +{ + using std::abs; + int const digits = boost::math::tools::digits(); + for (int m : {2, 3}) + { + auto f = [m](std::vector const& x) { + Real const d = x[0] - 1; + Real const d2m = m == 2 ? Real(1) : d; // d^(m-2) + std::vector F{d * d * d2m, x[1] - 2}; + std::vector J{m * d * d2m, 0, 0, 1}; + auto d2 = [m, d2m](std::vector const& v) { return std::vector{m * (m - 1) * d2m * v[0] * v[0], 0}; }; + return std::make_tuple(F, J, d2); + }; + for (bool halley : {false, true}) + { + for (Real start : {Real(1.001), Real(1.5), Real(3)}) + { + std::uintmax_t iters = 1000; + std::vector x0{start, 0}; + auto x = halley ? multivariate_halley(f, x0, digits, iters) : multivariate_newton(f, x0, digits, iters); + CHECK_LE(abs(x[0] - 1), 8 * boost::math::tools::epsilon()); + CHECK_ULP_CLOSE(Real(2), x[1], 1); + CHECK_LE(iters, std::uintmax_t(1000 - 1)); + } + } + } +} + +// cond(J) ~ 4e6, and F carries an erratic error of 1e-13 |x|, as a long computation of F might: near the root (1000, +// 1000) the Newton step is that error amplified by cond(J), about 1e-7 |x|. It is too big to be ignored and cannot +// decrease |F|, but x is as accurate as F allows, so that is convergence, not a failure. +void test_noisy_ill_conditioned() +{ + using std::abs; + using std::sin; + double const delta = 1e-6; + double const c = 1000; + auto f = [delta, c](std::vector const& x) { + double const u = x[0] - c; + double const w = x[1] - c; + double const p = u + w; + std::vector F{p + p * p / (10 * c), u + (1 + delta) * w}; + for (std::size_t i = 0; i < 2; ++i) + F[i] += 1e-13 * c * sin(1e6 * (x[0] + (2 + i) * x[1])); + std::vector> J{{1 + p / (5 * c), 1 + p / (5 * c)}, {1, 1 + delta}}; + auto d2 = [c](std::vector const& v) { + double const s = v[0] + v[1]; + return std::vector{s * s / (5 * c), 0}; + }; + return std::make_tuple(F, J, d2); + }; + for (bool halley : {false, true}) + { + for (auto const& x0 : std::vector>{{500, 1750}, {900, 1550}, {1250, 1375}, {2000, 1000}, {1100, 1000}}) + { + std::vector x{0, 0}; + CHECK_NO_THROW(x = halley ? multivariate_halley(f, x0, 53) : multivariate_newton(f, x0, 53)); + CHECK_LE(abs(x[0] / c - 1), 1e-6); + CHECK_LE(abs(x[1] / c - 1), 1e-6); + } + } +} + +// |F|^2 underflows long before F does; only F = 0 exactly stops the iteration. +void test_tiny_residual() +{ + auto f = [](std::array const& x) { + auto r = circle_line(x); + auto F = std::get<0>(r); + auto J = std::get<1>(r); + for (std::size_t i = 0; i < 2; ++i) + { + F[i] *= 1e-160; + for (std::size_t j = 0; j < 2; ++j) + J[i][j] *= 1e-160; + } + return std::make_pair(F, J); + }; + auto x = multivariate_newton(f, std::array{1, 1}, 53); + CHECK_ULP_CLOSE(4 / std::sqrt(5.0), x[0], 2); + CHECK_ULP_CLOSE(2 / std::sqrt(5.0), x[1], 2); +} + +void test_bad_arguments() +{ + auto f = [](std::vector const& x) { + auto d2 = [](std::vector const&) { return std::vector{0.0}; }; + return std::make_tuple(std::vector{x[0] * x[0] - 2, x[1] - 1}, std::vector{2 * x[0], 0, 0, 1}, d2); + }; + CHECK_THROW(multivariate_halley(f, std::vector{3, 0}, 53), boost::math::evaluation_error); + auto g = [](std::vector const& x) { return std::make_pair(x, x); }; + CHECK_THROW(multivariate_newton(g, std::vector{}, 53), std::domain_error); + CHECK_THROW(multivariate_newton(g, std::vector{1, 2}, 53), boost::math::evaluation_error); +} + +int main() +{ + test_circle_line(); + test_circle_line(); + test_rosenbrock(); + test_rosenbrock(); + test_broyden_tridiagonal(); + test_broyden_tridiagonal(); + test_powell_badly_scaled(); + test_powell_badly_scaled(); + test_singular(); + test_singular(); + test_convergence_order(); + test_convergence_order(); + test_convergence_order_multiprecision(); + test_multiple_roots(); + test_multiple_roots(); + test_noisy_ill_conditioned(); + test_tiny_residual(); + test_bad_arguments(); + return boost::math::test::report_errors(); +}