soupy.optimization

soupy.optimization.bfgs

class soupy.optimization.bfgs.BFGS(cost_functional, parameters=hippylib.ParameterList)

Bases: object

Implement BFGS technique with backtracking inexact line search and damped updating See Nocedal & Wright (06), ch.6.2, ch.7.3, ch.18.3 The user must provide a ControlCostFunctional that describes the forward problem, cost functional, and all the derivatives for the gradient and the Hessian.

Initialize the BFGS solver. Type BFGS_ParameterList().showMe() for default parameters and their description

Parameters:

• cost_functional – The cost functional

• parameters – The parameters for the BFGS solver

solve(z, H0inv=<soupy.optimization.bfgs.RescaledIdentity object>, box_bounds=None, constraint_projection=None)

Solve the constrained optimization problem with initial guess z

Parameters:

• z (dolfin.Vector) – The initial guess

• H0inv – Initial approximation of the inverse of the Hessian. Has optional method update(x) that will update the operator

• box_bounds (list) – Bound constraint. A list with two entries (min and max). Can be either a scalar value or a dolfin.Vector of the same size as z

• constraint_projection (ProjectableConstraint) – Alternative projectable constraint

Returns:

The optimization solution z and a dictionary of information

Note

The input z will be overwritten

termination_reasons = ['Maximum number of Iteration reached', 'Norm of the gradient less than tolerance', 'Maximum number of backtracking reached', 'Norm of (g, da) less than tolerance']

soupy.optimization.bfgs.BFGS_ParameterList()

class soupy.optimization.bfgs.BFGS_operator(parameters=hippylib.ParameterList)

Bases: object

set_H0inv(H0inv)

Set user-defined operator corresponding to H0inv

Parameters:

H0inv – Fenics operator with method solve()

solve(x, b)

Solve system: \(H_{\mathrm{bfgs}} x = b\) where \(H_{\mathrm{bfgs}}\) is the approximation to the Hessian build by BFGS. That is, we apply \(x = (H_{\mathrm{bfgs}})^{-1} b = H_k b\) where \(H_k\) matrix is BFGS approximation to the inverse of the Hessian. Computation done via double-loop algorithm.

Parameters:

• x (dolfin.Vector) – The solution to the system

• b (dolfin.Vector) – The right-hand side of the system

update(s, y)

Update BFGS operator with most recent gradient update. To handle potential break from secant condition, update done via damping

Parameters:

• s (dolfin.Vector) – The update in medium parameters

• y (dolfin.Vector) – The update in gradient

soupy.optimization.bfgs.BFGSoperator_ParameterList()

class soupy.optimization.bfgs.RescaledIdentity(init_vector=None)

Bases: object

Default operator for H0inv, corresponds to applying \(d0 I\)

init_vector(x, dim)

solve(x, b)

Applies the operator \(d0 I\)

Parameters:

• x (dolfin.Vector) – solution vector

• b (dolfin.Vector) – right-hand side vector

soupy.optimization.inexactNewtonCG

class soupy.optimization.inexactNewtonCG.InexactNewtonCG(cost_functional, parameters=hippylib.ParameterList, preconditioner=None, norm_weighting=None, callback=None)

Bases: object

Inexact Newton-CG method to solve constrained optimization problems in the reduced parameter space. The Newton system is solved inexactly by early termination of CG iterations via Eisenstat-Walker (to prevent oversolving) and Steihaug (to avoid negative curvature) criteria. Globalization is performed using one of the following methods:

• line search (LS) based on the armijo sufficient reduction condition; or

• trust region (TR) based on the prior preconditioned norm of the update direction.

The stopping criterion is based on a control on the norm of the gradient and a control of the inner product between the gradient and the Newton direction.

The user must provide a model that describes the forward problem, cost functionals, and all the derivatives for the gradient and the Hessian.

More specifically the model object should implement following methods:

• cost(z) -> evaluate the cost functional

• grad(g) -> evaluate the gradient and store to g

• hessian(zhat, Hzhat) -> apply the cost Hessian

Type help(Model) for additional information

Constructor for the SGD solver

Parameters:

