Skip to content

polykin.math.fixpoint¤

fixpoint_anderson ¤

fixpoint_anderson(
    g: Callable[[FloatVector], FloatVector],
    x0: FloatVector,
    m: int = 3,
    xtol: float = 1e-06,
    maxiter: int = 50,
) -> VectorRootResult

Find the solution of a N-dimensional fixed-point function using Anderson acceleration.

PARAMETER DESCRIPTION
g

Identity function for the fixed-point problem, i.e. g(x) = x.

TYPE: Callable[[FloatVector], FloatVector]

x0

Initial guess.

TYPE: FloatVector

m

Number of previous steps (m >= 1) to use in the acceleration.

TYPE: int DEFAULT: 3

xtol

Absolute tolerance for x value. The algorithm will terminate when ||g(x_k) - x_k||∞ <= xtol.

TYPE: float DEFAULT: 1e-06

maxiter

Maximum number of iterations.

TYPE: int DEFAULT: 50

RETURNS DESCRIPTION
VectorRootResult

Dataclass with root solution results.

Examples:

Find the solution of a 2D fixed-point function.

>>> from polykin.math import fixpoint_anderson
>>> import numpy as np
>>> def g(x):
...     x1, x2 = x
...     g1 = 0.5*np.cos(x1) + 0.1*x2 + 0.5
...     g2 = np.sin(x2) - 0.2*x1 + 1.2
...     return np.array([g1, g2])
>>> sol = fixpoint_anderson(g, x0=np.array([0.0, 0.0]))
>>> print(f"x={sol.x}")
x=[0.97458605 1.93830731]
>>> print(f"g(x)={g(sol.x)}")
g(x)=[0.97458605 1.93830731]
Source code in src/polykin/math/fixpoint/anderson.py 
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def fixpoint_anderson(
        g: Callable[[FloatVector], FloatVector],
        x0: FloatVector,
        m: int = 3,
        xtol: float = 1e-6,
        maxiter: int = 50
) -> VectorRootResult:
    r"""Find the solution of a N-dimensional fixed-point function using Anderson
    acceleration.

    Parameters
    ----------
    g : Callable[[FloatVector], FloatVector]
        Identity function for the fixed-point problem, i.e. `g(x) = x`.
    x0 : FloatVector
        Initial guess.
    m : int, optional
        Number of previous steps (`m >= 1`) to use in the acceleration.
    xtol : float, optional
        Absolute tolerance for `x` value. The algorithm will terminate when
        `||g(x_k) - x_k||∞ <= xtol`.
    maxiter : int, optional
        Maximum number of iterations.

    Returns
    -------
    VectorRootResult
        Dataclass with root solution results.

    Examples
    --------
    Find the solution of a 2D fixed-point function.
    >>> from polykin.math import fixpoint_anderson
    >>> import numpy as np
    >>> def g(x):
    ...     x1, x2 = x
    ...     g1 = 0.5*np.cos(x1) + 0.1*x2 + 0.5
    ...     g2 = np.sin(x2) - 0.2*x1 + 1.2
    ...     return np.array([g1, g2])
    >>> sol = fixpoint_anderson(g, x0=np.array([0.0, 0.0]))
    >>> print(f"x={sol.x}")
    x=[0.97458605 1.93830731]
    >>> print(f"g(x)={g(sol.x)}")
    g(x)=[0.97458605 1.93830731]
    """

    nfeval = 0
    message = ""

    # Different ordering of arrays to optimize memory access
    ΔG = np.zeros((m, x0.size))
    ΔF = np.zeros((x0.size, m))

    g0 = g(x0)
    nfeval += 1
    f0 = g0 - x0

    if np.linalg.norm(f0, np.inf) <= xtol:
        message = "||g(x0) - x0||∞ <= xtol"
        return VectorRootResult(True, message, nfeval, 0, x0, f0)

    x = g0
    gx = g0
    fx = f0
    success = False
    for k in range(1, maxiter):

        mk = min(m, k)

        ΔG[:-1, :] = ΔG[1:, :]
        ΔF[:, :-1] = ΔF[:, 1:]

        ΔG[-1, :] = -gx
        ΔF[:, -1] = -fx

        gx = g(x)
        nfeval += 1
        fx = gx - x

        ΔG[-1, :] += gx
        ΔF[:, -1] += fx

        if np.linalg.norm(fx, np.inf) <= xtol:
            message = "||g(x_k) - x_k||∞ <= xtol"
            success = True
            break

        # There are methods to reuse the QR decomposition from the previous
        # iteration, but this is more complex to implement.
        try:
            Q, R = np.linalg.qr(ΔF[:, -mk:])
            gamma = np.linalg.lstsq(R, Q.T @ fx, rcond=None)[0]
        except np.linalg.LinAlgError:
            message = "Error in QR decomposition or least-squares solution."
            break

        if k < maxiter:
            x = gx - np.dot(gamma, ΔG[-mk:, :])

    else:
        message = f"Maximum number of iterations ({maxiter}) reached."

    return VectorRootResult(success, message, nfeval, k+1, x, fx)