r"""Iterative solvers.
This module contains methods for solving iterative problems.
"""
from __future__ import division
__all__ = ['SimpleIteration', 'Broyden', 'ModifiedBroyden',
'LineSearch', 'extrapolate_range']
from zope.interface import implements
import numpy as np
from mmf.objects import process_vars
from mmf.objects import ClassVar
from mmf.solve import solver_interface, solver_utils
from mmf.solve.broyden_solvers import Broyden, ModifiedBroyden
from mmf.solve.solver_interface import IterationError, Mesg
_FINFO = np.finfo(float)
_TINY = _FINFO.tiny
_ABS_TOL = 1e-12
_REL_TOL = 1e-12
_MESGS = Mesg()
def _extrapolate(x, x0, x_controls, c_min=0.1):
r"""Provide a convex extrapolator."""
if x_controls is None:
return x
inds, vals = x_controls
inds = np.asarray(inds, dtype=int)
vals = np.asarray(vals)
if x0 is None:
# No dx, so we just enforce constraints.
x[inds] = vals
return x
dx = x - x0
# These are the values of c where v = c*dx
c = np.divide(vals - x0[inds], dx[inds])
# These are indices that must be constrained and the indices
# that cannot be constrained (otherwise there would be too
# short a step) and so who's values should be simply set.
constrained_inds = np.where(np.logical_and(
0 <= c, c < 1))[0]
set_inds = np.where(
np.logical_and(0 <= c, c < 1, c < c_min))[0]
if 0 == len(constrained_inds):
# Nothing
return x
c_step = max(c[constrained_inds].min(), c_min)
x = x0 + c_step*dx
x[inds[set_inds]] = vals[set_inds]
return x
[docs]class SimpleIteration(solver_interface.Solver):
r"""This solver simple performs direct iteration with a weight or
a simple weighted Newton's step (if a direction is supplied by the
problem).
Examples
--------
>>> from mmf.solve.solver_interface import Problem, BareState
>>> from mmf.objects import process_vars, Required, Computed
>>> class Sqrt(Problem):
... '''Iterative solution to sqrt(n).'''
... _state_vars = [
... ('n', Required, 'Array of n to compute sqrt(n)'),
... ('state0', Computed)]
... process_vars()
... def __init__(self, *varargin, **kwargs):
... self.n = np.array(self.n)
... self.state0 = BareState(x_=self.n)
... def F(self, x, iter_data={}, compute_dx=False,
... abs_tol=None, rel_tol=None, mesgs=_MESGS):
... '''Return one step of Newton's method `F(x)`.'''
... return x - (x*x - self.n)/x/2.0
>>> Sqrt.has_F
True
>>> solver = SimpleIteration(x_tol=0.01, G_tol=0.01, weight=0.5)
>>> problem = Sqrt(n=[1.0, 4.0, 9.0])
>>> state = problem.get_initial_state()
>>> state, niter = solver.solve(problem=problem, state=state)
>>> list(state.x_) # doctest: +ELLIPSIS
[1.0..., 2.0..., 3.0...]
"""
implements(solver_interface.ISolver)
_state_vars = [
('weight', 1.0, "New step `x1 + weight*dx`"),
('c_min', 0.1,
r"""Used by :meth:`_extrapolate`. When reducing the step
size to satisfy constraints, components are not
extrapolated to less than this value of their original
length to ensure that progress is always made."""),
]
process_vars()
[docs] def iterate(self, x, problem, solver_data,
x0=None, G0=None, dx0=None, F0=None,
mesgs=_MESGS):
r"""Perform a single iteration and return the triplet
(x_new, x_prev, G(x_prev)). The latter two entries are required
for convergence checking.
Mutate solver_data if used.
If there are any problems, then either an exception should be
raise, or a message should be appended to the list mesg.
"""
def extrapolate(x_controls):
return _extrapolate(x, x0, x_controls, self.c_min)
(x1, G1, dx1, F1) = self.problem_step(
problem=problem, solver_data=solver_data,
x=x, x0=x0, G0=G0, dx0=dx0, F0=F0,
compute_dx=True, compute_F=True,
extrapolate=extrapolate,
abs_tol=np.minimum(self.x_tol, self.G_tol),
rel_tol=self.x_tol_rel, mesgs=mesgs)
if dx1 is not None:
x_new = x1 + self.weight*dx1
elif F1 is not None:
x_new = self.weight*F1 + (1 - self.weight)*x1
else:
raise ValueError(
"%s solver requires either dx or F." %
(self.__class__.__name__))
return (x_new, (x1, G1, dx1, F1))
[docs]class LineSearch(solver_interface.Solver):
r"""This solver simple performs direct iteration with a line
search. If the problem supports it, then Newton's method is
applied.
Examples
--------
>>> from solver_interface import Problem, BareState
>>> from mmf.objects import process_vars, Required, Computed
>>> class Sqrt(Problem):
... '''Iterative solution to sqrt(n).'''
... _state_vars = [
... ('n', Required, 'Array of n to compute sqrt(n)'),
... ('state0', Computed)]
... process_vars()
... def __init__(self, *varargin, **kwargs):
... self.n = np.array(self.n)
... self.state0 = BareState(x_=self.n)
... def F(self, x, iter_data=None, compute_dx=False,
... abs_tol=None, rel_tol=None, mesgs=_MESGS):
... '''Return G(x) = x - F(x) where F(x) is one step of
... Newton's method.'''
... return x - (x*x - self.n)/x/2.0
>>> solver = LineSearch(x_tol=0.01, G_tol=0.01, weight=0.5)
>>> problem = Sqrt(n=[1.0, 4.0, 9.0])
>>> state = problem.get_initial_state()
>>> sol, niter = solver.solve(problem=problem, state=state)
>>> list(sol.x_) # doctest: +ELLIPSIS
[1.0..., 2.0..., 3.0...]
Now we try using Newton's method, but as implemented by LineSearch.
>>> class Sqrt(Problem):
... '''Iterative solution to sqrt(n).'''
... _state_vars = [
... ('n', Required, 'Array of n to compute sqrt(n)'),
... ('state0', Computed),
... ('has_dx', ClassVar(True))]
... process_vars()
... def __init__(self, *varargin, **kwargs):
... self.n = np.array(self.n)
... self.state0 = BareState(x_=self.n)
... def G(self, x, iter_data={}, compute_dx=False,
... abs_tol=None, rel_tol=None, mesgs=_MESGS):
... '''Return G(x) = x - F(x) where F(x) is one step of
... Newton's method.'''
... G = x*x - self.n
... dx = -G/x/2.0
... return G, dx
>>> solver = LineSearch(x_tol=0.01, G_tol=0.01, weight=0.5)
>>> problem = Sqrt(n=[1.0, 4.0, 9.0])
>>> state = problem.get_initial_state()
>>> sol, niter = solver.solve(problem=problem, state=state)
>>> list(sol.x_) # doctest: +ELLIPSIS
[1.0..., 2.0..., 3.0...]
"""
implements(solver_interface.ISolver)
_state_vars = [
('weight', 1.0,
r"""New step is a mixture with old step. (weight = 1.0 has
no old step, weight=0.0 does not move.)"""),
('ls_exact', False,
r"""If `True`, then the line search algorithm attempts to
find the exact minimum to tolerance
:attr:`ls_exact_rel_tol`"""),
('ls_strategy', 'default'),
('ls_exact_rel_tol', 1e-5),
('ls_avg_descent_factor', 1e-4,
r"""To ensure that we terminate, the line search algorithm
requires that each step descend with at least this slope,
otherwise the routine will terminate with an error
indicating that a local minimum was found rather than a
root."""),
('ls_min_scale_factor', 0.1,
r"""Minimum factor by which the step size is reduced in
the line search. See :func:`solver_utils.line_search`."""),
('ls_max_scale_factor', 0.5,
r"""Minimum factor by which the step size is reduced in
the line search. See
:func:`solver_utils.line_search`."""),
('ls_max_steps', 10,
r"""Maximum number of line search steps to take before
bailing out."""),
('min_step', 1e-12,
r"""To ensure termination, if the step size becomes smaller
than this, then the routine will terminate, raising an
:exc:`IterationError` exception."""),
('bad_step_factor', 1.1,
r"""For the first :attr:`N_bad_steps`, the algorithm
allows steps where the error increases by this factor.
This can allow and algorithm to "explore" a bit before
settling into convergence pattern. This can be useful in
a deep valley where requiring a strict decrease in
magnitude can prevent convergence.
"""),
('N_bad_steps', 0,
r"""For this number of steps, an increase in the error by
:attr:`bad_step_factor`is permitted."""),
('c_min', 0.1,
r"""Used by :meth:`_extrapolate`. When reducing the step
size to satisfy constraints, components are not
extrapolated to less than this value of their original
length to ensure that progress is always made."""),
]
process_vars()
[docs] def iterate(self, x, problem, solver_data,
x0=None, G0=None, dx0=None, F0=None, mesgs=_MESGS):
r"""Return `(x_new, (x_prev, G(x_prev), dx_prev, F(x_prev))`.
Mutate `solver_data` if used.
If there are any problems, then either an exception should be
raise, or a message should be appended to the list mesg.
Raises
------
:class:`solver_interface.IterationError`
If a step cannot be found because function does not decrease
fast enough, or step size is too small.
"""
if x0 is None:
# First step
x1 = x
else:
dx = x - x0
x_new = [None]
trust_dg = False
if problem.has_G:
def g(l):
x_new[0] = x0 + l*dx
_G = self.problem_G(
problem=problem, solver_data=solver_data,
x=x_new[0], compute_dx=False, mesgs=mesgs)[0]
if (hasattr(solver_data, 'Jinv_G') and
solver_data.Jinv_G is not None):
solver_data.Jinv_G.update(dx=l*dx, dg=_G - G0)
return np.linalg.norm(_G/problem.G_scale,
ord=self.norm_p)**2
g0 = np.linalg.norm(G0/problem.G_scale,
ord=self.norm_p)**2
if problem.has_dx:
trust_dg = True
else:
def g(l):
x_new[0] = x0 + l*dx
_F = self.problem_F(
problem=problem, solver_data=solver_data,
x=x_new[0], mesgs=mesgs)
_G = x_new[0] - _F
if (hasattr(solver_data, 'Jinv_F') and
solver_data.Jinv_F is not None):
solver_data.Jinv_F.update(dx=l*dx, dg=_F - F0)
return np.linalg.norm(_G/problem.G_scale,
ord=self.norm_p)**2
g0 = np.linalg.norm((x0 - F0)/problem.G_scale,
ord=self.norm_p)**2
if self.ls_exact:
l, g_l = solver_utils.exact_line_search(g=g, g0=g0,
rel_tol=self.ls_exact_rel_tol)
else:
try:
l, g_l = solver_utils.line_search(
g=g, g0=g0,
trust_dg=trust_dg,
avg_descent_factor=self.ls_avg_descent_factor,
min_step=self.min_step/np.linalg.norm(dx),
min_scale_factor=self.ls_min_scale_factor,
max_scale_factor=self.ls_max_scale_factor,
max_steps=self.ls_max_steps,
mesgs=mesgs)
except solver_utils.NoDescentError, err:
mesgs.warn(str(err) + "\nTrying without trusting dg.")
# Try once more without trusting dg
l, g_l = solver_utils.line_search(
g=g, g0=g0,
trust_dg=False,
avg_descent_factor=self.ls_avg_descent_factor,
min_step=self.min_step/np.linalg.norm(dx),
min_scale_factor=self.ls_min_scale_factor,
max_scale_factor=self.ls_max_scale_factor,
max_steps=self.ls_max_steps,
mesgs=mesgs)
x1 = x_new[0]
def extrapolate(x_controls):
return _extrapolate(x1, x0, x_controls, self.c_min)
compute_G = compute_F = False
compute_dx = True
if not problem.has_dx:
compute_F = True
x, G, dx, F = self.problem_step(
problem=problem, solver_data=solver_data,
x=x1, x0=x0, G0=G0, dx0=dx0, F0=F0,
compute_G=compute_G, compute_dx=compute_dx,
compute_F=compute_F,
extrapolate=extrapolate,
abs_tol=np.minimum(self.G_tol, self.x_tol),
rel_tol=self.x_tol_rel,
mesgs=mesgs)
if dx is not None:
x_new = x + dx
elif F is not None:
x_new = ((1.0 - self.weight)*x + self.weight*F)
else:
raise ValueError(
"%s solver requires either dx or F." %
(self.__class__.__name__))
dx = x_new - x
return (x_new, (x, G, dx, F))