• cost_functional (soupy.ControlCostFunctional or similar) – The cost functional object

• parameters (hippylib.ParameterList.) – The parameters of the solver. Type InexactNewtonCG_ParameterList().showMe() for list of default parameters and their descriptions.

• preconditioner – Optional preconditioner for the Hessian inverse with method mult and solve Has optional method setLinearizationPoint

• norm_weighting – Weighting matrix for the norm with method mult

solve(z)

Solve the constrained optimization problem with initial guess z

Parameters:

z (dolfin.Vector) – The initial guess

Returns:

The optimization solution z and a dictionary of information

Note

The input z will be overwritten

termination_reasons = ['Maximum number of Iteration reached', 'Norm of the gradient less than tolerance', 'Maximum number of backtracking reached', 'Norm of (g, dz) less than tolerance']

soupy.optimization.inexactNewtonCG.InexactNewtonCG_ParameterList()

Generate a ParameterList for InexactNewtonCG. type InexactNewtonCG_ParameterList().showMe() for default values and their descriptions

soupy.optimization.inexactNewtonCG.LS_ParameterList()

Generate a ParameterList for line search globalization. type LS_ParameterList().showMe() for default values and their descriptions

soupy.optimization.inexactNewtonCG.TR_ParameterList()

Generate a ParameterList for Trust Region globalization. type RT_ParameterList().showMe() for default values and their descriptions

soupy.optimization.steepestDescent

class soupy.optimization.steepestDescent.SteepestDescent(cost_functional, parameters=hippylib.ParameterList)

Bases: object

Gradient Descent to solve optimization under uncertainty problems Globalization is performed using the Armijo sufficient reduction condition (backtracking). The stopping criterion is based on a control on the norm of the gradient.

The user must provide a cost functional that provides the evaluation and gradient

More specifically the cost functional object should implement following methods:

• generate_vector() -> generate the object containing state, parameter, adjoint, and control.

• cost(z) -> evaluate the cost functional

• grad(g) -> evaluate the gradient of the cost functional

Constructor for the Steepest Descent solver.

Parameters:

• cost_functional (soupy.ControlCostFunctional or similar) – The cost functional object

• parameters (hippylib.ParameterList.) – The parameters of the solver. Type SteepestDescent_ParameterList().showMe() for list of default parameters and their descriptions.

solve(z, box_bounds=None, constraint_projection=None)

Solve the constrained optimization problem using steepest descent with initial guess z

Parameters:

• z (dolfin.Vector) – The initial guess

• box_bounds (list) – Bound constraint. A list with two entries (min and max). Can be either a scalar value or a dolfin.Vector of the same size as z

• constraint_projection (ProjectableConstraint) – Alternative projectable constraint

Returns:

The optimization solution z and a dictionary of information

Note

The input z is overwritten

termination_reasons = ['Maximum number of Iteration reached', 'Norm of the gradient less than tolerance', 'Maximum number of backtracking reached', 'Norm of the step less than tolerance']

soupy.optimization.steepestDescent.SteepestDescent_ParameterList()

soupy.optimization.sgd

class soupy.optimization.sgd.SGD(cost_functional, parameters=hippylib.ParameterList)

Bases: object

Stochastic Gradient Descent to solve optimization under uncertainty problems The stopping criterion is based on a control on the norm of the gradient. Line search is only possible when the option stochastic_approximation is False.

The user must provide a cost functional that provides the evaluation and gradient

More specifically the cost functional object should implement following methods:

• generate_vector() -> generate the object containing state, parameter, adjoint, and control.

• cost(z, order, rng) -> evaluate the cost functional which allows a given rng

• grad(g) -> evaluate the gradient of the cost functional

Constructor for the SGD solver

Parameters:

• cost_functional (soupy.ControlCostFunctional or similar) – The cost functional object

• parameters (hippylib.ParameterList.) – The parameters of the solver. Type SGD_ParameterList().showMe() for list of default parameters and their descriptions.

solve(z, box_bounds=None, constraint_projection=None)

Solve the constrained optimization problem using steepest descent with initial guess z

Parameters:

