Skip to content

polykin.math.fixpoint¤

fixpoint_wegstein ¤

fixpoint_wegstein(
    g: Callable[[FloatVector], FloatVector],
    x0: FloatVector,
    *,
    wait: int = 1,
    qmin: float = -5.0,
    qmax: float = 0.0,
    tolx: float = 1e-06,
    sclx: FloatVector | None = None,
    maxiter: int = 50,
    callback: (
        Callable[
            [int, FloatVector, FloatVector],
            tuple[bool, bool],
        ]
        | None
    ) = None
) -> VectorRootResult

Find the solution of a N-dimensional fixed-point problem using the bounded Wegstein acceleration method.

The bounded Wegstein acceleration method is an extrapolation technique to accelerate the convergence of fixed-point iterations. For N-dimensional problems, each component of the fixed-point vector is treated separately according to:

\[ x_{k+1} = x_k + (1 - q_k) \left( g(x_k) - x_k \right) \]

where \(q_{min} \leq q_k \leq q_{max}\) is the acceleration parameter determined by:

\[\begin{aligned} q_k &= \frac{s_k}{s_k - 1} \\ s_k &= \frac{g(x_k) - g(x_{k-1})}{x_k - x_{k-1}} \end{aligned}\]

When \(q=0\), the Wegstein method is equivalent to the standard fixed-point iteration. When \(q<0\), the convergence is accelerated, and when \(0<q<1\) the convergence is dampened.

References

  • J.H. Wegstein, "Accelerating convergence of iterative processes", Communications of the ACM, 1(6): 9-13, 1958.
PARAMETER DESCRIPTION
g

Fixed-point mapping defining the problem g(x) = x.

TYPE: Callable[[FloatVector], FloatVector]

x0

Initial guess.

TYPE: FloatVector

wait

Number of direct substitution iterations before the first acceleration iteration.

TYPE: int DEFAULT: 1

qmin

Minimum value for the acceleration parameter.

TYPE: float DEFAULT: -5.0

qmax

Maximum value for the acceleration parameter.

TYPE: float DEFAULT: 0.0

tolx

Absolute tolerance for x value. The algorithm will terminate when ||sclx*(g(x) - x)||∞ <= tolx.

TYPE: float DEFAULT: 1e-06

sclx

Positive scaling factors for the components of x. Ideally, these should be chosen so that sclx*x is of order 1 near the solution for all components. By default, scaling is determined automatically from x0.

TYPE: FloatVector | None DEFAULT: None

maxiter

Maximum number of iterations.

TYPE: int DEFAULT: 50

callback

Optional callback with signature callback(niter, x, fx)->(stop, success) called at each iteration. If stop is True, the iteration is terminated. If success is True, the optimization is considered successful.

TYPE: Callable[[int, FloatVector, FloatVector], tuple[bool, bool]] | None DEFAULT: None

RETURNS DESCRIPTION
VectorRootResult

Dataclass with root solution results.
See Also
  • fixpoint_anderson: Alternative method better suited for problems with coupling between components.
  • fixpoint_damped: Alternative method for problems with weak coupling between components.
  • fixpoint_dem: Alternative method for problems with weak coupling between components.

Examples:

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

