improvements to fixed-point solvers

Author: yoelcortes

flexsolve/__init__.py

@@ -6,26 +6,33 @@
 """
 from . import open_solvers
 from . import bounded_solvers
-from . import iterative_solvers
+from . import fixed_point_solvers
+from . import line_search
+from . import numerical_analysis
 from . import problem_list
 from . import profiler
 from . import problem
 from . import utils
 
-__all__ = (*open_solvers.__all__,
-           *bounded_solvers.__all__,
-           *iterative_solvers.__all__,
-           *problem_list.__all__,
-           *problem.__all__,
-           *profiler.__all__,
-           'utils',
+__all__ = (
+    *open_solvers.__all__,
+    *bounded_solvers.__all__,
+    *fixed_point_solvers.__all__,
+    *line_search.__all__,
+    *numerical_analysis.__all__,
+    *problem_list.__all__,
+    *problem.__all__,
+    *profiler.__all__,
+    'utils',
 )
 
 from .open_solvers import *
 from .bounded_solvers import *
-from .iterative_solvers import *
+from .fixed_point_solvers import *
+from .line_search import *
+from .numerical_analysis import *
 from .problem_list import *
 from .problem import *
 from .profiler import *
 
-__version__ = '0.5.7'
+__version__ = '0.5.8'

flexsolve/iterative_solvers.py flexsolve/fixed_point_solvers.pyflexsolve/iterative_solvers.py renamed to flexsolve/fixed_point_solvers.py

@@ -17,17 +17,60 @@
            'conditional_aitken',
 ) 
 
+iteration_history = []
+
+class CachedFunction:
+    __slots__ = ('f', 'inputs', 'outputs')
+    
+    def __init__(self, f, inputs=None, outputs=None):
+        self.f = f
+        self.inputs = [] if inputs is None else inputs
+        self.outputs = [] if inputs is None else outputs
+        
+    def __call__(self, x, *args, **kwargs):
+        self.inputs.append(x)
+        y = self.f(x, *args, **kwargs)
+        self.outputs.append(y)
+        return y
+
+def GPiteration():
+    pass
+
+def GPacceleration(
+        f, x, xtol=5e-8, args=(), maxiter=50, memory=10, checkiter=True,
+        checkconvergence=True, bounds=None, convergenceiter=0, subset=0, *, rtol=0
+    ):
+    """Iterative fixed-point solver with GP acceleration."""
+    x0 = x1 = x
+    errors = np.zeros(convergenceiter)
+    f = CachedFunction(f)
+    fixedpoint_converged = utils.fixedpoint_converged
+    for iter in range(maxiter):
+        x1 = GPiteration(f, x0, args, memory, bounds)
+        e = np.abs(x1 - x0)
+        if fixedpoint_converged(x0, e, xtol, rtol, subset): return x1
+        if convergenceiter:
+            mean = utils.mean(e)
+            if iter > convergenceiter and mean > errors.mean():
+                if checkconvergence: utils.raise_convergence_error()
+                else: return x1
+            errors = np.roll(errors, shift=1)
+            errors[-1] = mean
+        x0 = x1
+    if checkiter: utils.raise_iter_error()
+    return x1
+
 @register_jitable(cache=True)
 def fixed_point(f, x, xtol=5e-8, args=(), maxiter=50, checkiter=True,
-                checkconvergence=True, convergenceiter=0, subset=0):
+                checkconvergence=True, convergenceiter=0, subset=0, *, rtol=0):
     """Iterative fixed-point solver."""
     x0 = x1 = x
     errors = np.zeros(convergenceiter)
     fixedpoint_converged = utils.fixedpoint_converged
     for iter in range(maxiter):
         x1 = f(x0, *args)
         e = np.abs(x1 - x0)
-        if fixedpoint_converged(e, xtol, subset): return x1
+        if fixedpoint_converged(x0, e, xtol, rtol, subset): return x1
         if convergenceiter:
             mean = utils.mean(e)
             if iter > convergenceiter and mean > errors.mean():
@@ -51,7 +94,8 @@ def conditional_fixed_point(f, x):
 
 @register_jitable(cache=True)
 def wegstein(f, x, xtol=5e-8, args=(), maxiter=50, checkiter=True,
-             checkconvergence=True, convergenceiter=0, subset=0):
+             checkconvergence=True, convergenceiter=0, subset=0, *, 
+             lb=-float('inf'), ub=float('inf'), rtol=0, exp=1.0):
     """Iterative Wegstein solver."""
     errors = np.zeros(convergenceiter)
     x0 = x
@@ -65,9 +109,9 @@ def wegstein(f, x, xtol=5e-8, args=(), maxiter=50, checkiter=True,
             x1 = g0
             g1 = f(x1, *args)
         e = np.abs(g1 - x1)
-        if fixedpoint_converged(e, xtol, subset): return g1
+        if fixedpoint_converged(x1, e, xtol, rtol, subset): return g1
         x0 = x1
-        x1 = wegstein_iter(x1, dx, g1, g0)
+        x1 = wegstein_iter(x1, dx, g1, g0, lb, ub, exp)
         g0 = g1
         if convergenceiter:
             mean = utils.mean(e)
@@ -80,7 +124,7 @@ def wegstein(f, x, xtol=5e-8, args=(), maxiter=50, checkiter=True,
     return x1
 
 @register_jitable(cache=True)
-def conditional_wegstein(f, x):
+def conditional_wegstein(f, x, lb=-float('inf'), ub=float('inf'), exp=0.5):
     """Conditional iterative Wegstein solver."""
     x0 = x
     g0, condition = f(x0)
@@ -94,13 +138,13 @@ def conditional_wegstein(f, x):
         g1 = g1
         dx = x1-x0
         x0 = x1
-        x1 = wegstein_iter(x1, dx, g1, g0)
+        x1 = wegstein_iter(x1, dx, g1, g0, lb, ub, exp)
         g0 = g1
     return x1
 
 @register_jitable(cache=True)
 def aitken(f, x, xtol=5e-8, args=(), maxiter=50, checkiter=True,
-           checkconvergence=True, convergenceiter=0, subset=0):
+           checkconvergence=True, convergenceiter=0, subset=0, *, rtol=0):
     """Iterative Aitken solver."""
     gg = x
     errors = np.zeros(convergenceiter)
@@ -113,10 +157,10 @@ def aitken(f, x, xtol=5e-8, args=(), maxiter=50, checkiter=True,
             g = f(x, *args)
         dxg = x - g
         e = np.abs(dxg)
-        if fixedpoint_converged(e, xtol, subset): return g
+        if fixedpoint_converged(x, e, xtol, rtol, subset): return g
         gg = f(g, *args)
         dgg_g = gg - g
-        if fixedpoint_converged(np.abs(dgg_g), xtol, subset): return gg
+        if fixedpoint_converged(g, np.abs(dgg_g), xtol, rtol, subset): return gg
         x = aitken_iter(x, gg, dxg, dgg_g)
         if convergenceiter:
             mean = utils.mean(e)

flexsolve/line_search.py

@@ -0,0 +1,22 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jul 16 13:18:51 2025
+
+@author: yoelr
+"""
+from scipy.optimize import fminbound
+
+__all__ = ('line_search',)
+
+def line_search(f, x, descent, t0, t1, ttol=1e-3, maxiter=20):
+    return fminbound(
+        lambda t: f(x + t * descent),
+        x1=t0,
+        x2=t1,
+        xtol=ttol,
+        maxfun=maxiter,
+    )
+    
+    
+    
+    

