Optimization
============

Root-finding with the secant method (``secant``)
------------------------------------------------

The function ``secant`` locates a root of a given function using the secant method. A simple example use of the secant method is to compute pi as the root of sin(*x*) closest to *x* = 3::

    >>> from mpmath import *
    >>> mp.dps = 30
    >>> print secant(sin, 3)
    3.14159265358979323846264338328

The secant method can be used to find complex roots of analytic functions, although it must in that case generally be given a nonreal starting value (or else it will never leave the real line)::

    >>> mp.dps = 15
    >>> print secant(lambda x: x**3 + 2*x + 1, j)
    (0.226698825758202 + 1.46771150871022j)

A good initial guess for the location of the root is required for the method to be effective, so it is somewhat more appropriate to think of the secant method as a root-polishing method than a root-finding method. When the rough location of the root is known, the secant method can be used to refine it to very high precision in only a few steps. If the root is a first-order root, only roughly log(prec) iterations are required. (The secant method is far less efficient for double roots.) It may be worthwhile to compute the initial approximation to a root using a machine precision solver (for example using one of SciPy's many solvers), and then refining it to high precision using mpmath's ``secant`` method.

Neat examples
.............

A nice application is to compute nontrivial roots of the Riemann zeta function with many digits (good initial values are needed for convergence)::

    >>> mp.dps = 30
    >>> print secant(zeta, 0.5+14j)
    (0.5 + 14.1347251417346937904572519836j)

The secant method can also be used as an optimization algorithm, by passing it a derivative of a function. The following example locates the positive minimum of the gamma function::

    >>> mp.dps = 20
    >>> print secant(lambda x: diff(gamma, x), 1)
    1.4616321449683623413

Finally, a useful application is to compute inverse functions, such as the Lambert W function which is the inverse of *w* exp(*w*), given the first term of the solution's asymptotic expansion as the initial value::

    >>> def lambert(x):
    ...     return secant(lambda w: w*exp(w) - x, log(1+x))
    ...
    >>> mp.dps = 15
    >>> print lambert(1)
    0.567143290409784
    >>> print lambert(1000)
    5.2496028524016

Options
.......

Strictly speaking, the secant method requires two initial values. By default, you only have to provide the first point ``x0``; ``secant`` automatically sets the second point (somewhat arbitrarily) to ``x0 + 1/4``. Manually providing also the second point can help in some cases if ``secant`` fails to converge.

By default, ``secant`` performs a maximum of 20 steps, which can be increased or decreased using the ``maxsteps`` keyword argument. You can pass ``secant`` the option ``verbose=True`` to show detailed progress.

Finding all roots of a polynomial (``polyroots``)
-------------------------------------------------

The function ``polyroots`` computes all *n* real or complex roots of an *n*-th degree polynomial. It will for example successfully compute the two real roots of 3x^2 - 7x + 2

    >>> mp.dps = 15
    >>> for r in polyroots([2,-7,3]):
    ...     print r
    ...
    0.333333333333333
    2.0

or the three roots of x^3 - x^2 - 14x + 24:

    >>> polyroots([24,-14,-1,1])
    [mpf('-4.0'), mpf('2.0'), mpf('3.0')]

The following example computes all the 5th roots of unity; i.e. the roots of x^5 - 1:

    >>> mp.dps = 20
    >>> for r in polyroots([-1, 0, 0, 0, 0, 1]):
    ...     print r
    ...
    1.0
    (-0.8090169943749474241 + 0.58778525229247312917j)
    (-0.8090169943749474241 - 0.58778525229247312917j)
    (0.3090169943749474241 + 0.95105651629515357212j)
    (0.3090169943749474241 - 0.95105651629515357212j)

Although all roots are internally calculated using complex arithmetic, any root found to have an imaginary part smaller than the estimated numerical error is truncated to a real number. Real roots are placed first in the returned list, sorted by value. The remaining complex numbers are sorted by real their parts so that conjugate roots end up next to each other.

Provided there are no repeated roots, ``polyroots`` can typically compute all roots with high precision:

    >>> mp.dps = 70
    >>> for r in polyroots([1, 0, -10, 0, 1]):
    ...     print r
    ...
    -3.146264369941972342329135065715570445512477129187328701232486717442666
    -0.3178372451957822447257576172961742883731333784334325548791272414612005
    0.3178372451957822447257576172961742883731333784334325548791272414612005
    3.146264369941972342329135065715570445512477129187328701232486717442666
    >>>
    >>> print sqrt(3) + sqrt(2)
    3.146264369941972342329135065715570445512477129187328701232486717442665
    >>> print sqrt(3) - sqrt(2)
    0.3178372451957822447257576172961742883731333784334325548791272414612005

Ill-conditioned polynomials
...........................

If a root with multiplicity M is present, it can typically only be computed accurately to dps/M digits (an sometimes less, due to slow convergence). For example, a triple root can be expected to be computed to roughly 5 accurate digits at standard 15-digit precision:

    >>> mp.dps = 15
    >>> for r in polyroots([1,3,3,1]):
    ...     print r
    ...
    (-1.00000007263211 - 1.9031226565599e-7j)
    (-0.999999852005337 + 2.06642147088233e-7j)
    (-1.00000162463939 - 1.50181300875847e-6j)

An example of an extremely ill-conditioned polynomial is Wilkinson's polynomial which has roots at the integers 1 through 20.

    >>> W = [2432902008176640000, -8752948036761600000, 13803759753640704000,
    ...     -12870931245150988800, 8037811822645051776, -3599979517947607200,
    ...     1206647803780373360, -311333643161390640, 63030812099294896,
    ...     -10142299865511450, 1307535010540395, -135585182899530,
    ...     11310276995381, -756111184500, 40171771630, -1672280820, 53327946,
    ...     -1256850, 20615, -210, 1]
    ...
    >>> roots = polyroots(W)
    >>> for r in roots:
    ...     print r
    ...
    1.0
    2.0
    3.0
    4.00000000000184
    4.99999999994666
    6.00000000011394
    7.00000001128901
    7.99999998684772
    8.99999991431934
    9.99999897717367
    10.9999998763699
    11.9999947045394
    13.0000024841578
    13.9999973809886
    15.0000157682936
    16.000000498559
    17.0000017669947
    18.0000008200121
    18.999999711431
    20.0000000125722

All roots have been separated out, but ``polyroots`` fails to converge to high accuracy. The roots can however effectively be polished using the secant method:

    >>> for r in roots:
    ...     print secant(extraprec(50)(lambda x: polyval(W, x)), r)
    ...
    1.0
    2.0
    3.0
    4.0
    5.0
    6.0
    7.0
    8.0
    9.0
    10.0
    11.0
    12.0
    13.0
    14.0
    15.0
    16.0
    17.0
    18.0
    19.0
    20.0

The extra precision ensures that the evaluation of the polynomial is accurate.