mici.solvers module#

Solvers for non-linear systems of equations for implicit integrators.

class mici.solvers.FixedPointSolver(*args, **kwargs)[source]#

Bases: Protocol

Solver for fixed point equation func(x) = x.

__call__(func, x0, **kwargs)[source]#

Solve fixed point equation.

Parameters:

  • func (ArrayFunction) – Function to solve for fixed point of.

  • x0 (ArrayLike) – Point to initialize solver at.

Returns:

Fixed point solved for.

Return type:

ArrayLike

class mici.solvers.ProjectionSolver(*args, **kwargs)[source]#

Bases: Protocol

Solver for projection on to manifold step in constrained integrator.

Solves an equation of the form

\[r(\lambda) = c(\Phi_{2,1}(t)(q, p + \partial c(q)^T \lambda)) = 0,\]

for the vector of Lagrange multipliers \(\lambda\) to project a point on to the manifold defined by the zero level set of a constraint function \(c\), with \(\Phi_{2,1}\) the flow map for the \(h_2\) Hamiltonian component for the system restricted to the position component output. The map \(\Phi_{2,1}\) is assumed to be linear in its second (momentum) argument.

__call__(state, state_prev, time_step, system, **kwargs)[source]#

Solve for projection on to manifold step.

Parameters:

  • state (ChainState) – Current chain state after unconstrained step.

  • state_prev (ChainState) – Previous chain state on manifold.

  • time_step (float) – Integrator time step for unconstrained step.

  • system (ConstrainedTractableFlowSystem) – Hamiltonian system constrained dynamics are being simulated for.

Returns:

Chain state after projection on to manifold.

Return type:

ChainState

mici.solvers.euclidean_norm(vct)[source]#

Calculate the Euclidean (L-2) norm of a vector.

mici.solvers.maximum_norm(vct)[source]#

Calculate the maximum (L-infinity) norm of a vector.

mici.solvers.solve_fixed_point_direct(func, x0, convergence_tol=1e-09, divergence_tol=10000000000.0, max_iters=100, norm=<function maximum_norm>)[source]#

Solve fixed point equation func(x) = x using direct iteration.

Parameters:

  • func (ArrayFunction) – Function to find fixed point of.

  • x0 (ArrayLike) – Initial state (function argument).

  • convergence_tol (float) – Convergence tolerance - solver successfully terminates when norm(func(x) - x) < convergence_tol.

  • divergence_tol (float) – Divergence tolerance - solver aborts if norm(func(x) - x) > divergence_tol on any iteration.

  • max_iters (int) – Maximum number of iterations before raising exception.

  • norm (ScalarFunction) – Norm to use to assess convergence.

Returns:

Solution to fixed point equation with norm(func(x) - x) < convergence_tol.

Raises:

ConvergenceError – If solver does not converge within max_iters iterations, diverges or encounters a ValueError during the iteration.

Return type:

ArrayLike

mici.solvers.solve_fixed_point_steffensen(func, x0, convergence_tol=1e-09, divergence_tol=10000000000.0, max_iters=100, norm=<function maximum_norm>)[source]#

Solve fixed point equation func(x) = x using Steffensen’s method.

Steffennsen’s method achieves quadratic convergence but at the cost of two function evaluations per iteration so for functions where convergence is achieved in a small number of iterations, direct iteration may be cheaper.

Parameters:

  • func (ArrayFunction) – Function to find fixed point of.

  • x0 (ArrayLike) – Initial state (function argument).

  • convergence_tol (float) – Convergence tolerance - solver successfully terminates when norm(func(x) - x) < convergence_tol.

  • divergence_tol (float) – Divergence tolerance - solver aborts if norm(func(x) - x) > divergence_tol on any iteration.

  • max_iters (int) – Maximum number of iterations before raising exception.

  • norm (ScalarFunction) – Norm to use to assess convergence.

Returns:

Solution to fixed point equation with norm(func(x) - x) < convergence_tol.

Raises:

ConvergenceError – If solver does not converge within max_iters iterations, diverges or encounters a ValueError during the iteration.

Return type:

ArrayLike

mici.solvers.solve_projection_onto_manifold_newton(state, state_prev, time_step, system, constraint_tol=1e-09, position_tol=1e-08, divergence_tol=10000000000.0, max_iters=50, norm=<function maximum_norm>)[source]#

Solve constraint equation using Newton’s method.

Requires re-evaluating both the constraint function system.constr and constraint Jacobian system.jacob_constr within the solver loop and computation of matrix decompositions on each iteration.

Solves an equation of the form

\[r(\lambda) = c(\Phi_{2,1}(t)(q, p + \partial c(q)^T \lambda)) = 0,\]

for the vector of Lagrange multipliers \(\lambda\) to project a point on to the manifold defined by the zero level set of a constraint function \(c\), with \(\Phi_{2,1}\) the flow map for the \(h_2\) Hamiltonian component for the system restricted to the position component output.