flexsolve/numerical_analysis.py

@@ -0,0 +1,139 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Jul 13 16:48:47 2025
+
+@author: yoelr
+"""
+import numpy as np
+from scipy.differentiate import jacobian
+from numpy.linalg import cond as matrix_condition_number
+from scipy.linalg import block_diag
+
+__all__ = (
+    'function_condition_number',
+    'independent_functions_condition_number',
+    'matrix_condition_number',
+)
+
+def norm(x):
+    return np.linalg.norm(x, ord=2) if np.ndim(x) else x
+
+def independent_functions_condition_number(*fxs, **jac_kwargs):
+    """
+    Estimate the condition number of a set of nonlinear functions at a 
+    given point using finite-difference Jacobian from SciPy.
+
+    Parameters
+    ----------
+    fxs : callable
+        Pairs of functios and points at which to evaluate the condition number.
+        The function is from ℝⁿ to ℝᵐ and must accept a 1D NumPy array.         
+    **jac_kwargs : dict, optional
+        Additional keyword arguments passed to `scipy.differentiate.jacobian`,
+        such as `initial_step`, `order`, `step_direction`, etc.
+
+    Returns
+    -------
+    kappa : float
+        The estimated condition number at point `x`. If `f(x)` has zero norm,
+        returns `np.inf`.
+
+    Notes
+    -----
+    The condition number is computed as:
+
+        κ(x) = ||J(x)|| * ||x|| / ||f(x)||
+
+    where J(x) is the Jacobian of `f` at `x`, and all norms are 2-norms 
+    (by default).
+    """
+    xs = []
+    fs = []
+    Js = []
+    def add(lst, value):
+        if np.ndim(value):
+            lst.extend(value)
+        else:
+            lst.append(value)
+            
+    for f, x in fxs:
+        x = np.asarray(x, dtype=float)
+        res = jacobian(f, x, **jac_kwargs)
+        fx = f(x)
+        J = res.df
+        add(xs, x)
+        add(fs, fx)
+        add(Js, J)
+        
+    J = block_diag(*Js)
+    x = np.array(xs)
+    fx = np.array(fs)
+    norm_x = norm(x)
+    norm_fx = norm(fx)
+    norm_J = norm(J)
+    if norm_fx == 0:
+        kappa = np.inf
+    else:
+        kappa = norm_J * norm_x / norm_fx
+    return kappa
+
+def function_condition_number(f, x, **jac_kwargs):
+    """
+    Estimate the condition number of a nonlinear function at a given point
+    using finite-difference Jacobian from SciPy.
+
+    Parameters
+    ----------
+    f : callable
+        A function from ℝⁿ to ℝᵐ. Must accept a 1D NumPy array of shape (n,)
+        and return a 1D array of shape (m,).
+    x : array_like, shape (n,)
+        Point at which to evaluate the condition number.
+    **jac_kwargs : dict, optional
+        Additional keyword arguments passed to `scipy.differentiate.jacobian`,
+        such as `initial_step`, `order`, `step_direction`, etc.
+
+    Returns
+    -------
+    kappa : float
+        The estimated condition number at point `x`. If `f(x)` has zero norm,
+        returns `np.inf`.
+
+    Notes
+    -----
+    The condition number is computed as:
+
+        κ(x) = ||J(x)|| * ||x|| / ||f(x)||
+
+    where J(x) is the Jacobian of `f` at `x`, and all norms are 2-norms 
+    (by default).
+    """
+    x = np.asarray(x, dtype=float)
+    # Compute numerical Jacobian using SciPy's differentiate API
+    res = jacobian(f, x, **jac_kwargs)
+    J = res.df  # The Jacobian matrix (m, n)
+    fx = f(x)
+    norm_x = norm(x)
+    norm_fx = norm(fx)
+    norm_J = norm(J)
+    if norm_fx == 0:
+        kappa = np.inf
+    else:
+        kappa = norm_J * norm_x / norm_fx
+    return kappa
+
+
+if __name__ == '__main__':
+    
+    def f0(x0):
+        return x0[0] ** 2 - x0[1] ** 0.5 + 5  
+    
+    def f1(x1):
+        return x1[0] ** 0.8 - x1[1] ** 3
+    
+    def f(x):
+        return np.array([f0(x[:2]), f1(x[2:])])
+    
+    x = np.array([2, 3, 0.5, -1])
+    kappa = function_condition_number(f, x)
+    kappa_mix = independent_functions_condition_number((f0, x[:2]), (f1, x[2:]))

flexsolve/problem.py

@@ -4,7 +4,7 @@
 
 @author: yrc2
 """
-from .iterative_solvers import wegstein, aitken
+from .fixed_point_solvers import wegstein, aitken
 from .profiler import Profiler
 
 __all__ = ('Problem',)