>>> from polykin.math import fixpoint_wegstein
>>> 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_wegstein(g, x0=np.array([0.0, 0.0]), qmax=0.5)
>>> 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/wegstein.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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
def fixpoint_wegstein(
    g: Callable[[FloatVector], FloatVector],
    x0: FloatVector,
    *,
    wait: int = 1,
    qmin: float = -5.0,
    qmax: float = 0.0,
    tolx: float = 1e-6,
    sclx: FloatVector | None = None,
    maxiter: int = 50,
    callback: Callable[[int, FloatVector, FloatVector], tuple[bool, bool]] | None = None,
) -> VectorRootResult:
    r"""Find the solution of a N-dimensional fixed-point problem using the
    bounded Wegstein acceleration method.

    The bounded Wegstein acceleration method is an extrapolation technique to
    accelerate the convergence of fixed-point iterations. For N-dimensional
    problems, each component of the fixed-point vector is treated separately
    according to:

    $$ x_{k+1} = x_k + (1 - q_k) \left( g(x_k) - x_k \right) $$

    where $q_{min} \leq q_k \leq q_{max}$ is the acceleration parameter
    determined by:

    \begin{aligned}
    q_k &= \frac{s_k}{s_k - 1} \\
    s_k &= \frac{g(x_k) - g(x_{k-1})}{x_k - x_{k-1}}
    \end{aligned}

    When $q=0$, the Wegstein method is equivalent to the standard fixed-point
    iteration. When $q<0$, the convergence is accelerated, and when $0<q<1$ the
    convergence is dampened.

    **References**

    * J.H. Wegstein, "Accelerating convergence of iterative processes",
      Communications of the ACM, 1(6): 9-13, 1958.

    Parameters
    ----------
    g : Callable[[FloatVector], FloatVector]
        Fixed-point mapping defining the problem `g(x) = x`.
    x0 : FloatVector
        Initial guess.
    wait : int
        Number of direct substitution iterations before the first acceleration
        iteration.
    qmin : float
        Minimum value for the acceleration parameter.
    qmax : float
        Maximum value for the acceleration parameter.
    tolx : float
        Absolute tolerance for `x` value. The algorithm will terminate when
        `||sclx*(g(x) - x)||∞ <= tolx`.
    sclx : FloatVector | None
        Positive scaling factors for the components of `x`. Ideally, these
        should be chosen so that `sclx*x` is of order 1 near the solution for
        all components. By default, scaling is determined automatically from `x0`.
    maxiter : int
        Maximum number of iterations.
    callback : Callable[[int, FloatVector, FloatVector], tuple[bool, bool]] | None
        Optional callback with signature `callback(niter, x, fx)->(stop, success)` called
        at each iteration. If `stop` is `True`, the iteration is terminated. If `success`
        is `True`, the optimization is considered successful.

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

    See Also
    --------
    * [`fixpoint_anderson`](fixpoint_anderson.md):
      Alternative method better suited for problems with coupling between components.
    * [`fixpoint_damped`](fixpoint_damped.md):
      Alternative method for problems with weak coupling between components.
    * [`fixpoint_dem`](fixpoint_dem.md):
      Alternative method for problems with weak coupling between components.

    Examples
    --------
    Find the solution of a 2D fixed-point function.
    >>> from polykin.math import fixpoint_wegstein
    >>> 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_wegstein(g, x0=np.array([0.0, 0.0]), qmax=0.5)
    >>> print(f"x = {sol.x}")
    x = [0.97458605 1.93830731]
    >>> print(f"g(x) = {g(sol.x)}")
    g(x) = [0.97458605 1.93830731]
    """
    method = "Wegstein fixed-point"
    success = False
    message = ""
    nfeval = 0

    sclx = np.abs(sclx) if sclx is not None else scalex(x0)

    x = x0.copy()
    n = x.size
    wait = max(wait, 1)

    gx = np.full(n, np.nan)
    fx = np.full(n, np.nan)
    xm = np.full(n, np.nan)

    niter = 0

    for niter in range(1, maxiter + 1):
        gxm = gx
        gx = g(x)
        nfeval += 1
        fx = gx - x

        if callback is not None:
            stop, _success = callback(niter, x, fx)
            if stop:
                message = "Terminated by user callback."
                success = _success
                break

        if np.linalg.norm(sclx * fx, np.inf) <= tolx:
            message = "||sclx*(g(x) - x)||∞ ≤ tolx"
            success = True
            break

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

        if niter < wait + 1:
            xm = x
            x = gx
        else:
            Δx = x - xm
            Δg = gx - gxm
            s = np.zeros(n)
            mask_s = np.abs(Δx) > eps
            np.divide(Δg, Δx, out=s, where=mask_s)
            q = s / (s - 1)
            q = np.clip(q, qmin, qmax)
            xm = x
            x = x + (1 - q) * fx

    return VectorRootResult(method, success, message, nfeval, None, niter, x, fx, None)