The Jacobian of the residual function \(r\) is

\[\partial r(\lambda) = \partial c(\Phi_{2,1}(t)(q, p + \partial c(q)^T \lambda)) \partial_2 (\Phi_{2,1}(t))(q, p + \partial c(q)^T \lambda) \partial c(q)^T\]

where \(\partial_2 (\Phi_{2,1}(t))\) is the Jacobian of the (position restricted) flow-map \(\Phi_{2,1}\) with respect to its second (momentum) argument. We assume here that \(\Phi_{2,1}(t)\) is linear in its second (momentum) argument, such that \(\partial_2 (\Phi_{2,1}(t))(q, p + \partial c(q)^T \lambda)\) is constant with respect to \(\lambda\).

The full Newton update is

\[\lambda' = \lambda - \partial r(\lambda)^{-1} r(\lambda)\]

which requires evaluating \(\partial c\) on each iteration and solving a linear system in the residual Jacobian \(\partial r(\lambda)\).

Parameters:

  • state (ChainState) – Post h2_flow update state to project.

  • state_prev (ChainState) – Previous state in co-tangent bundle before h2_flow update which defines the co-tangent space to perform projection in.

  • time_step (float) – Integrator time step used in h2_flow update.

  • system (ConstrainedEuclideanMetricSystem) – Hamiltonian system defining h2_flow and constr functions used to define constraint equation to solve.

  • constraint_tol (float) – Convergence tolerance in constraint space. Iteration will continue until norm(constr(pos)) < constraint_tol where pos is the position at the current iteration.

  • position_tol (float) – Convergence tolerance in position space. Iteration will continue until norm(delta_pos) < position_tol where delta_pos is the change in the position in the current iteration.

  • divergence_tol (float) – Divergence tolerance - solver aborts if norm(constr(pos)) > divergence_tol on any iteration where pos is the position at the current iteration and raises mici.errors.ConvergenceError.

  • max_iters (int) – Maximum number of iterations to perform before aborting and raising mici.errors.ConvergenceError.

  • norm (ScalarFunction) – Norm to use to test for convergence.

Returns:

Updated state object with position component satisfying constraint equation to within constraint_tol, i.e. norm(system.constr(state.pos)) < constraint_tol.

Raises:

ConvergenceError – If solver does not converge within max_iters iterations, diverges or encounters a ValueError during the iteration.

Return type:

ChainState

Solve constraint equation using Newton’s method with backtracking line-search.

Requires re-evaluating both the constraint function system.constr and constraint Jacobian system.jacob_constr within the solver loop and computation of matrix decompositions on each iteration.

Solves an equation of the form

\[r(\lambda) = c(\Phi_{2,1}(t)(q, p + \partial c(q)^T \lambda)) = 0,\]

for the vector of Lagrange multipliers \(\lambda\) to project a point on to the manifold defined by the zero level set of a constraint function \(c\), with \(\Phi_{2,1}\) the flow map for the \(h_2\) Hamiltonian component for the system restricted to the position component output.

The Jacobian of the residual function \(r\) is

\[\partial r(\lambda) = \partial c(\Phi_{2,1}(t)(q, p + \partial c(q)^T \lambda)) \partial_2 (\Phi_{2,1}(t))(q, p + \partial c(q)^T \lambda) \partial c(q)^T\]

where \(\partial_2 (\Phi_{2,1}(t))\) is the Jacobian of the (position restricted) flow-map \(\Phi_{2,1}\) with respect to its second (momentum) argument. We assume here that \(\Phi_{2,1}(t)\) is linear in its second (momentum) argument, such that \(\partial_2 (\Phi_{2,1}(t))(q, p + \partial c(q)^T \lambda)\) is constant with respect to \(\lambda\).

The scaled Newton update is

\[\lambda'(\alpha) = \lambda - \alpha \partial r(\lambda)^{-1} r(\lambda)\]

which requires evaluating \(\partial c\) on each iteration and solving a linear system in the residual Jacobian \(\partial r(\lambda)\).

The step size \(\alpha \in [0, 1]\) is initialised at \(\alpha = 1\) with a backtracking line search performed by multiplying \(\alpha\) by 0.5 until \(|r(\lambda'(\alpha))| < |r(\lambda)\).

Parameters:

  • state (ChainState) – Post h2_flow update state to project.

  • state_prev (ChainState) – Previous state in co-tangent bundle before h2_flow update which defines the co-tangent space to perform projection in.

  • time_step (float) – Integrator time step used in h2_flow update.

  • system (ConstrainedEuclideanMetricSystem) – Hamiltonian system defining h2_flow and constr functions used to define constraint equation to solve.

  • constraint_tol (float) – Convergence tolerance in constraint space. Iteration will continue until norm(constr(pos)) < constraint_tol where pos is the position at the current iteration.

  • position_tol (float) – Convergence tolerance in position space. Iteration will continue until norm(delta_pos) < position_tol where delta_pos is the change in the position in the current iteration.

  • divergence_tol (float) – Divergence tolerance - solver aborts if norm(constr(pos)) > divergence_tol on any iteration where pos is the position at the current iteration and raises mici.errors.ConvergenceError.

  • max_iters (int) – Maximum number of iterations to perform before aborting and raising mici.errors.ConvergenceError.

  • max_line_search_iters (int) – Maximum number of ‘inner’ line search iterations to perform to try to find step size along search direction which decreases residual norm.

  • norm (ScalarFunction) – Norm to use to test for convergence.

Returns:

Updated state object with position component satisfying constraint equation to within constraint_tol, i.e. norm(system.constr(state.pos)) < constraint_tol.

Raises:

ConvergenceError – If solver does not converge within max_iters iterations, diverges or encounters a ValueError during the iteration.

Return type:

ChainState

mici.solvers.solve_projection_onto_manifold_quasi_newton(state, state_prev, time_step, system, constraint_tol=1e-09, position_tol=1e-08, divergence_tol=10000000000.0, max_iters=50, norm=<function maximum_norm>)[source]#

Solve constraint equation using symmetric quasi-Newton method.

Only requires re-evaluating the constraint function system.constr within the solver loop and no recomputation of matrix decompositions on each iteration.

Solves an equation of the form

\[r(\lambda) = c(\Phi_{2,1}(t)(q, p + \partial c(q)^T \lambda)) = 0,\]

for the vector of Lagrange multipliers \(\lambda\) to project a point on to the manifold defined by the zero level set of a constraint function \(c\), with \(\Phi_{2,1}\) the flow map for the \(h_2\) Hamiltonian component for the system restricted to the position component output.

The Jacobian of the residual function \(r\) is

\[\partial r(\lambda) = \partial c(\Phi_{2,1}(t)(q, p + \partial c(q)^T \lambda)) \partial_2 (\Phi_{2,1}(t))(q, p + \partial c(q)^T \lambda) \partial c(q)^T\]

where \(\partial_2 (\Phi_{2,1}(t))\) is the Jacobian of the (position restricted) flow-map \(\Phi_{2,1}\) with respect to its second (momentum) argument. We assume here that \(\Phi_{2,1}(t)\) is linear in its second (momentum) argument, such that \(\partial_2 (\Phi_{2,1}(t))(q, p + \partial c(q)^T \lambda)\) is constant with respect to \(\lambda\), and henceforth use the shorthand \(\partial_2 (\Phi_{2,1}(t))\) to refer to this matrix.

The full Newton update is

\[\lambda' = \lambda - \partial r(\lambda)^{-1} r(\lambda)\]

which requires evaluating \(\partial c\) on each iteration and solving a linear system in the residual Jacobian \(\partial r(\lambda)\).

The symmetric quasi-Newton iteration instead uses the approximation

\[\partial c(\Phi_{2,1}(t)(q, p + \partial c(q)^T\lambda)) \partial_2 (\Phi_{2,1}(t)) \partial c(q)^T \approx \partial c(q) \partial_2 (\Phi_{2,1}(t)) \partial c(q)^T\]

with the corresponding update

\[\lambda' = \lambda - (\partial c(q) \partial_2 (\Phi_{2,1}(t)) \partial c(q)^T)^{-1} r(\lambda)\]

allowing a previously computed decomposition of the matrix

\[\partial c(q) \partial_2 (\Phi_{2,1}(t)) \partial c(q)^T,\]

to be used to solve the linear system in each iteration with no requirement to evaluate \(\partial c\) (system.jacob_constr) on each iteration.

Parameters:

  • state (ChainState) – Post h2_flow update state to project.

  • state_prev (ChainState) – Previous state in co-tangent bundle before h2_flow update which defines the co-tangent space to perform projection in.

  • time_step (float) – Integrator time step used in h2_flow update.

  • system (ConstrainedEuclideanMetricSystem) – Hamiltonian system defining h2_flow and constr functions used to define constraint equation to solve.

  • constraint_tol (float) – Convergence tolerance in constraint space. Iteration will continue until norm(constr(pos)) < constraint_tol where pos is the position at the current iteration.

  • position_tol (float) – Convergence tolerance in position space. Iteration will continue until norm(delta_pos) < position_tol where delta_pos is the change in the position in the current iteration.

  • divergence_tol (float) – Divergence tolerance - solver aborts if norm(constr(pos)) > divergence_tol on any iteration where pos is the position at the current iteration and raises mici.errors.ConvergenceError.

  • max_iters (int) – Maximum number of iterations to perform before aborting and raising mici.errors.ConvergenceError.

  • norm (ScalarFunction) – Norm to use to test for convergence.

Returns:

Updated state object with position component satisfying constraint equation to within constraint_tol, i.e. norm(system.constr(state.pos)) < constraint_tol.

Raises:

ConvergenceError – If solver does not converge within max_iters iterations, diverges or encounters a ValueError during the iteration.

Return type:

ChainState