fixed_point_zhou.py

import numpy as np
from scipy.optimize.minpack import _lazywhere, _relerr

def fixed_point_zhou(func, x0, args=(), xtol=1E-14, maxiter=100000, p=6, rcond=1E-6):
    """
    Find a fixed point of the function.

    Given a function of one or more variables and a starting point, find a
    fixed-point of the function: i.e. where ``func(x0) == x0``.

    Parameters
    ----------
    func : function
        Function to evaluate.
    x0 : array_like
        Fixed point of function.
    args : tuple, optional
        Extra arguments to `func`.
    xtol : float, optional
        Convergence tolerance, defaults to 1E-14
    maxiter : int, optional
        Maximum number of iterations, defaults to 100000.
    p : int, optional, default=6
        Length of history to consider when constructing quasi-newton correction.
        A (p,p) matrix will be pseudo-inverted, so performance depends on finding
        a reasonable choice of p.
    rcond : float, optional, default=1E-6
        Used to discard poorly conditioned subspaces in the pseudoinverse.

    References
    ----------
    .. [1] Hua Zhou, David Alexander, Kenneth Lange. Dec. 2009.  Stat. Comp.

    Examples
    --------
    >>> from fixed_point_zhou import fixed_point_zhou
    >>> def func(x, c1, c2):
    ...    return np.sqrt(c1/(x+c2))
    >>> c1 = np.array([10,12.])
    >>> c2 = np.array([3, 5.])
    >>> fixed_point_zhou(func, np.array([1.2, 1.3]), args=(c1,c2))
    array([ 1.4920333 ,  1.37228132])
    """

    K = len(x0)
    U = np.zeros((K, p))  # History of F(xn) - xn, with most recent at U[:, -1]
    V = np.zeros((K, p))  # History of F(F(xn)) - F(xn), with most recent at V[:, -1]
    xn = x0.copy()

    for i in range(maxiter):
        xold = xn
        Fxn = func(xn, *args)
        FFxn = func(Fxn, *args)
        u = Fxn - xn
        v = FFxn - Fxn
        U = np.roll(U, -1, axis=1)  # Cycle the oldest entry into the -1 entry for replacement with the current value of u
        V = np.roll(V, -1, axis=1)  # Cycle the oldest entry into the -1 entry for replacement with the current value of v
        U[:, -1] = u
        V[:, -1] = v
        minv = np.linalg.pinv(U.T.dot(U) - U.T.dot(V), rcond=rcond)  # Be conservative about conditioning.  E.g. it's better to zero out poorly conditioned subspaces.
        #rhs are the rightmost terms in the paper equation immediately before eqn (2).  e.g. x^{n+1} = F(x^n) - [...]
        rhs = xn - Fxn
        rhs = U.T.dot(rhs)
        rhs = minv.dot(rhs)
        rhs = V.dot(rhs)
        xn = Fxn - rhs
        relerr = _lazywhere(xold != 0, (xn, xold), f=_relerr, fillvalue=xn)
        if np.all(np.abs(relerr) < xtol):
            break

    return xn