• z (dolfin.Vector) – The initial guess

• box_bounds (list) – Bound constraint. A list with two entries (min and max). Can be either a scalar value or a dolfin.Vector of the same size as z

• constraint_projection (ProjectableConstraint) – Alternative projectable constraint

Returns:

The optimization solution z and summary of iterates

Note

The input z is overwritten

termination_reasons = ['Maximum number of Iteration reached', 'Norm of the gradient less than tolerance', 'Maximum number of backtracking reached', 'Norm of the step less than tolerance']

soupy.optimization.sgd.SGD_ParameterList()

soupy.optimization.projectableConstraint

class soupy.optimization.projectableConstraint.InnerProductEqualityConstraint(c, a)

Bases: ProjectableConstraint

Class implements the constraint \(c^T z - a = 0\)

Parameters:

• c (dolfin.Vector) – Constraint vector

• a (float) – Value of the constraint

cost(z)

Returns constraint violation math:c^T z - a for input z :param z: :type z: dolfin.Vector

Returns:

\(c^T z - a\) for input z

project(z)

Projects the vector z onto the hyperplane \(\{z : c^T z = a\}\)

Parameters:

z (dolfin.Vector)

class soupy.optimization.projectableConstraint.ProjectableConstraint

Bases: object

Base class for a constraint for which a projection operator into the feasible set is available

cost(z)

Returns the amount of violation of the constraint by z

Parameters:

z (dolfin.Vector)

project(z)

Projects the vector z onto the constraint set

Parameters:

z (dolfin.Vector)

soupy.optimization.cgSolverSteihaug

class soupy.optimization.cgSolverSteihaug.CGSolverSteihaug(parameters=hippylib.ParameterList)

Bases: object

Solve the linear system \(A x = b\) using preconditioned conjugate gradient ( \(B\) preconditioner) and the Steihaug stopping criterion:

• reason of termination 0: we reached the maximum number of iterations (no convergence)

• reason of termination 1: we reduced the residual up to the given tolerance (convergence)

• reason of termination 2: we reached a negative direction (premature termination due to not spd matrix)

• reason of termination 3: we reached the boundary of the trust region

The stopping criterion is based on either

• the absolute preconditioned residual norm check: \(|| r^* ||_{B^{-1}} < atol\)

• the relative preconditioned residual norm check: \(|| r^* ||_{B^{-1}}/|| r^0 ||_{B^{-1}} < rtol,\)

where \(r^* = b - Ax^*\) is the residual at convergence and \(r^0 = b - Ax^0\) is the initial residual.

The operator A is set using the method set_operator(A). A must provide the following two methods:

• A.mult(x,y): y = Ax

• A.init_vector(x, dim): initialize the vector x so that it is compatible with the range (dim = 0) or the domain (dim = 1) of A.

The preconditioner B is set using the method set_preconditioner(B). B must provide the following method: - B.solve(z,r): z is the action of the preconditioner B on the vector r

To solve the linear system \(Ax = b\) call self.solve(x,b). Here x and b are assumed to be dolfin.Vector objects.

Type CGSolverSteihaug_ParameterList().showMe() for default parameters and their descriptions

reason = ['Maximum Number of Iterations Reached', 'Relative/Absolute residual less than tol', 'Reached a negative direction', 'Reached trust region boundary']

set_TR(radius, B_op=None)

set_operator(A)

Set the operator \(A\).

set_preconditioner(B_solver)

Set the preconditioner \(B\).

solve(x, b)

Solve the linear system \(Ax = b\)

Parameters:

• x (dolfin.Vector) – Solution vector

• b (dolfin.Vector) – Right hand side vector

update_x_with_TR(x, alpha, d)

update_x_without_TR(x, alpha, d)

soupy.optimization.cgSolverSteihaug.CGSolverSteihaug_ParameterList()

Generate a :code:hippylib.`ParameterList` for CGSolverSteihaug. Type CGSolverSteihaug_ParameterList().showMe() for default values and their descriptions

class soupy.optimization.cgSolverSteihaug.IdentityOperator

Bases: object

mult(x, b)

solve(x, b)