|
| 1 | +<!DOCTYPE html> |
| 2 | +<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>Efficient Halley's method for nonlinear solving · TaylorDiff.jl</title><meta name="title" content="Efficient Halley's method for nonlinear solving · TaylorDiff.jl"/><meta property="og:title" content="Efficient Halley's method for nonlinear solving · TaylorDiff.jl"/><meta property="twitter:title" content="Efficient Halley's method for nonlinear solving · TaylorDiff.jl"/><meta name="description" content="Documentation for TaylorDiff.jl."/><meta property="og:description" content="Documentation for TaylorDiff.jl."/><meta property="twitter:description" content="Documentation for TaylorDiff.jl."/><meta property="og:url" content="https://juliadiff.org/TaylorDiff.jl/examples/halley/"/><meta property="twitter:url" content="https://juliadiff.org/TaylorDiff.jl/examples/halley/"/><link rel="canonical" href="https://juliadiff.org/TaylorDiff.jl/examples/halley/"/><script data-outdated-warner src="../../assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.050/juliamono.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.16.8/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="../.."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="../../assets/documenter.js"></script><script src="../../search_index.js"></script><script src="../../siteinfo.js"></script><script src="../../../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/catppuccin-mocha.css" data-theme-name="catppuccin-mocha"/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/catppuccin-macchiato.css" data-theme-name="catppuccin-macchiato"/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/catppuccin-frappe.css" data-theme-name="catppuccin-frappe"/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/catppuccin-latte.css" data-theme-name="catppuccin-latte"/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="../../assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href="../../">TaylorDiff.jl</a></span></div><button class="docs-search-query input is-rounded is-small is-clickable my-2 mx-auto py-1 px-2" id="documenter-search-query">Search docs (Ctrl + /)</button><ul class="docs-menu"><li><a class="tocitem" href="../../">Home</a></li><li><span class="tocitem">Examples</span><ul><li class="is-active"><a class="tocitem" href>Efficient Halley's method for nonlinear solving</a><ul class="internal"><li><a class="tocitem" href="#Introduction"><span>Introduction</span></a></li><li><a class="tocitem" href="#Implementation"><span>Implementation</span></a></li></ul></li></ul></li><li><a class="tocitem" href="../../theory/">Theory</a></li><li><a class="tocitem" href="../../api/">API</a></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><a class="docs-sidebar-button docs-navbar-link fa-solid fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a><nav class="breadcrumb"><ul class="is-hidden-mobile"><li><a class="is-disabled">Examples</a></li><li class="is-active"><a href>Efficient Halley's method for nonlinear solving</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href>Efficient Halley's method for nonlinear solving</a></li></ul></nav><div class="docs-right"><a class="docs-navbar-link" href="https://github.com/JuliaDiff/TaylorDiff.jl/blob/main/docs/src/examples/halley.md#" title="Edit source on GitHub"><span class="docs-icon fa-solid"></span></a><a class="docs-settings-button docs-navbar-link fa-solid fa-gear" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-article-toggle-button fa-solid fa-chevron-up" id="documenter-article-toggle-button" href="javascript:;" title="Collapse all docstrings"></a></div></header><article class="content" id="documenter-page"><h1 id="Efficient-Halley's-method-for-nonlinear-solving"><a class="docs-heading-anchor" href="#Efficient-Halley's-method-for-nonlinear-solving">Efficient Halley's method for nonlinear solving</a><a id="Efficient-Halley's-method-for-nonlinear-solving-1"></a><a class="docs-heading-anchor-permalink" href="#Efficient-Halley's-method-for-nonlinear-solving" title="Permalink"></a></h1><h2 id="Introduction"><a class="docs-heading-anchor" href="#Introduction">Introduction</a><a id="Introduction-1"></a><a class="docs-heading-anchor-permalink" href="#Introduction" title="Permalink"></a></h2><p>Say we have a system of <span>$n$</span> equations with <span>$n$</span> unknowns <span>$f(x)=0$</span>, and <span>$f\in \mathbb R^n\to\mathbb R^n$</span> is sufficiently smooth.</p><p>Given a initial guess <span>$x_0$</span>, Newton's method finds a solution by iterating like</p><p class="math-container">\[x_{i+1}=x_i-J(x_i)^{-1}f(x_i)\]</p><p>and this method converges quadratically.</p><p>We can make it converge faster using higher-order derivative information. For example, Halley's method iterates like</p><p class="math-container">\[x_{i+1}=x_i-(a_i\odot a_i)\oslash(a_i-b_i/2)\]</p><p>where the vector multiplication and division <span>$\odot,\oslash$</span> are defined element-wise, and term <span>$a_i$</span> and <span>$b_i$</span> are defined by <span>$J(x_i)a_i = f(x_i)$</span> and <span>$J(x_i)b_i = H(x_i)a_ia_i$</span>.</p><p>Halley's method is proved to converge cubically, which is faster than Newton's method. Here, we demonstrate that with TaylorDiff.jl, you can compute the Hessian-vector-vector product <span>$H(x_i)a_ia_i$</span> very efficiently, such that the Halley's method is almost as cheap as Newton's method per iteration.</p><h2 id="Implementation"><a class="docs-heading-anchor" href="#Implementation">Implementation</a><a id="Implementation-1"></a><a class="docs-heading-anchor-permalink" href="#Implementation" title="Permalink"></a></h2><p>We first define the two iteration schemes mentioned above:</p><pre><code class="language-julia hljs">using TaylorDiff, LinearAlgebra |
| 3 | +import ForwardDiff |
| 4 | + |
| 5 | +function newton(f, x, p; tol = 1e-12, maxiter = 100) |
| 6 | + fp = Base.Fix2(f, p) |
| 7 | + for i in 1:maxiter |
| 8 | + fx = fp(x) |
| 9 | + error = norm(fx) |
| 10 | + println("Iteration $i: x = $x, f(x) = $fx, error = $error") |
| 11 | + error < tol && return |
| 12 | + J = ForwardDiff.jacobian(fp, x) |
| 13 | + a = J \ fx |
| 14 | + @. x -= a |
| 15 | + end |
| 16 | +end |
| 17 | + |
| 18 | +function halley(f, x, p; tol = 1e-12, maxiter = 100) |
| 19 | + fp = Base.Fix2(f, p) |
| 20 | + for i in 1:maxiter |
| 21 | + fx = f(x, p) |
| 22 | + error = norm(fx) |
| 23 | + println("Iteration $i: x = $x, f(x) = $fx, error = $error") |
| 24 | + error < tol && return |
| 25 | + J = ForwardDiff.jacobian(fp, x) |
| 26 | + a = J \ fx |
| 27 | + hvvp = derivative(fp, x, a, Val(2)) |
| 28 | + b = J \ hvvp |
| 29 | + @. x -= (a * a) / (a - b / 2) |
| 30 | + end |
| 31 | +end</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">halley (generic function with 1 method)</code></pre><p>Note that in Halley's method, the hessian-vector-vector product is computed with <code>derivative(fp, x, a, Val(2))</code>. It is guaranteed that asymptotically this is only taking 2x more time compared to evaluating <code>fp(x)</code> itself.</p><p>Now we define some test function:</p><pre><code class="language-julia hljs">f(x, p) = x .* x - p</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">f (generic function with 1 method)</code></pre><p>The Newton's method takes 6 iterations to converge:</p><pre><code class="language-julia hljs">newton(f, [1., 1.], [2., 2.])</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">Iteration 1: x = [1.0, 1.0], f(x) = [-1.0, -1.0], error = 1.4142135623730951 |
| 32 | +Iteration 2: x = [1.5, 1.5], f(x) = [0.25, 0.25], error = 0.3535533905932738 |
| 33 | +Iteration 3: x = [1.4166666666666667, 1.4166666666666667], f(x) = [0.006944444444444642, 0.006944444444444642], error = 0.009820927516480105 |
| 34 | +Iteration 4: x = [1.4142156862745099, 1.4142156862745099], f(x) = [6.007304882871267e-6, 6.007304882871267e-6], error = 8.495612038666664e-6 |
| 35 | +Iteration 5: x = [1.4142135623746899, 1.4142135623746899], f(x) = [4.510614104447086e-12, 4.510614104447086e-12], error = 6.378971641140442e-12 |
| 36 | +Iteration 6: x = [1.4142135623730951, 1.4142135623730951], f(x) = [4.440892098500626e-16, 4.440892098500626e-16], error = 6.280369834735101e-16</code></pre><p>While the Halley's method takes 4 iterations to converge:</p><pre><code class="language-julia hljs">halley(f, [1., 1.], [2., 2.])</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">Iteration 1: x = [1.0, 1.0], f(x) = [-1.0, -1.0], error = 1.4142135623730951 |
| 37 | +Iteration 2: x = [1.4, 1.4], f(x) = [-0.04000000000000026, -0.04000000000000026], error = 0.05656854249492417 |
| 38 | +Iteration 3: x = [1.4142131979695431, 1.4142131979695431], f(x) = [-1.0306887576749801e-6, -1.0306887576749801e-6], error = 1.4576140196894333e-6 |
| 39 | +Iteration 4: x = [1.414213562373095, 1.414213562373095], f(x) = [-4.440892098500626e-16, -4.440892098500626e-16], error = 6.280369834735101e-16</code></pre></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../../">« Home</a><a class="docs-footer-nextpage" href="../../theory/">Theory »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option><option value="catppuccin-latte">catppuccin-latte</option><option value="catppuccin-frappe">catppuccin-frappe</option><option value="catppuccin-macchiato">catppuccin-macchiato</option><option value="catppuccin-mocha">catppuccin-mocha</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.8.0 on <span class="colophon-date" title="Tuesday 26 November 2024 23:06">Tuesday 26 November 2024</span>. Using Julia version 1.11.1.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html> |
0 commit comments