From 887f1a07436707d847101fcb93fc6d0d15108850 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Thu, 29 Jun 2017 15:00:20 -0400 Subject: [PATCH 01/18] added initial version of PrimalDualVector --- src/kona/linalg/vectors/composite.py | 42 ++++++- src/kona/test/test_primal_dual_vector.py | 143 +++++++++++++++++++++++ 2 files changed, 183 insertions(+), 2 deletions(-) create mode 100644 src/kona/test/test_primal_dual_vector.py diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index 1a4c04f..21f1baa 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -207,6 +207,44 @@ def infty(self): norms.append(self._vectors[i].infty) return max(norms) +class PrimalDualVector(CompositeVector): + """ + A composite vector made up of primal, dual equality, and dual inequality vectors. + + Parameters + ---------- + _memory : KonaMemory + All-knowing Kona memory manager. + _primal : DesignVector + Primal component of the composite vector. + _dual : DualVector + Dual components of the composite vector. + """ + + init_dual = 0.0 # default initial value for multipliers + + def __init__(self, primal_vec, dual_vec): + assert isinstance(primal_vec, DesignVector), \ + 'PrimalDualVector() >> Mismatched primal vector. ' + \ + 'Must be DesignVector!' + assert isinstance(dual_vec, DualVectorEQ) or \ + isinstance(dual_vec, DualVectorINEQ) or \ + isinstance(dual_vec, CompositeDualVector), \ + 'PrimalDualVector() >> Mismatched dual vector. ' + \ + 'Must be DualVectorEQ, DualVectorINEQ CompositeDualVector!' + + self.primal = primal_vec + self.dual = dual_vec + + super(PrimalDualVector, self).__init__([primal_vec, dual_vec]) + + def equals_init_guess(self): + """ + Sets the primal-dual vector to the initial guess, using the initial design. + """ + self.primal.equals_init_design() + self.dual.equals(self.init_dual) + class ReducedKKTVector(CompositeVector): """ A composite vector representing a combined primal and dual vectors. @@ -369,7 +407,7 @@ def equals_constraints(self, at_primal, at_state, scale=1.0): class CompositePrimalVector(CompositeVector): """ - A composite vector representing a combined design and slack vectors.. + A composite vector representing a combined design and slack vectors. Parameters ---------- @@ -466,4 +504,4 @@ def equals_lagrangian_total_gradient( # package imports at the bottom to prevent import errors import numpy as np from kona.linalg.vectors.common import DesignVector -from kona.linalg.vectors.common import DualVectorEQ, DualVectorINEQ \ No newline at end of file +from kona.linalg.vectors.common import DualVectorEQ, DualVectorINEQ diff --git a/src/kona/test/test_primal_dual_vector.py b/src/kona/test/test_primal_dual_vector.py new file mode 100644 index 0000000..abb8a83 --- /dev/null +++ b/src/kona/test/test_primal_dual_vector.py @@ -0,0 +1,143 @@ +import unittest +import numpy as np + +from kona.linalg.memory import KonaMemory +from dummy_solver import DummySolver +from kona.linalg.matrices.common import * +from kona.linalg.vectors.composite import CompositePrimalVector +from kona.linalg.vectors.composite import CompositeDualVector +from kona.linalg.vectors.composite import PrimalDualVector + +class PrimalDualVectorTestCase(unittest.TestCase): + + def setUp(self): + solver = DummySolver(10, 5, 5) # 10 DVs, 5 equality, 5 inequality + self.km = km = KonaMemory(solver) + + km.primal_factory.request_num_vectors(5) + km.state_factory.request_num_vectors(3) + km.eq_factory.request_num_vectors(4) + km.ineq_factory.request_num_vectors(8) # doe we need this many? + km.allocate_memory() + + # get some work vectors + self.primal_work = km.primal_factory.generate() + self.slack_work = km.ineq_factory.generate() # only used for error testing + self.state_work = km.state_factory.generate() + self.eq_work = km.eq_factory.generate() + self.ineq_work = km.ineq_factory.generate() + + # get static eval point vectors + self.design = km.primal_factory.generate() + self.state = km.state_factory.generate() + self.adjoint = km.state_factory.generate() + self.dual_eq = km.eq_factory.generate() + self.dual_ineq = km.ineq_factory.generate() + self.dual = CompositeDualVector(self.dual_eq, self.dual_ineq) + self.at_kkt1 = PrimalDualVector(self.design, self.dual) + self.at_kkt2 = PrimalDualVector(self.design, self.dual_ineq) + self.at_kkt3 = PrimalDualVector(self.design, self.dual_eq) + + # set the evaluation point + self.design.equals_init_design() + self.state.equals_primal_solution(self.design) + self.dual.equals(1.0) + + # case 1: primal-dual vector is a complete vector with design, eq and ineq + self.pv1 = km.primal_factory.generate() + self.eq1 = km.eq_factory.generate() + self.ineq1 = km.ineq_factory.generate() + self.pv1.base.data = 2*np.ones(10) + self.eq1.base.data = 4*np.ones(5) + self.ineq1.base.data = 5*np.ones(5) + primal1 = self.pv1 + dual1 = CompositeDualVector(self.eq1, self.ineq1) + self.pd_vec1 = PrimalDualVector(primal1, dual1) + + # case 2: primal-dual vector has design and ineq + self.pv2 = km.primal_factory.generate() + self.ineq2 = km.ineq_factory.generate() + self.pv2.base.data = 2*np.ones(10) + self.ineq2.base.data = 4*np.ones(5) + primal2 = self.pv2 + dual2 = self.ineq2 + self.pd_vec2 = PrimalDualVector(primal2, dual2) + + # case 3: primal-dual vector only has design and eq + self.pv3 = km.primal_factory.generate() + self.eq3 = km.eq_factory.generate() + primal3 = self.pv3 + dual3 = self.eq3 + self.pd_vec3 = PrimalDualVector(primal3, dual3) + + def test_bad_init_args(self): + try: + PrimalDualVector(self.eq1, self.eq1) + except AssertionError as err: + self.assertEqual( + str(err), + "PrimalDualVector() >> Mismatched primal vector. " + + "Must be DesignVector!") + else: + self.fail('AssertionError expected') + + try: + PrimalDualVector(self.pv1, self.pv1) + except AssertionError as err: + self.assertEqual( + str(err), + "PrimalDualVector() >> Mismatched dual vector. " + + "Must be DualVectorEQ, DualVectorINEQ CompositeDualVector!") + else: + self.fail('AssertionError expected') + + try: + PrimalDualVector(self.ineq1, self.ineq1) + except AssertionError as err: + self.assertEqual( + str(err), + "PrimalDualVector() >> Mismatched primal vector. " + + "Must be DesignVector!") + else: + self.fail('AssertionError expected') + + try: + PrimalDualVector( + CompositePrimalVector(self.pv1, self.slack_work), self.eq1) + except AssertionError as err: + self.assertEqual( + str(err), + "PrimalDualVector() >> Mismatched primal vector. " + + "Must be DesignVector!") + else: + self.fail('AssertionError expected') + + def test_initial_guess_case1(self): + self.pd_vec1.equals_init_guess() + + err = self.pv1.base.data - 10 * np.ones(10) + self.assertEqual(np.linalg.norm(err), 0) + + err = self.eq1.base.data - self.pd_vec1.init_dual * (np.ones(5)) + self.assertEqual(np.linalg.norm(err), 0) + + err = self.ineq1.base.data - self.pd_vec1.init_dual * (np.ones(5)) + self.assertEqual(np.linalg.norm(err), 0) + + def test_init_guess_case2(self): + self.pd_vec2.equals_init_guess() + + err = self.pv2.base.data - 10 * np.ones(10) + self.assertEqual(np.linalg.norm(err), 0) + + err = self.ineq2.base.data - self.pd_vec2.init_dual * (np.ones(5)) + self.assertEqual(np.linalg.norm(err), 0) + + def test_init_guess_case3(self): + self.pd_vec3.equals_init_guess() + + err = self.pv3.base.data - 10*np.ones(10) + self.assertEqual(np.linalg.norm(err), 0) + + err = self.eq3.base.data - self.pd_vec3.init_dual*(np.ones(5)) + self.assertEqual(np.linalg.norm(err), 0) From 1e04ecd78606e8512eb3580d4b3dae2c4c0a5062 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Fri, 30 Jun 2017 10:37:52 -0400 Subject: [PATCH 02/18] added equals_opt_residual for PrimalDualVector type --- src/kona/linalg/vectors/common.py | 24 +++++++- src/kona/linalg/vectors/composite.py | 78 ++++++++++++++++++++++++ src/kona/test/test_dual_vectors.py | 19 +++++- src/kona/test/test_primal_dual_vector.py | 26 +++++++- 4 files changed, 142 insertions(+), 5 deletions(-) diff --git a/src/kona/linalg/vectors/common.py b/src/kona/linalg/vectors/common.py index 95fd618..7e26c66 100644 --- a/src/kona/linalg/vectors/common.py +++ b/src/kona/linalg/vectors/common.py @@ -682,9 +682,31 @@ def equals_constraints(self, at_primal, at_state, scale=1.0): at_design.base.data, at_state.base) self.times(scale) + def equals_mangasarian(self, ineq, mult): + """ + Evaluate the nonlinear function, due to Mangasarian, that is equivalent to the + complimentarity conditions, multiplier bounds and inequality constraints: + :math:`G(g,\\lambda_g) = -|g - \\lambda_g|^3 + (g)^3 + (\\lambda_g)^3`. + Stores the result in-place. + + Parameters + ---------- + ineq : DualVectorINEQ + Inequality constraint at the current design and state. + mult : DualVectorINEQ + The multiplier corresponding to the inequality constraints. + """ + assert isinstance(ineq, DualVectorINEQ), \ + 'DualVectorINEQ >> Invalid type for ineq. ' + \ + 'Must be DualVecINEQ!' + assert isinstance(mult, DualVectorINEQ), \ + 'DualVectorINEQ >> Invalid type for mult. ' + \ + 'Must be DualVecINEQ!' + self.base.data[:] = -np.fabs(ineq.base.data[:] - mult.base.data[:])**3 \ + + ineq.base.data[:]**3 + mult.base.data[:]**3 # package imports at the bottom to prevent import errors import numpy as np from kona.linalg.vectors.composite import ReducedKKTVector from kona.linalg.vectors.composite import CompositePrimalVector -from kona.linalg.matrices.common import * \ No newline at end of file +from kona.linalg.matrices.common import * diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index 21f1baa..3a6b66e 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -245,6 +245,84 @@ def equals_init_guess(self): self.primal.equals_init_design() self.dual.equals(self.init_dual) + def equals_opt_residual(self, x, state, adjoint, barrier=None, + obj_scale=1.0, cnstr_scale=1.0): + """ + Calculates the following nonlinear vector function: + + .. math:: + r(x,\\lambda_h,\\lambda_g) = + \\begin{bmatrix} + \\nabla_x f(x, u) - \\nabla_x h(x, u)^T \\lambda_{h} - \\nabla_x g(x, u)^T \\lambda_{g} \\\\ + h(x,u) \\\\ + -\\theta(|g(x,u) - \\lambda_g|) + \\theta(g(x,u)) + \theta(\\lambda_g) + \\end{bmatrix} + + where :math:`h(x,u)` are the equality constraints, and :math:`g(x,u)` are the + inequality constraints. The vectors :math:`\\lambda_h` and :math:`\\lambda_g` + are the associated Lagrange multipliers. The function :math:`\\theta(z)` + is any strictly increasing function; typically :math:`\\theta(z) = z^3`. The + solution to :math:`r(x,\\lambda_h,\\lambda_g) = 0` is equivalent to the + first-order optimality conditions for the optimization problem + + .. math:: + \\begin{align*} + \\min_x &f(x,u(x)) \\\\ + \\textsf{s.t.} &h(x,u(x)) = 0, \\\\ + &g(x,u(x)) \geq 0. + \\end{align*} + + Parameters + ---------- + x : PrimalDualVector + Evaluate first-order optimality conditions at this primal-dual point. + state : StateVector + Evaluate first-order optimality conditions at this state point. + adjoint : StateVector + Evaluate first-order optimality conditions using this adjoint vector. + obj_scale : float, optional + Scaling for the objective function. + cnstr_scale : float, optional + Scaling for the constraints. + """ + assert isinstance(x, PrimalDualVector), \ + "PrimalDualVector() >> invalid type x in equals_opt_residual. " + \ + "x vector must be a PrimalDualVector!" + if isinstance(x.dual, CompositeDualVector): + dual_eq = x.dual.eq + dual_ineq = x.dual.ineq + elif isinstance(x.dual, DualVectorINEQ): + dual_eq = None + dual_ineq = x.dual + elif isinstance(x.dual, DualVectorEQ): + dual_eq = x.dual + dual_ineq = None + else: + raise AssertionError("PrimalDualVector() >> invalid dual vector type") + design = x.primal + + design_opt = self.primal + + # first include the objective partial and adjoint contribution + design_opt.equals_total_gradient(design, state, adjoint, obj_scale) + if dual_eq is not None: + # add the Lagrange multiplier products for equality constraints + design_opt.base.data[:] += design_opt._memory.solver.multiply_dCEQdX_T( + design.base.data, state.base, dual_eq.base.data) * \ + cnstr_scale + if dual_ineq is not None: + # add the Lagrange multiplier products for inequality constraints + design_opt.base.data[:] += design_opt._memory.solver.multiply_dCINdX_T( + design.base.data, state.base, dual_ineq.base.data) * \ + cnstr_scale + # include constraint terms + self.dual.equals_constraints(design, state, cnstr_scale) + if isinstance(self.dual, DualVectorINEQ): + self.dual.equals_mangasarian(self.dual, dual_ineq) + elif isinstance(self.dual, CompositeDualVector): + self.dual.ineq.equals_mangasarian(self.dual.ineq, dual_ineq) + print self.dual.ineq.base.data + class ReducedKKTVector(CompositeVector): """ A composite vector representing a combined primal and dual vectors. diff --git a/src/kona/test/test_dual_vectors.py b/src/kona/test/test_dual_vectors.py index 47370fd..65ad066 100644 --- a/src/kona/test/test_dual_vectors.py +++ b/src/kona/test/test_dual_vectors.py @@ -35,12 +35,14 @@ def setUp(self): km.primal_factory.request_num_vectors(1) km.state_factory.request_num_vectors(1) - km.ineq_factory.request_num_vectors(1) + km.ineq_factory.request_num_vectors(3) km.allocate_memory() self.pv = km.primal_factory.generate() self.sv = km.state_factory.generate() self.dv = km.ineq_factory.generate() + self.mult = km.ineq_factory.generate() + self.mang = km.ineq_factory.generate() def test_equals_constraints(self): at_design = self.pv @@ -50,6 +52,21 @@ def test_equals_constraints(self): self.dv.equals_constraints(at_design, at_state) self.assertEqual(self.dv.inner(self.dv), 1000) + def test_equals_mangasarian(self): + at_design = self.pv + at_design.equals(1) + at_state = self.sv + at_state.equals(2) + self.dv.equals_constraints(at_design, at_state) + self.mult.equals(-9) + self.mang.equals_mangasarian(self.dv, self.mult) + for i in xrange(len(self.mang.base.data)): + self.assertEqual(self.mang.base.data[i], -1730) + self.mult.equals(-11) + self.mang.equals_mangasarian(self.dv, self.mult) + for i in xrange(len(self.mang.base.data)): + self.assertEqual(self.mang.base.data[i], -2332) + class DualVectorEQTestCaseIDF(unittest.TestCase): def setUp(self): diff --git a/src/kona/test/test_primal_dual_vector.py b/src/kona/test/test_primal_dual_vector.py index abb8a83..010c2dc 100644 --- a/src/kona/test/test_primal_dual_vector.py +++ b/src/kona/test/test_primal_dual_vector.py @@ -34,9 +34,9 @@ def setUp(self): self.dual_eq = km.eq_factory.generate() self.dual_ineq = km.ineq_factory.generate() self.dual = CompositeDualVector(self.dual_eq, self.dual_ineq) - self.at_kkt1 = PrimalDualVector(self.design, self.dual) - self.at_kkt2 = PrimalDualVector(self.design, self.dual_ineq) - self.at_kkt3 = PrimalDualVector(self.design, self.dual_eq) + self.at_pd1 = PrimalDualVector(self.design, self.dual) + self.at_pd2 = PrimalDualVector(self.design, self.dual_ineq) + self.at_pd3 = PrimalDualVector(self.design, self.dual_eq) # set the evaluation point self.design.equals_init_design() @@ -141,3 +141,23 @@ def test_init_guess_case3(self): err = self.eq3.base.data - self.pd_vec3.init_dual*(np.ones(5)) self.assertEqual(np.linalg.norm(err), 0) + + def test_opt_residual_case1(self): + # case 1 has both equality and inequality constraints + # recall: self.dual = 1, so dCeqdU^T*dual.eq + dCineqdU^Tdual.ineq = 5 + 5 = 10 + dCdU(self.design, self.state).T.product(self.dual, self.adjoint, self.state_work) + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.plus(self.adjoint) + self.state_work.times(-1.) + # We get adjoint = 9 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + self.pd_vec1.equals_opt_residual(self.at_pd1, self.state, self.adjoint) + + # check results + exp_dLdX_norm = np.sqrt(10. * 20.**2) + self.assertEqual(self.pd_vec1.primal.norm2, exp_dLdX_norm) + exp_dLdEq_norm = np.sqrt(5. * 200.**2) + self.assertEqual(self.pd_vec1.dual.eq.norm2, exp_dLdEq_norm) + exp_dLdIn_norm = np.sqrt(5. * 119402**2) + self.assertEqual(self.pd_vec1.dual.ineq.norm2, exp_dLdIn_norm) From 1725bc8ba0af6ae37e3cfb3d6a9848887acc0c7f Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Fri, 30 Jun 2017 11:23:24 -0400 Subject: [PATCH 03/18] added MultiSecantApprox abstract base class --- src/kona/linalg/matrices/hessian/basic.py | 58 ++++++++++++++++++++++- 1 file changed, 56 insertions(+), 2 deletions(-) diff --git a/src/kona/linalg/matrices/hessian/basic.py b/src/kona/linalg/matrices/hessian/basic.py index 6427eb6..60ed2e3 100644 --- a/src/kona/linalg/matrices/hessian/basic.py +++ b/src/kona/linalg/matrices/hessian/basic.py @@ -100,7 +100,7 @@ class QuasiNewtonApprox(BaseHessian): def __init__(self, vector_factory, optns={}): assert isinstance(vector_factory, VectorFactory), \ - "LimitedMemoryBFGS() >> Invalid vector factory!" + "QuasiNewtonApprox() >> Invalid vector factory!" super(QuasiNewtonApprox, self).__init__(vector_factory, optns) self.max_stored = get_opt(optns, 10, 'max_stored') @@ -121,9 +121,63 @@ def add_correction(self, s_new, y_new): """ raise NotImplementedError # pragma: no cover +class MultiSecantApprox(BaseHessian): + """ Base class for multi-secant approximations of Jacobians + + Attributes + ---------- + max_stored : int + Maximum number of previous steps and residuals stored. + ptr : list of integers + Integer "pointers" that indicate the oldest to newest data + x_hist : list of KonaVector + The sequence of solutions used to build x_diff + r_hist : list of KonaVector + The sequence of residuals used to build r_diff + x_diff : list of KonaVector + Difference between subsequent solutions: :math:`\\Delta x_k = x_{k+1} - x_k` + r_diff : list of KonaVector + Difference between subsequent residuals: :math:`\\Delta r_k = r_{k+1} - r_k` + """ + + def __init__(self, vector_factory, optns={}): + assert isinstance(vector_factory, VectorFactory), \ + "MultiSecantApprox() >> Invalid vector factory!" + super(MultiSecantApprox, self).__init__(vector_factory, optns) + self.max_stored = get_opt(optns, 10, 'max_stored') + self.ptr = [] + self.x_hist = [] + self.r_hist = [] + self.x_diff = [] + self.r_diff = [] + + def add_to_history(self, x_new, r_new): + """ + Adds to the solution and residual history + + Parameters + ---------- + x_new : KonaVector + New solution vector. + r_new : KonaVector + New residual vector, or data that can be used to build residual. + """ + raise NotImplementedError # pragma: no cover + + def build_difference_matrices(mu=0.0): + """ + Constructs the solution and residual differences from the history + + Parameters + ---------- + mu : float + Homotopy continuation parameter + """ + raise NotImplementedError # pragma: no cover + # imports at the bottom to prevent circular import errors import sys from kona.options import get_opt from kona.linalg.memory import VectorFactory from kona.linalg.vectors.common import DesignVector, StateVector -from kona.linalg.vectors.common import DualVectorEQ, DualVectorINEQ \ No newline at end of file +from kona.linalg.vectors.common import DualVectorEQ, DualVectorINEQ From fd3ff4f41da1bd15a9ac9c3ed6896d2fc18f7846 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Mon, 3 Jul 2017 12:38:35 -0400 Subject: [PATCH 04/18] refactored implementation of PrimalDualVector --- src/kona/linalg/vectors/common.py | 24 ++- src/kona/linalg/vectors/composite.py | 153 +++++++++++++-- src/kona/test/test_primal_dual_vector.py | 239 ++++++++++++++++++++--- 3 files changed, 360 insertions(+), 56 deletions(-) diff --git a/src/kona/linalg/vectors/common.py b/src/kona/linalg/vectors/common.py index 7e26c66..3e330d3 100644 --- a/src/kona/linalg/vectors/common.py +++ b/src/kona/linalg/vectors/common.py @@ -59,11 +59,15 @@ def plus(self, vector): Parameters ---------- - vector : KonaVector + vector : float or KonaVector Vector to be added. """ - assert isinstance(vector, type(self)) - self.base.plus(vector.base) + if isinstance(vector, + (float, np.float32, np.float64, int, np.int32, np.int64)): + self.base.data[:] += vector + else: + assert isinstance(vector, type(self)) + self.base.plus(vector.base) def minus(self, vector): """ @@ -73,16 +77,20 @@ def minus(self, vector): Parameters ---------- - vector : KonaVector + vector : float or KonaVector Vector to be subtracted. """ if vector == self: # special case... self.equals(0) - assert isinstance(vector, type(self)) - self.base.times_scalar(-1.) - self.base.plus(vector.base) - self.base.times_scalar(-1.) + if isinstance(vector, + (float, np.float32, np.float64, int, np.int32, np.int64)): + self.base.data[:] -= vector + else: + assert isinstance(vector, type(self)) + self.base.times_scalar(-1.) + self.base.plus(vector.base) + self.base.times_scalar(-1.) def times(self, factor): """ diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index 3a6b66e..d9c075a 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -18,6 +18,7 @@ def _check_type(self, vec): except TypeError: raise TypeError("CompositeVector() >> " + "Mismatched internal vectors!") + def equals(self, rhs): """ Used as the assignment operator. @@ -215,38 +216,153 @@ class PrimalDualVector(CompositeVector): ---------- _memory : KonaMemory All-knowing Kona memory manager. - _primal : DesignVector + primal : DesignVector Primal component of the composite vector. - _dual : DualVector - Dual components of the composite vector. + eq : DualVectorEQ + Dual component corresponding to the equality constraints. + ineq: DualVectorINEQ + Dual component corresponding to the inequality constraints. """ init_dual = 0.0 # default initial value for multipliers - def __init__(self, primal_vec, dual_vec): + def __init__(self, primal_vec, eq_vec=None, ineq_vec=None): assert isinstance(primal_vec, DesignVector), \ - 'PrimalDualVector() >> Mismatched primal vector. ' + \ - 'Must be DesignVector!' - assert isinstance(dual_vec, DualVectorEQ) or \ - isinstance(dual_vec, DualVectorINEQ) or \ - isinstance(dual_vec, CompositeDualVector), \ - 'PrimalDualVector() >> Mismatched dual vector. ' + \ - 'Must be DualVectorEQ, DualVectorINEQ CompositeDualVector!' - + 'PrimalDualVector() >> primal_vec must be a DesignVector!' + if eq_vec is not None: + assert isinstance(eq_vec, DualVectorEQ), \ + 'PrimalDualVector() >> eq_vec must be a DualVectorEQ!' + if ineq_vec is not None: + assert isinstance(ineq_vec, DualVectorINEQ), \ + 'PrimalDualVector() >> ineq_vec must be a DualVectorINEQ!' self.primal = primal_vec - self.dual = dual_vec - - super(PrimalDualVector, self).__init__([primal_vec, dual_vec]) + self.eq = eq_vec + self.ineq = ineq_vec + super(PrimalDualVector, self).__init__([primal_vec, eq_vec, ineq_vec]) def equals_init_guess(self): """ Sets the primal-dual vector to the initial guess, using the initial design. """ self.primal.equals_init_design() - self.dual.equals(self.init_dual) + if self.eq is not None: + self.eq.equals(self.init_dual) + if self.ineq is not None: + self.ineq.equals(self.init_dual) + + def equals_KKT_conditions(self, x, state, adjoint, obj_scale=1.0, cnstr_scale=1.0): + """ + Calculates the total derivative of the Lagrangian + :math:`\\mathcal{L}(x, u) = f(x, u)+ \\lambda_{h}^T h(x, u) + \\lambda_{g}^T g(x, u)` with + respect to :math:`\\begin{pmatrix}x && \\lambda_{h} && \\lambda_{g}\\end{pmatrix}^T`, + where :math:`h` denotes the equality constraints (if any) and :math:`g` denotes + the inequality constraints (if any). Note that these (total) derivatives do not + represent the complete set of first-order optimality conditions in the case of + inequality constraints. - def equals_opt_residual(self, x, state, adjoint, barrier=None, - obj_scale=1.0, cnstr_scale=1.0): + Parameters + ---------- + x : PrimalDualVector + Evaluate derivatives at this primal-dual point. + state : StateVector + Evaluate derivatives at this state point. + adjoint : StateVector + Evaluate derivatives using this adjoint vector. + obj_scale : float, optional + Scaling for the objective function. + cnstr_scale : float, optional + Scaling for the constraints. + """ + assert isinstance(x, PrimalDualVector), \ + "PrimalDualVector() >> invalid type x in equals_opt_residual. " + \ + "x vector must be a PrimalDualVector!" + if x.eq is None or self.eq is None: + assert self.eq is None and x.eq is None, \ + "PrimalDualVector() >> inconsistency with x.eq and self.eq!" + if x.ineq is None or self.ineq is None: + assert self.ineq is None and x.ineq is None, \ + "PrimalDualVector() >> inconsistency with x.ineq and self.ineq!" + # set some aliases + design = x.primal + dLdx = self.primal + ceq = self.eq + cineq = self.ineq + # first include the objective partial and adjoint contribution + dLdx.equals_total_gradient(design, state, adjoint, obj_scale) + if x.eq is not None: + # add the Lagrange multiplier products for equality constraints + dLdx.base.data[:] += dLdx._memory.solver.multiply_dCEQdX_T( + design.base.data, state.base, x.eq.base.data) * \ + cnstr_scale + if x.ineq is not None: + # add the Lagrange multiplier products for inequality constraints + dLdx.base.data[:] += dLdx._memory.solver.multiply_dCINdX_T( + design.base.data, state.base, x.ineq.base.data) * \ + cnstr_scale + # include constraint terms + if ceq is not None: + ceq.equals_constraints(design, state, cnstr_scale) + if cineq is not None: + cineq.equals_constraints(design, state ,cnstr_scale) + + def equals_homotopy_residual(self, dLdx, x, init, mu=1.0): + """ + Using dLdx=:math:`\\begin{pmatrix} \\nabla_x L && h && g \\end{pmatrix}`, which can be + obtained from the method equals_KKT_conditions, as well as the initial values + init=:math:`\\begin{pmatrix} x_0 && h(x_0,u_0) && g(x_0,u_0) \\end{pmatrix}` and the + current point :math:`\\begin{pmatrix} x && \lambda_h && \lambda_g \end{pmatrix}`, this + method computes the following nonlinear vector function: + + .. math:: + r(x,\\lambda_h,\\lambda_g;\\mu) = + \\begin{bmatrix} + \\mu\\left[\\nabla_x f(x, u) - \\nabla_x h(x, u)^T \\lambda_{h} - \\nabla_x g(x, u)^T \\lambda_{g}\\right] + (1 - \\mu)(x - x_0) \\\\ + h(x,u) - (1-\mu)h(x_0,u_0) \\\\ + -|g(x,u) - (1-\\mu)*g_0 - \\lambda_g|^3 + (g(x,u) - (1-\\mu)g_0)^3 + \\lambda_g^3 - (1-\\mu)\\hat{g} + \\end{bmatrix} + + where :math:`h(x,u)` are the equality constraints, and :math:`g(x,u)` are the + inequality constraints. The vectors :math:`\\lambda_h` and :math:`\\lambda_g` + are the associated Lagrange multipliers. When mu=1.0, we recover a set of nonlinear + algebraic equations equivalent to the first-order optimality conditions. + + Parameters + ---------- + dLdx : PrimalDualVector + The total derivative of the Lagranginan with respect to the primal and dual variables. + x : PrimalDualVector + The current solution vector value corresponding to dLdx. + init: PrimalDualVector + The initial primal variable, as well as the initial constraint values. + mu : float + Homotopy parameter; must be between 0 and 1. + """ + assert isinstance(dLdx, PrimalDualVector), \ + "PrimalDualVector() >> dLdx must be a PrimalDualVector!" + assert isinstance(init, PrimalDualVector), \ + "PrimalDualVector() >> init must be a PrimalDualVector!" + if dLdx.eq is None or self.eq is None or x.eq is None or init.eq is None: + assert dLdx.eq is None and self.eq is None and x.eq is None and init.eq is None, \ + "PrimalDualVector() >> inconsistent eq component in self, dLdx, x, and/or init!" + if dLdx.ineq is None or self.ineq is None or x.ineq is None or init.ineq is None: + assert dLdx.ineq is None and self.ineq is None and x.ineq is None \ + and init.ineq is None, \ + "PrimalDualVector() >> inconsistent ineq component in self, dLdx, x, and/or init!" + if dLdx == self or x == self or init == self: + "PrimalDualVector() >> equals_homotopy_residual is not in-place!" + # construct the primal part of the residual: mu*dLdx + (1-mu)(x - x_0) + self.primal.equals_ax_p_by(1.0, x.primal, -1.0, init.primal) + self.primal.equals_ax_p_by(mu, dLdx.primal, (1.0 - mu), self.primal) + if dLdx.eq is not None: + # include the equality constraint part: h - (1-mu)h_0 + self.eq.equals_ax_p_by(1.0, dLdx.eq, -(1.0 - mu), init.eq) + if dLdx.ineq is not None: + # include the inequality constraint part + self.ineq.equals_ax_p_by(1.0, dLdx.ineq, -(1.0 - mu), init.ineq) + self.ineq.equals_mangasarian(self.ineq, x.ineq) + self.ineq.minus((1.0 - mu)*0.1) + + def equals_opt_residual(self, x, state, adjoint, obj_scale=1.0, cnstr_scale=1.0): """ Calculates the following nonlinear vector function: @@ -321,7 +437,6 @@ def equals_opt_residual(self, x, state, adjoint, barrier=None, self.dual.equals_mangasarian(self.dual, dual_ineq) elif isinstance(self.dual, CompositeDualVector): self.dual.ineq.equals_mangasarian(self.dual.ineq, dual_ineq) - print self.dual.ineq.base.data class ReducedKKTVector(CompositeVector): """ diff --git a/src/kona/test/test_primal_dual_vector.py b/src/kona/test/test_primal_dual_vector.py index 010c2dc..466f95b 100644 --- a/src/kona/test/test_primal_dual_vector.py +++ b/src/kona/test/test_primal_dual_vector.py @@ -14,18 +14,21 @@ def setUp(self): solver = DummySolver(10, 5, 5) # 10 DVs, 5 equality, 5 inequality self.km = km = KonaMemory(solver) - km.primal_factory.request_num_vectors(5) + km.primal_factory.request_num_vectors(10) km.state_factory.request_num_vectors(3) - km.eq_factory.request_num_vectors(4) - km.ineq_factory.request_num_vectors(8) # doe we need this many? + km.eq_factory.request_num_vectors(6) + km.ineq_factory.request_num_vectors(10) # doe we need this many? km.allocate_memory() # get some work vectors - self.primal_work = km.primal_factory.generate() + self.primal_work1 = km.primal_factory.generate() + self.primal_work2 = km.primal_factory.generate() self.slack_work = km.ineq_factory.generate() # only used for error testing self.state_work = km.state_factory.generate() - self.eq_work = km.eq_factory.generate() - self.ineq_work = km.ineq_factory.generate() + self.eq_work1 = km.eq_factory.generate() + self.eq_work2 = km.eq_factory.generate() + self.ineq_work1 = km.ineq_factory.generate() + self.ineq_work2 = km.ineq_factory.generate() # get static eval point vectors self.design = km.primal_factory.generate() @@ -33,15 +36,17 @@ def setUp(self): self.adjoint = km.state_factory.generate() self.dual_eq = km.eq_factory.generate() self.dual_ineq = km.ineq_factory.generate() - self.dual = CompositeDualVector(self.dual_eq, self.dual_ineq) - self.at_pd1 = PrimalDualVector(self.design, self.dual) - self.at_pd2 = PrimalDualVector(self.design, self.dual_ineq) - self.at_pd3 = PrimalDualVector(self.design, self.dual_eq) + self.at_pd1 = PrimalDualVector(self.design, eq_vec=self.dual_eq, + ineq_vec=self.dual_ineq) + self.at_pd2 = PrimalDualVector(self.design, ineq_vec=self.dual_ineq) + self.at_pd3 = PrimalDualVector(self.design, eq_vec=self.dual_eq) + self.at_pd4 = PrimalDualVector(self.design) # set the evaluation point self.design.equals_init_design() self.state.equals_primal_solution(self.design) - self.dual.equals(1.0) + self.dual_eq.equals(1.0) + self.dual_ineq.equals(1.0) # case 1: primal-dual vector is a complete vector with design, eq and ineq self.pv1 = km.primal_factory.generate() @@ -51,8 +56,13 @@ def setUp(self): self.eq1.base.data = 4*np.ones(5) self.ineq1.base.data = 5*np.ones(5) primal1 = self.pv1 - dual1 = CompositeDualVector(self.eq1, self.ineq1) - self.pd_vec1 = PrimalDualVector(primal1, dual1) + dual_eq1 = self.eq1 + dual_ineq1 = self.ineq1 + self.pd_vec1 = PrimalDualVector(primal1, eq_vec=dual_eq1, ineq_vec=dual_ineq1) + self.pd_work11 = PrimalDualVector(self.primal_work1, eq_vec=self.eq_work1, + ineq_vec=self.ineq_work1) + self.pd_work12 = PrimalDualVector(self.primal_work2, eq_vec=self.eq_work2, + ineq_vec=self.ineq_work2) # case 2: primal-dual vector has design and ineq self.pv2 = km.primal_factory.generate() @@ -61,14 +71,25 @@ def setUp(self): self.ineq2.base.data = 4*np.ones(5) primal2 = self.pv2 dual2 = self.ineq2 - self.pd_vec2 = PrimalDualVector(primal2, dual2) + self.pd_vec2 = PrimalDualVector(primal2, ineq_vec=dual2) + self.pd_work21 = PrimalDualVector(self.primal_work1, ineq_vec=self.ineq_work1) + self.pd_work22 = PrimalDualVector(self.primal_work2, ineq_vec=self.ineq_work2) # case 3: primal-dual vector only has design and eq self.pv3 = km.primal_factory.generate() self.eq3 = km.eq_factory.generate() primal3 = self.pv3 dual3 = self.eq3 - self.pd_vec3 = PrimalDualVector(primal3, dual3) + self.pd_vec3 = PrimalDualVector(primal3, eq_vec=dual3) + self.pd_work31 = PrimalDualVector(self.primal_work1, eq_vec=self.eq_work1) + self.pd_work32 = PrimalDualVector(self.primal_work2, eq_vec=self.eq_work2) + + # case 4: primal-dual vector has design only + self.pv4 = km.primal_factory.generate() + primal4 = self.pv4 + self.pd_vec4 = PrimalDualVector(primal4) + self.pd_work41 = PrimalDualVector(self.primal_work1) + self.pd_work42 = PrimalDualVector(self.primal_work2) def test_bad_init_args(self): try: @@ -76,8 +97,7 @@ def test_bad_init_args(self): except AssertionError as err: self.assertEqual( str(err), - "PrimalDualVector() >> Mismatched primal vector. " + - "Must be DesignVector!") + "PrimalDualVector() >> primal_vec must be a DesignVector!") else: self.fail('AssertionError expected') @@ -86,8 +106,7 @@ def test_bad_init_args(self): except AssertionError as err: self.assertEqual( str(err), - "PrimalDualVector() >> Mismatched dual vector. " + - "Must be DualVectorEQ, DualVectorINEQ CompositeDualVector!") + "PrimalDualVector() >> eq_vec must be a DualVectorEQ!") else: self.fail('AssertionError expected') @@ -96,8 +115,16 @@ def test_bad_init_args(self): except AssertionError as err: self.assertEqual( str(err), - "PrimalDualVector() >> Mismatched primal vector. " + - "Must be DesignVector!") + "PrimalDualVector() >> primal_vec must be a DesignVector!") + else: + self.fail('AssertionError expected') + + try: + PrimalDualVector(self.pv1, self.eq1, self.eq1) + except AssertionError as err: + self.assertEqual( + str(err), + "PrimalDualVector() >> ineq_vec must be a DualVectorINEQ!") else: self.fail('AssertionError expected') @@ -107,8 +134,7 @@ def test_bad_init_args(self): except AssertionError as err: self.assertEqual( str(err), - "PrimalDualVector() >> Mismatched primal vector. " + - "Must be DesignVector!") + "PrimalDualVector() >> primal_vec must be a DesignVector!") else: self.fail('AssertionError expected') @@ -142,22 +168,177 @@ def test_init_guess_case3(self): err = self.eq3.base.data - self.pd_vec3.init_dual*(np.ones(5)) self.assertEqual(np.linalg.norm(err), 0) - def test_opt_residual_case1(self): + def test_init_guess_case4(self): + self.pd_vec4.equals_init_guess() + + err = self.pv4.base.data - 10*np.ones(10) + self.assertEqual(np.linalg.norm(err), 0) + + def test_kkt_conditions_case1(self): # case 1 has both equality and inequality constraints # recall: self.dual = 1, so dCeqdU^T*dual.eq + dCineqdU^Tdual.ineq = 5 + 5 = 10 - dCdU(self.design, self.state).T.product(self.dual, self.adjoint, self.state_work) + dCdU(self.design, self.state).T.product(CompositeDualVector(self.dual_eq, self.dual_ineq), + self.adjoint, self.state_work) # recall: dFdU = -1 self.state_work.equals_objective_partial(self.design, self.state) self.state_work.plus(self.adjoint) self.state_work.times(-1.) # We get adjoint = 9 dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) - self.pd_vec1.equals_opt_residual(self.at_pd1, self.state, self.adjoint) + self.pd_vec1.equals_KKT_conditions(self.at_pd1, self.state, self.adjoint) # check results exp_dLdX_norm = np.sqrt(10. * 20.**2) self.assertEqual(self.pd_vec1.primal.norm2, exp_dLdX_norm) exp_dLdEq_norm = np.sqrt(5. * 200.**2) - self.assertEqual(self.pd_vec1.dual.eq.norm2, exp_dLdEq_norm) - exp_dLdIn_norm = np.sqrt(5. * 119402**2) - self.assertEqual(self.pd_vec1.dual.ineq.norm2, exp_dLdIn_norm) + self.assertEqual(self.pd_vec1.eq.norm2, exp_dLdEq_norm) + exp_dLdIn_norm = np.sqrt(5. * 200**2) + self.assertEqual(self.pd_vec1.ineq.norm2, exp_dLdIn_norm) + + def test_kkt_conditions_case2(self): + # case 2 has inequality constraints only + # recall: self.dual = 1, so dCineqdU^Tdual.ineq = 5 = 5 + dCdU(self.design, self.state).T.product(self.dual_ineq, self.adjoint, self.state_work) + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.plus(self.adjoint) + self.state_work.times(-1.) + # We get adjoint = 4 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + self.pd_vec2.equals_KKT_conditions(self.at_pd2, self.state, self.adjoint) + + # check results + exp_dLdX_norm = np.sqrt(10. * 10.**2) + self.assertEqual(self.pd_vec2.primal.norm2, exp_dLdX_norm) + exp_dLdEq_norm = np.sqrt(5. * 200.**2) + self.assertEqual(self.pd_vec2.ineq.norm2, exp_dLdEq_norm) + + def test_kkt_conditions_case3(self): + # case 3 has equality constraints only + # recall: self.dual = 1, so dCeqdU^Tdual.eq = 5 = 5 + dCdU(self.design, self.state).T.product(self.dual_eq, self.adjoint, self.state_work) + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.plus(self.adjoint) + self.state_work.times(-1.) + # We get adjoint = 4 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + self.pd_vec3.equals_KKT_conditions(self.at_pd3, self.state, self.adjoint) + + # check results + exp_dLdX_norm = np.sqrt(10. * 10.**2) + self.assertEqual(self.pd_vec3.primal.norm2, exp_dLdX_norm) + exp_dLdEq_norm = np.sqrt(5. * 200**2) + self.assertEqual(self.pd_vec3.eq.norm2, exp_dLdEq_norm) + + def test_kkt_conditions_case4(self): + # case 4 has no constraints + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.times(-1.) + # We get adjoint = -1 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + self.pd_vec4.equals_KKT_conditions(self.at_pd4, self.state, self.adjoint) + + # check results + exp_dLdX_norm = np.sqrt(10. * 0.**2) + self.assertEqual(self.pd_vec4.primal.norm2, exp_dLdX_norm) + + def test_homotopy_residual_case1(self): + # case 1 has both equality and inequality constraints + # recall: self.dual = 1, so dCeqdU^T*dual.eq + dCineqdU^Tdual.ineq = 5 + 5 = 10 + dCdU(self.design, self.state).T.product(CompositeDualVector(self.dual_eq, self.dual_ineq), + self.adjoint, self.state_work) + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.plus(self.adjoint) + self.state_work.times(-1.) + # We get adjoint = 9 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + dLdx = self.pd_work11 + dLdx.equals_KKT_conditions(self.at_pd1, self.state, self.adjoint) + + init = self.pd_work12 + init.equals_init_guess() + init.eq.equals_constraints(init.primal, self.state) + init.ineq.equals_constraints(init.primal, self.state) + + x = self.at_pd1 + self.pd_vec1.equals_homotopy_residual(dLdx, x, init, mu=0.5) + + # check results + exp_dLdX_norm = np.sqrt(10. * 10.**2) + self.assertAlmostEqual(self.pd_vec1.primal.norm2, exp_dLdX_norm, places=10) + exp_dLdEq_norm = np.sqrt(5. * 100.**2) + self.assertAlmostEqual(self.pd_vec1.eq.norm2, exp_dLdEq_norm, places=10) + exp_dLdIn_norm = np.sqrt(5. * 29701.95**2) + self.assertAlmostEqual(self.pd_vec1.ineq.norm2, exp_dLdIn_norm, places=10) + + def test_homotopy_residual_case2(self): + # case 2 has inequality constraints only + # recall: self.dual = 1, so dCineqdU^Tdual.ineq = 5 + dCdU(self.design, self.state).T.product(self.dual_ineq, self.adjoint, self.state_work) + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.plus(self.adjoint) + self.state_work.times(-1.) + # We get adjoint = 9 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + dLdx = self.pd_work21 + dLdx.equals_KKT_conditions(self.at_pd2, self.state, self.adjoint) + init = self.pd_work22 + init.equals_init_guess() + init.ineq.equals_constraints(init.primal, self.state) + + x = self.at_pd2 + self.pd_vec2.equals_homotopy_residual(dLdx, x, init, mu=0.5) + + # check results + exp_dLdX_norm = np.sqrt(10. * 5**2) + self.assertAlmostEqual(self.pd_vec2.primal.norm2, exp_dLdX_norm, places=10) + exp_dLdIn_norm = np.sqrt(5. * 29701.95**2) + self.assertAlmostEqual(self.pd_vec2.ineq.norm2, exp_dLdIn_norm, places=10) + + def test_homotopy_residual_case3(self): + # case 3 has equality constraints only + # recall: self.dual = 1, so dCeqdU^Tdual.ineq = 5 + dCdU(self.design, self.state).T.product(self.dual_eq, self.adjoint, self.state_work) + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.plus(self.adjoint) + self.state_work.times(-1.) + # We get adjoint = 4 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + dLdx = self.pd_work31 + dLdx.equals_KKT_conditions(self.at_pd3, self.state, self.adjoint) + init = self.pd_work32 + init.equals_init_guess() + init.eq.equals_constraints(init.primal, self.state) + + x = self.at_pd3 + self.pd_vec3.equals_homotopy_residual(dLdx, x, init, mu=0.5) + + # check results + exp_dLdX_norm = np.sqrt(10. * 5**2) + self.assertAlmostEqual(self.pd_vec3.primal.norm2, exp_dLdX_norm, places=10) + exp_dLdEq_norm = np.sqrt(5. * 100.**2) + self.assertAlmostEqual(self.pd_vec3.eq.norm2, exp_dLdEq_norm, places=10) + + def test_homotopy_residual_case4(self): + # case 4 has no constraints + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.times(-1.) + # We get adjoint = -1 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + dLdx = self.pd_work41 + dLdx.equals_KKT_conditions(self.at_pd4, self.state, self.adjoint) + init = self.pd_work42 + init.equals_init_guess() + + x = self.at_pd4 + self.pd_vec4.equals_homotopy_residual(dLdx, x, init, mu=0.5) + + # check results + exp_dLdX_norm = np.sqrt(10. * 0**2) + self.assertAlmostEqual(self.pd_vec3.primal.norm2, exp_dLdX_norm, places=10) From 2afa6163126af0defa91212a790074531760e450 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Mon, 3 Jul 2017 17:16:09 -0400 Subject: [PATCH 05/18] added get_base_data method to extract numpy data from PrimalDualVectors --- src/kona/linalg/vectors/composite.py | 85 +++++------------------- src/kona/test/test_primal_dual_vector.py | 36 ++++++++++ 2 files changed, 51 insertions(+), 70 deletions(-) diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index d9c075a..5497b94 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -362,81 +362,26 @@ def equals_homotopy_residual(self, dLdx, x, init, mu=1.0): self.ineq.equals_mangasarian(self.ineq, x.ineq) self.ineq.minus((1.0 - mu)*0.1) - def equals_opt_residual(self, x, state, adjoint, obj_scale=1.0, cnstr_scale=1.0): + def get_base_data(self, A): """ - Calculates the following nonlinear vector function: - - .. math:: - r(x,\\lambda_h,\\lambda_g) = - \\begin{bmatrix} - \\nabla_x f(x, u) - \\nabla_x h(x, u)^T \\lambda_{h} - \\nabla_x g(x, u)^T \\lambda_{g} \\\\ - h(x,u) \\\\ - -\\theta(|g(x,u) - \\lambda_g|) + \\theta(g(x,u)) + \theta(\\lambda_g) - \\end{bmatrix} - - where :math:`h(x,u)` are the equality constraints, and :math:`g(x,u)` are the - inequality constraints. The vectors :math:`\\lambda_h` and :math:`\\lambda_g` - are the associated Lagrange multipliers. The function :math:`\\theta(z)` - is any strictly increasing function; typically :math:`\\theta(z) = z^3`. The - solution to :math:`r(x,\\lambda_h,\\lambda_g) = 0` is equivalent to the - first-order optimality conditions for the optimization problem - - .. math:: - \\begin{align*} - \\min_x &f(x,u(x)) \\\\ - \\textsf{s.t.} &h(x,u(x)) = 0, \\\\ - &g(x,u(x)) \geq 0. - \\end{align*} + Inserts the PrimalDualVector's underlying data into the given array Parameters ---------- - x : PrimalDualVector - Evaluate first-order optimality conditions at this primal-dual point. - state : StateVector - Evaluate first-order optimality conditions at this state point. - adjoint : StateVector - Evaluate first-order optimality conditions using this adjoint vector. - obj_scale : float, optional - Scaling for the objective function. - cnstr_scale : float, optional - Scaling for the constraints. + A : numpy array + Array into which data is inserted. """ - assert isinstance(x, PrimalDualVector), \ - "PrimalDualVector() >> invalid type x in equals_opt_residual. " + \ - "x vector must be a PrimalDualVector!" - if isinstance(x.dual, CompositeDualVector): - dual_eq = x.dual.eq - dual_ineq = x.dual.ineq - elif isinstance(x.dual, DualVectorINEQ): - dual_eq = None - dual_ineq = x.dual - elif isinstance(x.dual, DualVectorEQ): - dual_eq = x.dual - dual_ineq = None - else: - raise AssertionError("PrimalDualVector() >> invalid dual vector type") - design = x.primal - - design_opt = self.primal - - # first include the objective partial and adjoint contribution - design_opt.equals_total_gradient(design, state, adjoint, obj_scale) - if dual_eq is not None: - # add the Lagrange multiplier products for equality constraints - design_opt.base.data[:] += design_opt._memory.solver.multiply_dCEQdX_T( - design.base.data, state.base, dual_eq.base.data) * \ - cnstr_scale - if dual_ineq is not None: - # add the Lagrange multiplier products for inequality constraints - design_opt.base.data[:] += design_opt._memory.solver.multiply_dCINdX_T( - design.base.data, state.base, dual_ineq.base.data) * \ - cnstr_scale - # include constraint terms - self.dual.equals_constraints(design, state, cnstr_scale) - if isinstance(self.dual, DualVectorINEQ): - self.dual.equals_mangasarian(self.dual, dual_ineq) - elif isinstance(self.dual, CompositeDualVector): - self.dual.ineq.equals_mangasarian(self.dual.ineq, dual_ineq) + # check the length of A against the number of design variables and constraints + ptr = 0 + A[ptr:ptr+self.primal._memory.ndv] = self.primal.base.data[:] + ptr += self.primal._memory.ndv + if self.eq is not None: + print self.eq._memory.neq + print self.eq.base.data.shape + A[ptr:ptr+self.eq._memory.neq] = self.eq.base.data[:] + ptr += self.eq._memory.neq + if self.ineq is not None: + A[ptr:ptr+self.ineq._memory.nineq] = self.ineq.base.data[:] class ReducedKKTVector(CompositeVector): """ diff --git a/src/kona/test/test_primal_dual_vector.py b/src/kona/test/test_primal_dual_vector.py index 466f95b..c2a2cc5 100644 --- a/src/kona/test/test_primal_dual_vector.py +++ b/src/kona/test/test_primal_dual_vector.py @@ -342,3 +342,39 @@ def test_homotopy_residual_case4(self): # check results exp_dLdX_norm = np.sqrt(10. * 0**2) self.assertAlmostEqual(self.pd_vec3.primal.norm2, exp_dLdX_norm, places=10) + + def test_get_base_data_case1(self): + # case 1 has both equality and inequality constraints + A = np.zeros((10+5+5,1)) + self.at_pd1.get_base_data(A[:,0]) + B = np.ones_like(A) + B[0:10] *= 10 + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) + + def test_get_base_data_case2(self): + # case 2 has inequality constraints only + A = np.zeros((10+5,1)) + self.at_pd2.get_base_data(A[:,0]) + B = np.ones_like(A) + B[0:10] *= 10 + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) + + def test_get_base_data_case3(self): + # case 3 has equality constraints only + A = np.zeros((10+5,1)) + self.at_pd3.get_base_data(A[:,0]) + B = np.ones_like(A) + B[0:10] *= 10 + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) + + def test_get_base_data_case4(self): + # case 4 has no constraints + A = np.zeros((10,1)) + self.at_pd4.get_base_data(A[:,0]) + B = np.ones_like(A) + B[0:10] *= 10 + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) From 21439d859980e762cb9f431c2ea8bac4cd1dde76 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Wed, 5 Jul 2017 14:52:33 -0400 Subject: [PATCH 06/18] added a set_base_data method to PrimalDualVector --- src/kona/linalg/vectors/composite.py | 40 ++++++++++--- src/kona/test/test_primal_dual_vector.py | 73 ++++++++++++++++++++---- 2 files changed, 94 insertions(+), 19 deletions(-) diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index 5497b94..844b1de 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -240,6 +240,17 @@ def __init__(self, primal_vec, eq_vec=None, ineq_vec=None): self.ineq = ineq_vec super(PrimalDualVector, self).__init__([primal_vec, eq_vec, ineq_vec]) + def get_num_var(self): + """ + Returns the total number of variables in the PrimalDualVector + """ + nvar = self.primal._memory.ndv + if self.eq is not None: + nvar += self.eq._memory.neq + if self.ineq is not None: + nvar += self.ineq._memory.nineq + return nvar + def equals_init_guess(self): """ Sets the primal-dual vector to the initial guess, using the initial design. @@ -290,13 +301,13 @@ def equals_KKT_conditions(self, x, state, adjoint, obj_scale=1.0, cnstr_scale=1. # first include the objective partial and adjoint contribution dLdx.equals_total_gradient(design, state, adjoint, obj_scale) if x.eq is not None: - # add the Lagrange multiplier products for equality constraints - dLdx.base.data[:] += dLdx._memory.solver.multiply_dCEQdX_T( + # subtract the Lagrange multiplier products for equality constraints + dLdx.base.data[:] -= dLdx._memory.solver.multiply_dCEQdX_T( design.base.data, state.base, x.eq.base.data) * \ cnstr_scale if x.ineq is not None: - # add the Lagrange multiplier products for inequality constraints - dLdx.base.data[:] += dLdx._memory.solver.multiply_dCINdX_T( + # subract the Lagrange multiplier products for inequality constraints + dLdx.base.data[:] -= dLdx._memory.solver.multiply_dCINdX_T( design.base.data, state.base, x.ineq.base.data) * \ cnstr_scale # include constraint terms @@ -371,18 +382,33 @@ def get_base_data(self, A): A : numpy array Array into which data is inserted. """ - # check the length of A against the number of design variables and constraints ptr = 0 A[ptr:ptr+self.primal._memory.ndv] = self.primal.base.data[:] ptr += self.primal._memory.ndv if self.eq is not None: - print self.eq._memory.neq - print self.eq.base.data.shape A[ptr:ptr+self.eq._memory.neq] = self.eq.base.data[:] ptr += self.eq._memory.neq if self.ineq is not None: A[ptr:ptr+self.ineq._memory.nineq] = self.ineq.base.data[:] + def set_base_data(self, A): + """ + Copies the given array into the PrimalDualVector's underlying data + + Parameters + ---------- + A : numpy array + Array that is copied into the PrimalDualVector. + """ + ptr = 0 + self.primal.base.data[:] = A[ptr:ptr+self.primal._memory.ndv] + ptr += self.primal._memory.ndv + if self.eq is not None: + self.eq.base.data[:] = A[ptr:ptr+self.eq._memory.neq] + ptr += self.eq._memory.neq + if self.ineq is not None: + self.ineq.base.data[:] = A[ptr:ptr+self.ineq._memory.nineq] + class ReducedKKTVector(CompositeVector): """ A composite vector representing a combined primal and dual vectors. diff --git a/src/kona/test/test_primal_dual_vector.py b/src/kona/test/test_primal_dual_vector.py index c2a2cc5..33b1be9 100644 --- a/src/kona/test/test_primal_dual_vector.py +++ b/src/kona/test/test_primal_dual_vector.py @@ -174,6 +174,18 @@ def test_init_guess_case4(self): err = self.pv4.base.data - 10*np.ones(10) self.assertEqual(np.linalg.norm(err), 0) + def test_get_num_var_case1(self): + self.assertEqual(self.pd_vec1.get_num_var(), 20) + + def test_get_num_var_case2(self): + self.assertEqual(self.pd_vec2.get_num_var(), 15) + + def test_get_num_var_case3(self): + self.assertEqual(self.pd_vec3.get_num_var(), 15) + + def test_get_num_var_case4(self): + self.assertEqual(self.pd_vec4.get_num_var(), 10) + def test_kkt_conditions_case1(self): # case 1 has both equality and inequality constraints # recall: self.dual = 1, so dCeqdU^T*dual.eq + dCineqdU^Tdual.ineq = 5 + 5 = 10 @@ -181,9 +193,9 @@ def test_kkt_conditions_case1(self): self.adjoint, self.state_work) # recall: dFdU = -1 self.state_work.equals_objective_partial(self.design, self.state) - self.state_work.plus(self.adjoint) + self.state_work.minus(self.adjoint) # L = f - ceq^T*lambda_eq - cineq^T*lambda_ineq self.state_work.times(-1.) - # We get adjoint = 9 + # We get adjoint = -11 dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) self.pd_vec1.equals_KKT_conditions(self.at_pd1, self.state, self.adjoint) @@ -201,9 +213,9 @@ def test_kkt_conditions_case2(self): dCdU(self.design, self.state).T.product(self.dual_ineq, self.adjoint, self.state_work) # recall: dFdU = -1 self.state_work.equals_objective_partial(self.design, self.state) - self.state_work.plus(self.adjoint) + self.state_work.minus(self.adjoint) # L = f - cineq^T*lambda_ineq self.state_work.times(-1.) - # We get adjoint = 4 + # We get adjoint = -6 dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) self.pd_vec2.equals_KKT_conditions(self.at_pd2, self.state, self.adjoint) @@ -219,9 +231,9 @@ def test_kkt_conditions_case3(self): dCdU(self.design, self.state).T.product(self.dual_eq, self.adjoint, self.state_work) # recall: dFdU = -1 self.state_work.equals_objective_partial(self.design, self.state) - self.state_work.plus(self.adjoint) + self.state_work.minus(self.adjoint) # L = f - cineq^T*lambda_ineq self.state_work.times(-1.) - # We get adjoint = 4 + # We get adjoint = -6 dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) self.pd_vec3.equals_KKT_conditions(self.at_pd3, self.state, self.adjoint) @@ -251,9 +263,9 @@ def test_homotopy_residual_case1(self): self.adjoint, self.state_work) # recall: dFdU = -1 self.state_work.equals_objective_partial(self.design, self.state) - self.state_work.plus(self.adjoint) + self.state_work.minus(self.adjoint) # L = f - ceq^T*lambda_eq - cineq^T*lambda_ineq self.state_work.times(-1.) - # We get adjoint = 9 + # We get adjoint = -11 dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) dLdx = self.pd_work11 dLdx.equals_KKT_conditions(self.at_pd1, self.state, self.adjoint) @@ -280,9 +292,9 @@ def test_homotopy_residual_case2(self): dCdU(self.design, self.state).T.product(self.dual_ineq, self.adjoint, self.state_work) # recall: dFdU = -1 self.state_work.equals_objective_partial(self.design, self.state) - self.state_work.plus(self.adjoint) + self.state_work.minus(self.adjoint) # L = f - cineq^T*lambda_ineq self.state_work.times(-1.) - # We get adjoint = 9 + # We get adjoint = -6 dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) dLdx = self.pd_work21 dLdx.equals_KKT_conditions(self.at_pd2, self.state, self.adjoint) @@ -305,9 +317,9 @@ def test_homotopy_residual_case3(self): dCdU(self.design, self.state).T.product(self.dual_eq, self.adjoint, self.state_work) # recall: dFdU = -1 self.state_work.equals_objective_partial(self.design, self.state) - self.state_work.plus(self.adjoint) + self.state_work.minus(self.adjoint) # L = f - ceq^T*lambda_eq self.state_work.times(-1.) - # We get adjoint = 4 + # We get adjoint = -6 dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) dLdx = self.pd_work31 dLdx.equals_KKT_conditions(self.at_pd3, self.state, self.adjoint) @@ -373,8 +385,45 @@ def test_get_base_data_case3(self): def test_get_base_data_case4(self): # case 4 has no constraints A = np.zeros((10,1)) + print A.shape self.at_pd4.get_base_data(A[:,0]) B = np.ones_like(A) B[0:10] *= 10 for i in range(A.shape[0]): self.assertEqual(A[i], B[i]) + + def test_set_base_data_case1(self): + # case 1 has both equality and inequality constraints + A = np.random.random((10+5+5,1)) + self.pd_vec1.set_base_data(A[:,0]) + B = np.zeros_like(A) + self.pd_vec1.get_base_data(B[:,0]) + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) + + def test_set_base_data_case2(self): + # case 2 has only inequality constraints + A = np.random.random((10+5,1)) + self.pd_vec2.set_base_data(A[:,0]) + B = np.zeros_like(A) + self.pd_vec2.get_base_data(B[:,0]) + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) + + def test_set_base_data_case3(self): + # case 3 has only equality constraints + A = np.random.random((10+5,1)) + self.pd_vec3.set_base_data(A[:,0]) + B = np.zeros_like(A) + self.pd_vec3.get_base_data(B[:,0]) + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) + + def test_set_base_data_case4(self): + # case 4 has no constraints + A = np.random.random((10,1)) + self.pd_vec4.set_base_data(A[:,0]) + B = np.zeros_like(A) + self.pd_vec4.get_base_data(B[:,0]) + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) From 0b1038827bd875fa3f33d6a8c554fa0826bb994a Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Thu, 6 Jul 2017 13:51:13 -0400 Subject: [PATCH 07/18] added class (and tests) for Anderson acceleration --- src/kona/linalg/matrices/hessian/__init__.py | 1 + .../matrices/hessian/anderson_multisecant.py | 104 ++++ src/kona/linalg/matrices/hessian/basic.py | 22 +- src/kona/linalg/vectors/composite.py | 3 +- src/kona/test/test_multi_secant.py | 540 ++++++++++++++++++ 5 files changed, 659 insertions(+), 11 deletions(-) create mode 100644 src/kona/linalg/matrices/hessian/anderson_multisecant.py create mode 100644 src/kona/test/test_multi_secant.py diff --git a/src/kona/linalg/matrices/hessian/__init__.py b/src/kona/linalg/matrices/hessian/__init__.py index aa11d4d..17c6898 100644 --- a/src/kona/linalg/matrices/hessian/__init__.py +++ b/src/kona/linalg/matrices/hessian/__init__.py @@ -2,6 +2,7 @@ from lbfgs import LimitedMemoryBFGS from lsr1 import LimitedMemorySR1 +from anderson_multisecant import AndersonMultiSecant from reduced_hessian import ReducedHessian from reduced_kkt import ReducedKKTMatrix diff --git a/src/kona/linalg/matrices/hessian/anderson_multisecant.py b/src/kona/linalg/matrices/hessian/anderson_multisecant.py new file mode 100644 index 0000000..5b14fe3 --- /dev/null +++ b/src/kona/linalg/matrices/hessian/anderson_multisecant.py @@ -0,0 +1,104 @@ +import numpy +from kona.linalg.matrices.hessian.basic import MultiSecantApprox +from kona.linalg.vectors.composite import CompositeDualVector, PrimalDualVector + +class AndersonMultiSecant(MultiSecantApprox): + """ + Anderson-Acceleration multi-secant approximation. + + Attributes + ---------- + init : PrimalDualVector + init.primal is the initial design, while init.eq and init.ineq are the initial constraints + work : PrimalDualVector + Work vector for constructing the difference matrices + """ + + def __init__(self, vector_factory, optns=None): + super(AndersonMultiSecant, self).__init__(vector_factory, optns) + # the number of vectors needed is num(x_hist) + num(r_hist) + num(x_diff) + num(r_diff) + + # num(init + work) = 2*max_stored + 2*(max_stored - 1) + 2 = 4*max_stored + self.primal_factory.request_num_vectors(4*self.max_stored) + if self.eq_factory is not None: + self.eq_factory.request_num_vectors(4*self.max_stored) + if self.ineq_factory is not None: + self.ineq_factory.request_num_vectors(4*self.max_stored) + self.init = None + self.work = None + + def _generate_vector(self): + """ + Create appropriate vector based on vector factory + """ + assert self.primal_factory is not None, \ + 'AndersonMultiSecant() >> primal_factory is not defined!' + dual_eq = None + dual_ineq = None + primal = self.primal_factory.generate() + if self.eq_factory is not None: + dual_eq = self.eq_factory.generate() + if self.ineq_factory is not None: + dual_ineq = self.ineq_factory.generate() + return PrimalDualVector(primal, eq_vec=dual_eq, ineq_vec=dual_ineq) + + def set_initial_data(self, init_data): + """ + Defines initial design (self.init.primal) and constraints (self.init.eq, self.init.ineq) + """ + self.init = self._generate_vector() + self.init.equals(init_data) + self.work = self._generate_vector() + + def add_to_history(self, x_new, r_new): + assert len(self.x_hist) == len(self.r_hist), \ + 'AndersonMultiSecant() >> inconsistent list sizes!' + if len(self.x_hist) == self.max_stored: + # we are out of memory, so pop the oldest vectors + self.x_hist.popleft() + self.r_hist.popleft() + # get new vectors to store the correction + x = self._generate_vector() + r = self._generate_vector() + # copy the input vectors + x.equals(x_new) + r.equals(r_new) + # insert into the lists + self.x_hist.append(x) + self.r_hist.append(r) + + def build_difference_matrices(self, mu=0.0): + assert len(self.x_hist) == len(self.r_hist), \ + 'AndersonMultiSecant() >> inconsistent list sizes!' + # clear the previous x_diff and r_diff lists + del self.x_diff[:], self.r_diff[:] + self.work.equals_homotopy_residual(self.r_hist[0], self.x_hist[0], self.init, mu=mu) + for k in range(len(self.x_hist)-1): + # generate new vectors for x_diff and r_diff lists + dx = self._generate_vector() + dr = self._generate_vector() + self.x_diff.append(dx) + self.r_diff.append(dr) + # now get the differences + self.x_diff[k].equals_ax_p_by(1.0, self.x_hist[k+1], -1.0, self.x_hist[k]) + self.r_diff[k].equals(self.work) + self.work.equals_homotopy_residual(self.r_hist[k+1], self.x_hist[k+1], self.init, mu=mu) + self.r_diff[k].minus(self.work) + self.r_diff[k].times(-1.0) + + def solve(self, in_vec, out_vec, rel_tol=1e-15): + # store the difference matrices and rhs in numpy array format + nvar = in_vec.get_num_var() + dR = numpy.empty((nvar, len(self.x_diff))) + dX = numpy.empty_like(dR) + for k, vec in enumerate(self.r_diff): + vec.get_base_data(dR[:,k]) + for k, vec in enumerate(self.x_diff): + vec.get_base_data(dX[:,k]) + rhs = numpy.empty((nvar)) + in_vec.get_base_data(rhs) + dRinv = numpy.linalg.pinv(dR, rcond=1e-14) + sol = numpy.empty_like(rhs) + sol[:] = -rhs - numpy.matmul(dX - dR, numpy.matmul(dRinv,rhs)) + out_vec.set_base_data(sol) + +# imports at the bottom to prevent circular import errors diff --git a/src/kona/linalg/matrices/hessian/basic.py b/src/kona/linalg/matrices/hessian/basic.py index 60ed2e3..2711d57 100644 --- a/src/kona/linalg/matrices/hessian/basic.py +++ b/src/kona/linalg/matrices/hessian/basic.py @@ -1,3 +1,4 @@ +from collections import deque class BaseHessian(object): """ @@ -39,7 +40,11 @@ def __init__(self, vector_factory, optns=None): self.state_factory = None self.eq_factory = None self.ineq_factory = None - if type(self.vec_fac) is list: + if type(self.vec_fac) is not list: + assert self.vec_fac._vec_type is DesignVector, \ + 'BaseHessian() >> non-list factory is not for DesignVectors!' + self.primal_factory = vector_factory + else: for factory in self.vec_fac: if factory is not None: if factory._vec_type is DesignVector: @@ -128,8 +133,6 @@ class MultiSecantApprox(BaseHessian): ---------- max_stored : int Maximum number of previous steps and residuals stored. - ptr : list of integers - Integer "pointers" that indicate the oldest to newest data x_hist : list of KonaVector The sequence of solutions used to build x_diff r_hist : list of KonaVector @@ -141,13 +144,12 @@ class MultiSecantApprox(BaseHessian): """ def __init__(self, vector_factory, optns={}): - assert isinstance(vector_factory, VectorFactory), \ - "MultiSecantApprox() >> Invalid vector factory!" + #assert isinstance(vector_factory, VectorFactory), \ + # "MultiSecantApprox() >> Invalid vector factory!" super(MultiSecantApprox, self).__init__(vector_factory, optns) self.max_stored = get_opt(optns, 10, 'max_stored') - self.ptr = [] - self.x_hist = [] - self.r_hist = [] + self.x_hist = deque() + self.r_hist = deque() self.x_diff = [] self.r_diff = [] @@ -164,7 +166,7 @@ def add_to_history(self, x_new, r_new): """ raise NotImplementedError # pragma: no cover - def build_difference_matrices(mu=0.0): + def build_difference_matrices(self, mu=0.0): """ Constructs the solution and residual differences from the history @@ -174,7 +176,7 @@ def build_difference_matrices(mu=0.0): Homotopy continuation parameter """ raise NotImplementedError # pragma: no cover - + # imports at the bottom to prevent circular import errors import sys from kona.options import get_opt diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index 844b1de..7ff3098 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -238,7 +238,8 @@ def __init__(self, primal_vec, eq_vec=None, ineq_vec=None): self.primal = primal_vec self.eq = eq_vec self.ineq = ineq_vec - super(PrimalDualVector, self).__init__([primal_vec, eq_vec, ineq_vec]) + vec_list = [primal_vec, eq_vec, ineq_vec] + super(PrimalDualVector, self).__init__([vec for vec in vec_list if vec is not None]) def get_num_var(self): """ diff --git a/src/kona/test/test_multi_secant.py b/src/kona/test/test_multi_secant.py new file mode 100644 index 0000000..18b8d4d --- /dev/null +++ b/src/kona/test/test_multi_secant.py @@ -0,0 +1,540 @@ +import unittest +from collections import deque + +import numpy as np + +from kona.linalg.memory import KonaMemory +from kona.user import UserSolver +from kona.linalg.vectors.composite import PrimalDualVector +from kona.linalg.matrices.hessian import AndersonMultiSecant + +class AndersonMultiSecantTestCase(unittest.TestCase): + '''Test cases for Anderson-Acceleration multi-secant method''' + + def test_generate_vector_case1(self): + # case 1: unconstrained + max_stored = 3 + solver = UserSolver(3) + km = KonaMemory(solver) + pf = km.primal_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant(pf, optns) + km.allocate_memory() + + vec = aa._generate_vector() + self.assertTrue(isinstance(vec, PrimalDualVector)) + self.assertTrue(vec.primal is not None) + self.assertTrue(vec.eq is None) + self.assertTrue(vec.ineq is None) + + def test_generate_vector_case2(self): + # case 2: constrained with equality constraints only + max_stored = 3 + solver = UserSolver(10, num_eq=3) + km = KonaMemory(solver) + pf = km.primal_factory + eqf = km.eq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, eqf], optns) + km.allocate_memory() + + vec = aa._generate_vector() + self.assertTrue(isinstance(vec, PrimalDualVector)) + self.assertTrue(vec.primal is not None) + self.assertTrue(vec.eq is not None) + self.assertTrue(vec.ineq is None) + + def test_generate_vector_case3(self): + # case 3: constrained with inequality constraints only + max_stored = 3 + solver = UserSolver(10, num_ineq=5) + km = KonaMemory(solver) + pf = km.primal_factory + ineqf = km.ineq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, ineqf], optns) + km.allocate_memory() + + vec = aa._generate_vector() + self.assertTrue(isinstance(vec, PrimalDualVector)) + self.assertTrue(vec.primal is not None) + self.assertTrue(vec.eq is None) + self.assertTrue(vec.ineq is not None) + + def test_generate_vector_case4(self): + # case 4: constrained with both equality and inequality constraints + max_stored = 3 + solver = UserSolver(10, num_eq=3, num_ineq=5) + km = KonaMemory(solver) + pf = km.primal_factory + eqf = km.eq_factory + ineqf = km.ineq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, eqf, ineqf], optns) + km.allocate_memory() + + vec = aa._generate_vector() + self.assertTrue(isinstance(vec, PrimalDualVector)) + self.assertTrue(vec.primal is not None) + self.assertTrue(vec.eq is not None) + self.assertTrue(vec.ineq is not None) + + def test_set_initial_data_case1(self): + # case 1: unconstrained + max_stored = 3 + solver = UserSolver(3) + km = KonaMemory(solver) + pf = km.primal_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant(pf, optns) + # add request that will be used to create pd_vec + km.primal_factory.request_num_vectors(1) + km.allocate_memory() + + p_vec = km.primal_factory.generate() + pd_vec = PrimalDualVector(p_vec) + A = np.random.random((3,1)) + pd_vec.set_base_data(A[:,0]) + aa.set_initial_data(pd_vec) + B = np.zeros_like(A) + aa.init.get_base_data(B[:,0]) + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) + + def test_set_initial_data_case2(self): + # case 2: constrained with equality constraints only + max_stored = 3 + solver = UserSolver(10, num_eq=3) + km = KonaMemory(solver) + pf = km.primal_factory + eqf = km.eq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, eqf], optns) + # add request that will be used to create pd_vec + km.primal_factory.request_num_vectors(1) + km.eq_factory.request_num_vectors(1) + km.allocate_memory() + + p_vec = km.primal_factory.generate() + eq_vec = km.eq_factory.generate() + pd_vec = PrimalDualVector(p_vec, eq_vec=eq_vec) + A = np.random.random((13,1)) + pd_vec.set_base_data(A[:,0]) + aa.set_initial_data(pd_vec) + B = np.zeros_like(A) + aa.init.get_base_data(B[:,0]) + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) + + def test_set_initial_data_case3(self): + # case 3: constrained with inequality constraints only + max_stored = 3 + solver = UserSolver(10, num_ineq=5) + km = KonaMemory(solver) + pf = km.primal_factory + ineqf = km.ineq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, ineqf], optns) + # add request that will be used to create pd_vec + km.primal_factory.request_num_vectors(1) + km.ineq_factory.request_num_vectors(1) + km.allocate_memory() + + p_vec = km.primal_factory.generate() + ineq_vec = km.ineq_factory.generate() + pd_vec = PrimalDualVector(p_vec, ineq_vec=ineq_vec) + A = np.random.random((15,1)) + pd_vec.set_base_data(A[:,0]) + aa.set_initial_data(pd_vec) + B = np.zeros_like(A) + aa.init.get_base_data(B[:,0]) + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) + + def test_set_initial_data_case4(self): + # case 4: constrained with both equality and inequality constraints + max_stored = 3 + solver = UserSolver(10, num_eq=3, num_ineq=5) + km = KonaMemory(solver) + pf = km.primal_factory + eqf = km.eq_factory + ineqf = km.ineq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, eqf, ineqf], optns) + # add request that will be used to create pd_vec + km.primal_factory.request_num_vectors(1) + km.eq_factory.request_num_vectors(1) + km.ineq_factory.request_num_vectors(1) + km.allocate_memory() + + p_vec = km.primal_factory.generate() + eq_vec = km.eq_factory.generate() + ineq_vec = km.ineq_factory.generate() + pd_vec = PrimalDualVector(p_vec, eq_vec=eq_vec, ineq_vec=ineq_vec) + A = np.random.random((18,1)) + pd_vec.set_base_data(A[:,0]) + aa.set_initial_data(pd_vec) + B = np.zeros_like(A) + aa.init.get_base_data(B[:,0]) + for i in range(A.shape[0]): + self.assertEqual(A[i], B[i]) + + def test_add_to_history_case1(self): + # case 1: unconstrained + max_stored = 3 + solver = UserSolver(3) + km = KonaMemory(solver) + pf = km.primal_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant(pf, optns) + # add request that will be used to create x and r + km.primal_factory.request_num_vectors(2) + km.allocate_memory() + + # create x and r vectors + p_vec1 = km.primal_factory.generate() + x = PrimalDualVector(p_vec1) + p_vec2 = km.primal_factory.generate() + r = PrimalDualVector(p_vec2) + + # Create synthetic solution and residual data and feed into aa history + Xrand = np.random.random((3,max_stored+1)) + Rrand = np.random.random((3,max_stored+1)) + for k in range(max_stored+1): + x.set_base_data(Xrand[:,k]) + r.set_base_data(Rrand[:,k]) + aa.add_to_history(x, r) + # At this point, the first vector from Xrand and Rrand should have been discardded + A = np.zeros(Xrand.shape[0]) + for k in range(max_stored): + aa.x_hist[k].get_base_data(A[:]) + for i in range(Xrand.shape[0]): + self.assertEqual(A[i],Xrand[i,k+1]) + aa.r_hist[k].get_base_data(A[:]) + for i in range(Xrand.shape[0]): + self.assertEqual(A[i],Rrand[i,k+1]) + + def test_add_to_history_case2(self): + # case 2: constrained with equality constraints only + max_stored = 3 + solver = UserSolver(10, num_eq=3) + km = KonaMemory(solver) + pf = km.primal_factory + eqf = km.eq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, eqf], optns) + # add request that will be used to create x and r + km.primal_factory.request_num_vectors(2) + km.eq_factory.request_num_vectors(2) + km.allocate_memory() + + # create x and r vectors + p_vec1 = km.primal_factory.generate() + eq_vec1 = km.eq_factory.generate() + x = PrimalDualVector(p_vec1, eq_vec=eq_vec1) + p_vec2 = km.primal_factory.generate() + eq_vec2 = km.eq_factory.generate() + r = PrimalDualVector(p_vec2, eq_vec=eq_vec2) + + # Create synthetic solution and residual data and feed into aa history + Xrand = np.random.random((13,max_stored+1)) + Rrand = np.random.random((13,max_stored+1)) + for k in range(max_stored+1): + x.set_base_data(Xrand[:,k]) + r.set_base_data(Rrand[:,k]) + aa.add_to_history(x, r) + # At this point, the first vector from Xrand and Rrand should have been discardded + A = np.zeros(Xrand.shape[0]) + for k in range(max_stored): + aa.x_hist[k].get_base_data(A[:]) + for i in range(Xrand.shape[0]): + self.assertEqual(A[i],Xrand[i,k+1]) + aa.r_hist[k].get_base_data(A[:]) + for i in range(Xrand.shape[0]): + self.assertEqual(A[i],Rrand[i,k+1]) + + def test_add_to_history_case3(self): + # case 3: constrained with inequality constraints only + max_stored = 3 + solver = UserSolver(10, num_ineq=5) + km = KonaMemory(solver) + pf = km.primal_factory + ineqf = km.ineq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, ineqf], optns) + # add request that will be used to create x and r + km.primal_factory.request_num_vectors(2) + km.ineq_factory.request_num_vectors(2) + km.allocate_memory() + + # create x and r vectors + p_vec1 = km.primal_factory.generate() + ineq_vec1 = km.ineq_factory.generate() + x = PrimalDualVector(p_vec1, ineq_vec=ineq_vec1) + p_vec2 = km.primal_factory.generate() + ineq_vec2 = km.ineq_factory.generate() + r = PrimalDualVector(p_vec2, ineq_vec=ineq_vec2) + + # Create synthetic solution and residual data and feed into aa history + Xrand = np.random.random((15,max_stored+1)) + Rrand = np.random.random((15,max_stored+1)) + for k in range(max_stored+1): + x.set_base_data(Xrand[:,k]) + r.set_base_data(Rrand[:,k]) + aa.add_to_history(x, r) + # At this point, the first vector from Xrand and Rrand should have been discardded + A = np.zeros(Xrand.shape[0]) + for k in range(max_stored): + aa.x_hist[k].get_base_data(A[:]) + for i in range(Xrand.shape[0]): + self.assertEqual(A[i],Xrand[i,k+1]) + aa.r_hist[k].get_base_data(A[:]) + for i in range(Xrand.shape[0]): + self.assertEqual(A[i],Rrand[i,k+1]) + + def test_add_to_history_case4(self): + # case 4: constrained with both equality and inequality constraints + max_stored = 3 + solver = UserSolver(10, num_eq=3, num_ineq=5) + km = KonaMemory(solver) + pf = km.primal_factory + eqf = km.eq_factory + ineqf = km.ineq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, eqf, ineqf], optns) + # add request that will be used to create x and r + km.primal_factory.request_num_vectors(2) + km.eq_factory.request_num_vectors(2) + km.ineq_factory.request_num_vectors(2) + km.allocate_memory() + + # create x and r vectors + p_vec1 = km.primal_factory.generate() + eq_vec1 = km.eq_factory.generate() + ineq_vec1 = km.ineq_factory.generate() + x = PrimalDualVector(p_vec1, eq_vec=eq_vec1, ineq_vec=ineq_vec1) + p_vec2 = km.primal_factory.generate() + eq_vec2 = km.eq_factory.generate() + ineq_vec2 = km.ineq_factory.generate() + r = PrimalDualVector(p_vec2, eq_vec=eq_vec2, ineq_vec=ineq_vec2) + + # Create synthetic solution and residual data and feed into aa history + Xrand = np.random.random((18,max_stored+1)) + Rrand = np.random.random((18,max_stored+1)) + for k in range(max_stored+1): + x.set_base_data(Xrand[:,k]) + r.set_base_data(Rrand[:,k]) + aa.add_to_history(x, r) + # At this point, the first vector from Xrand and Rrand should have been discardded + A = np.zeros(Xrand.shape[0]) + for k in range(max_stored): + aa.x_hist[k].get_base_data(A[:]) + for i in range(Xrand.shape[0]): + self.assertEqual(A[i],Xrand[i,k+1]) + aa.r_hist[k].get_base_data(A[:]) + for i in range(Xrand.shape[0]): + self.assertEqual(A[i],Rrand[i,k+1]) + + def test_build_difference_matrices_case1(self): + # case 1: unconstrained + max_stored = 3 + solver = UserSolver(3) + km = KonaMemory(solver) + pf = km.primal_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant(pf, optns) + # add request that will be used to create x and r + km.primal_factory.request_num_vectors(2) + km.allocate_memory() + + # create x and r vectors + p_vec1 = km.primal_factory.generate() + x = PrimalDualVector(p_vec1) + p_vec2 = km.primal_factory.generate() + r = PrimalDualVector(p_vec2) + aa.set_initial_data(x) + + # Create synthetic solution and residual data and feed into aa history + x_scalar = np.linspace(1, max_stored, num=max_stored) + r_scalar = np.ones_like(x_scalar) + for k in range(max_stored): + x.equals(x_scalar[k]) + r.equals(r_scalar[k]) + aa.add_to_history(x, r) + aa.build_difference_matrices() + A = np.zeros(3) + for k in range(max_stored-1): + aa.x_diff[k].get_base_data(A[:]) + for i in range(A.shape[0]): + self.assertAlmostEqual(A[i], 1.0, places=14) + + def test_build_difference_matrices_case2(self): + # case 2: constrained with equality constraints only + max_stored = 3 + solver = UserSolver(10, num_eq=3) + km = KonaMemory(solver) + pf = km.primal_factory + eqf = km.eq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, eqf], optns) + # add request that will be used to create x and r + km.primal_factory.request_num_vectors(2) + km.eq_factory.request_num_vectors(2) + km.allocate_memory() + + # create x and r vectors + p_vec1 = km.primal_factory.generate() + eq_vec1 = km.eq_factory.generate() + x = PrimalDualVector(p_vec1, eq_vec=eq_vec1) + p_vec2 = km.primal_factory.generate() + eq_vec2 = km.eq_factory.generate() + r = PrimalDualVector(p_vec2, eq_vec=eq_vec2) + aa.set_initial_data(x) + + # Create synthetic solution and residual data and feed into aa history + x_scalar = np.linspace(1, max_stored, num=max_stored) + r_scalar = np.ones_like(x_scalar) + for k in range(max_stored): + x.equals(x_scalar[k]) + r.equals(r_scalar[k]) + aa.add_to_history(x, r) + aa.build_difference_matrices() + A = np.zeros(13) + for k in range(max_stored-1): + aa.x_diff[k].get_base_data(A[:]) + for i in range(A.shape[0]): + self.assertAlmostEqual(A[i], 1.0, places=14) + + def test_build_difference_matrices_case3(self): + # case 3: constrained with inequality constraints only + max_stored = 3 + solver = UserSolver(10, num_ineq=5) + km = KonaMemory(solver) + pf = km.primal_factory + ineqf = km.ineq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, ineqf], optns) + # add request that will be used to create x and r + km.primal_factory.request_num_vectors(2) + km.ineq_factory.request_num_vectors(2) + km.allocate_memory() + + # create x and r vectors + p_vec1 = km.primal_factory.generate() + ineq_vec1 = km.ineq_factory.generate() + x = PrimalDualVector(p_vec1, ineq_vec=ineq_vec1) + p_vec2 = km.primal_factory.generate() + ineq_vec2 = km.ineq_factory.generate() + r = PrimalDualVector(p_vec2, ineq_vec=ineq_vec2) + aa.set_initial_data(x) + + # Create synthetic solution and residual data and feed into aa history + x_scalar = np.linspace(1, max_stored, num=max_stored) + r_scalar = np.ones_like(x_scalar) + for k in range(max_stored): + x.equals(x_scalar[k]) + r.equals(r_scalar[k]) + aa.add_to_history(x, r) + aa.build_difference_matrices() + A = np.zeros(15) + for k in range(max_stored-1): + aa.x_diff[k].get_base_data(A[:]) + for i in range(A.shape[0]): + self.assertAlmostEqual(A[i], 1.0, places=14) + + def test_build_difference_matrices_case4(self): + # case 4: constrained with both equality and inequality constraints + max_stored = 3 + solver = UserSolver(10, num_eq=3, num_ineq=5) + km = KonaMemory(solver) + pf = km.primal_factory + eqf = km.eq_factory + ineqf = km.ineq_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant([pf, eqf, ineqf], optns) + # add request that will be used to create x and r + km.primal_factory.request_num_vectors(2) + km.eq_factory.request_num_vectors(2) + km.ineq_factory.request_num_vectors(2) + km.allocate_memory() + + # create x and r vectors + p_vec1 = km.primal_factory.generate() + eq_vec1 = km.eq_factory.generate() + ineq_vec1 = km.ineq_factory.generate() + x = PrimalDualVector(p_vec1, eq_vec=eq_vec1, ineq_vec=ineq_vec1) + p_vec2 = km.primal_factory.generate() + eq_vec2 = km.eq_factory.generate() + ineq_vec2 = km.ineq_factory.generate() + r = PrimalDualVector(p_vec2, eq_vec=eq_vec2, ineq_vec=ineq_vec2) + aa.set_initial_data(x) + + # Create synthetic solution and residual data and feed into aa history + x_scalar = np.linspace(1, max_stored, num=max_stored) + r_scalar = np.ones_like(x_scalar) + for k in range(max_stored): + x.equals(x_scalar[k]) + r.equals(r_scalar[k]) + aa.add_to_history(x, r) + aa.build_difference_matrices() + A = np.zeros(18) + for k in range(max_stored-1): + aa.x_diff[k].get_base_data(A[:]) + for i in range(A.shape[0]): + self.assertAlmostEqual(A[i], 1.0, places=14) + + def test_solve_case1(self): + # case 1: unconstrained + max_stored = 4 + solver = UserSolver(3) + km = KonaMemory(solver) + pf = km.primal_factory + optns = {'max_stored': max_stored} + aa = AndersonMultiSecant(pf, optns) + # add request that will be used to create x and r + km.primal_factory.request_num_vectors(2) + km.allocate_memory() + + # create x and r vectors + p_vec1 = km.primal_factory.generate() + x = PrimalDualVector(p_vec1) + p_vec2 = km.primal_factory.generate() + r = PrimalDualVector(p_vec2) + aa.set_initial_data(x) + + # For the following synthetic iterates: + # Hessian matrix is [1 0 0; 0 100 0; 0 0 10] + # initial iterate is [1 1 1] + x.set_base_data(np.array([1.,1.,1.])) + r.set_base_data(np.array([1.,100.,10.])) + aa.add_to_history(x, r) + x.set_base_data(np.array([0.,1.,1.])) + r.set_base_data(np.array([0.,100.,10.])) + aa.add_to_history(x, r) + x.set_base_data(np.array([0.,0.,1.])) + r.set_base_data(np.array([0.,0.,10.])) + aa.add_to_history(x, r) + x.set_base_data(np.array([0.,0.,0.])) + r.set_base_data(np.array([0.,0.,0.])) + aa.add_to_history(x, r) + aa.build_difference_matrices(mu=1.0) + + # AA should be able to recover exact inverse + x.set_base_data(np.array([1.,0.,0.])) + aa.solve(x, r) + A = np.zeros(3) + r.get_base_data(A[:]) + self.assertAlmostEqual(np.linalg.norm(A - np.array([-1.,0.,0.])), 0.0, places=14) + + x.set_base_data(np.array([0.,100.,0.])) + aa.solve(x, r) + r.get_base_data(A[:]) + self.assertAlmostEqual(np.linalg.norm(A - np.array([0.,-1.,0.])), 0.0, places=14) + + x.set_base_data(np.array([0.,0.,10.])) + aa.solve(x, r) + r.get_base_data(A[:]) + self.assertAlmostEqual(np.linalg.norm(A - np.array([0.,0.,-1.])), 0.0, places=14) + +if __name__ == "__main__": + unittest.main() From fc83c01363e50a5b87e6ba6ff57071935328ff97 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Fri, 7 Jul 2017 16:56:01 -0400 Subject: [PATCH 08/18] preliminary version of multi-secant method --- src/kona/algorithms/__init__.py | 1 + src/kona/algorithms/homotopy_multi_secant.py | 287 +++++++++++++++++++ src/kona/linalg/common.py | 15 +- src/kona/linalg/vectors/common.py | 68 +++++ src/kona/linalg/vectors/composite.py | 95 +++++- src/kona/test/test_primal_dual_vector.py | 84 +++++- 6 files changed, 537 insertions(+), 13 deletions(-) create mode 100644 src/kona/algorithms/homotopy_multi_secant.py diff --git a/src/kona/algorithms/__init__.py b/src/kona/algorithms/__init__.py index 5ed6bd3..b4def12 100644 --- a/src/kona/algorithms/__init__.py +++ b/src/kona/algorithms/__init__.py @@ -5,5 +5,6 @@ from composite_step_rsnk import CompositeStepRSNK from predictor_corrector import PredictorCorrector from predictor_corrector_cnstr import PredictorCorrectorCnstr +from homotopy_multi_secant import HomotopyMultiSecant from verifier import Verifier import util diff --git a/src/kona/algorithms/homotopy_multi_secant.py b/src/kona/algorithms/homotopy_multi_secant.py new file mode 100644 index 0000000..fcb489e --- /dev/null +++ b/src/kona/algorithms/homotopy_multi_secant.py @@ -0,0 +1,287 @@ +from kona.algorithms.base_algorithm import OptimizationAlgorithm + + +class HomotopyMultiSecant(OptimizationAlgorithm): + """ + Algorithm for generic nonlinearly constrained optiimzation problems. + + Parameters + ---------- + primal_factory : VectorFactory + state_factory : VectorFactory + eq_factory : VectorFactory + ineq_factory: VectorFactory + optns : dict, optional + """ + + def __init__(self, primal_factory, state_factory, + eq_factory, ineq_factory, optns=None): + # trigger base class initialization + super(HomotopyMultiSecant, self).__init__( + primal_factory, state_factory, eq_factory, ineq_factory, optns + ) + + # number of vectors required in solve() method + # PrimalDualVectors = X, X_init, dX, R, dLdX, dXdmu, dXdmu_old + num_pd = 7 + self.primal_factory.request_num_vectors(num_pd) + if self.eq_factory is not None: + self.eq_factory.request_num_vectors(num_pd) + if self.ineq_factory is not None: + self.ineq_factory.request_num_vectors(num_pd) + self.state_factory.request_num_vectors(3) + + # iteration counters + self.iter = 0 + self.inner = 0 + + # set the type of multi-secant method + try: + multisecant = get_opt( + self.optns, AndersonMultiSecant, 'multi_secant', 'type') + hessian_optns = get_opt(self.optns, {}, 'multi_secant') + hessian_optns['out_file'] = self.info_file + self.multisecant = multisecant([self.primal_factory, self.eq_factory, + self.ineq_factory], hessian_optns) + except Exception: + raise BadKonaOption(self.optns, 'multi_secant','type') + + # homotopy options + self.mu = get_opt(self.optns, 0.0, 'homotopy', 'mu_init') + self.inner_grad_tol = get_opt(self.optns, 1e-2, 'homotopy', 'inner_grad_tol') + self.inner_feas_tol = get_opt(self.optns, 1e-2, 'homotopy', 'inner_feas_tol') + self.inner_max_iter = get_opt(self.optns, 10, 'homotopy', 'inner_maxiter') + self.dmu = get_opt(self.optns, 0.1, 'homotopy', 'init_step') + self.target_angle = get_opt(self.optns, 5.0*np.pi/180., 'homotopy', 'nominal_angle') + self.dmu_max = get_opt(self.optns, 0.1, 'homotopy', 'max_step') + self.dmu_min = get_opt(self.optns, 0.001, 'homotopy', 'min_step') + self.radius_max = get_opt(self.optns, 1.0, 'homotopy', 'radius_max') + + # The following data members are set by super class + # self.primal_tol + # self.cnstr_tol + # self.max_iter + # self.info_file + # self.hist_file + + def _write_header(self): + self.hist_file.write( + '# Kona homotopy-globalized multi-secant convergence history file\n' + + '# outer' + ' '*5 + + ' inner' + ' '*5 + + ' cost' + ' '*5 + + 'optimality ' + ' '*5 + + 'feasibility ' + ' '*5 + + 'objective ' + ' '*5 + + 'mu param ' + '\n' + ) + + def _write_history(self, opt, feas, obj): + self.hist_file.write( + '%7i'%self.iter + ' '*5 + + '%7i'%self.inner + ' '*5 + + '%7i'%self.primal_factory._memory.cost + ' '*5 + + '%11e'%opt + ' '*5 + + '%11e'%feas + ' '*5 + + '%11e'%obj + ' '*5 + + '%11e'%self.mu + '\n' + ) + + def _generate_vector(self): + """ + Create appropriate vector based on vector factory + """ + assert self.primal_factory is not None, \ + 'AndersonMultiSecant() >> primal_factory is not defined!' + dual_eq = None + dual_ineq = None + primal = self.primal_factory.generate() + if self.eq_factory is not None: + dual_eq = self.eq_factory.generate() + if self.ineq_factory is not None: + dual_ineq = self.ineq_factory.generate() + return PrimalDualVector(primal, eq_vec=dual_eq, ineq_vec=dual_ineq) + + def solve(self): + self._write_header() + self.info_file.write( + '\n' + + '**************************************************\n' + + '*** Using Multi-Secant Algorithm ***\n' + + '**************************************************\n' + + '\n') + + # generate primal-dual vectors + X = self._generate_vector() + X_init = self._generate_vector() + dX = self._generate_vector() + R = self._generate_vector() + dLdX = self._generate_vector() + dXdmu = self._generate_vector() + dXdmu_old = self._generate_vector() + + # generate state vectors + state = self.state_factory.generate() + state_work = self.state_factory.generate() + adjoint = self.state_factory.generate() + + # evaluate the initial design, state, and adjoint before starting outer iterations + X.equals_init_guess() + X_init.equals(X) + state.equals_primal_solution(X.primal) + adjoint.equals_homotopy_adjoint(X, state, state_work) + + # send initial point info to the user + solver_info = current_solution(self.iter, X.primal, state, adjoint, X.get_dual()) + if isinstance(solver_info, str): + self.info_file.write('\n' + solver_info + '\n') + + # initialize the multi-secant method + dLdX.equals_KKT_conditions(X, state, adjoint) + self.multisecant.add_to_history(X, dLdX) + dXdmu.equals(dLdX) + dXdmu.primal.equals(X_init.primal) + self.multisecant.set_initial_data(dXdmu) + + # initialize the tangent vectors + dXdmu.equals(0.) + dXdmu_old.equals(0.) + + # get initial optimality and feasibility + R.equals_homotopy_residual(dLdX, X, X_init, mu=1.0) + opt, feas = R.get_optimality_and_feasiblity() + grad_norm0 = max(opt, EPS) + feas_norm0 = max(feas, EPS) + self.info_file.write( + 'grad_norm0 = %e\n'%grad_norm0 + + 'feas_norm0 = %e\n'%feas_norm0) + + # Begin outer (homotopy) loop here + converged = False + for i in xrange(self.max_iter): + self.iter += 1 + self.inner = 0 + + # evaluate optimality and feasibility, and output as necessary + # Note: dLdX should be up-to-date at this point + R.equals_homotopy_residual(dLdX, X, X_init, mu=1.0) + grad_norm, feas_norm = R.get_optimality_and_feasiblity() + self.info_file.write( + '==========================================================\n' + + 'Outer (Homotopy) Iteration %i: mu = %e\n\n'%(self.iter,self.mu)) + + self.info_file.write( + 'grad_norm = %e (%e <-- tolerance)\n'%( + grad_norm, self.primal_tol) + + 'feas_norm = %e (%e <-- tolerance)\n'%( + feas_norm, self.cnstr_tol)) + + # write convergence history + obj_val = objective_value(X.primal, state) + self._write_history(grad_norm, feas_norm, obj_val) + + # check for convergence + if grad_norm < self.primal_tol*grad_norm0 and \ + feas_norm < self.cnstr_tol*feas_norm0: + converged = True + break + + # find the predictor step + R.equals_predictor_rhs(dLdX, X, X_init, mu=self.mu) + self.multisecant.build_difference_matrices(mu=self.mu) + dXdmu_old.equals(dXdmu) + self.multisecant.solve(R, dXdmu) + # adjust the step size + if self.iter > 1: + denom = np.sqrt(dXdmu.norm2**2 + 1)*np.sqrt(dXdmu_old.norm2**2 + 1) + numer = (dXdmu.inner(dXdmu_old) + 1.) + phi = math.acos(max(-1., min(1., numer/denom))) + if phi > EPS: + self.dmu /= (phi/self.target_angle) + self.dmu = min(min(self.dmu_max, self.radius_max/dXdmu.norm2), self.dmu) + self.dmu = max(self.dmu_min, self.dmu) + # take the predictor step + X.equals_ax_p_by(1.0, X, self.dmu, dXdmu) + self.mu = min(self.mu + self.dmu, 1.0) + + # evaluate the optimality residual and add to multisecant + state.equals_primal_solution(X.primal) + adjoint.equals_homotopy_adjoint(X, state, state_work) + dLdX.equals_KKT_conditions(X, state, adjoint) + self.multisecant.add_to_history(X, dLdX) + + R.equals_homotopy_residual(dLdX, X, X_init, mu=self.mu) + opt, feas = R.get_optimality_and_feasiblity() + inner_grad_norm0 = opt + inner_feas_norm0 = max(feas, EPS) + inner_grad_norm = opt + inner_feas_norm = feas + self.info_file.write( + '-------------------- Corrector Iterations\n' + + 'inner_grad_norm0 = %e'%inner_grad_norm0 + + ': inner_feas_norm0 = %e\n'%inner_feas_norm0) + + # inner (corrector) iterations + for k in xrange(self.inner_max_iter): + self.inner += 1 + self.info_file.write( + 'Corr. Iter. %i'%self.inner + + ': inner_grad_norm = %e'%inner_grad_norm + + ': inner_feas_norm = %e\n'%inner_feas_norm) + + # check for convergence + if np.fabs(self.mu - 1.0) < EPS: + # check for outer convergence if mu = 1.0 + if grad_norm < self.primal_tol*grad_norm0 and \ + feas_norm < self.cnstr_tol*feas_norm0: + converged = True + break + else: # check for inner convergence + if inner_grad_norm < self.inner_grad_tol*inner_grad_norm0 and \ + inner_feas_norm < self.inner_feas_tol*inner_feas_norm0: + break + + # solve for corrector step + self.multisecant.build_difference_matrices(mu=self.mu) + self.multisecant.solve(R, dX) + # safe-guard against large steps + if dX.norm2 > self.radius_max: + dX.times(self.radius_max/dX.norm2) + X.plus(dX) + + # evaluate at new X and store data + state.equals_primal_solution(X.primal) + obj_val = objective_value(X.primal, state) + adjoint.equals_homotopy_adjoint(X, state, state_work) + dLdX.equals_KKT_conditions(X, state, adjoint) + self.multisecant.add_to_history(X, dLdX) + R.equals_homotopy_residual(dLdX, X, X_init, mu=1.0) + grad_norm, feas_norm = R.get_optimality_and_feasiblity() + R.equals_homotopy_residual(dLdX, X, X_init, mu=self.mu) + inner_grad_norm, inner_feas_norm = R.get_optimality_and_feasiblity() + # write convergence history + self._write_history(grad_norm, feas_norm, obj_val) + + # end of inner (corrector) iterations + + # send initial point info to the user + solver_info = current_solution(self.iter, X.primal, state, adjoint, X.get_dual()) + if isinstance(solver_info, str): + self.info_file.write('\n' + solver_info + '\n') + + # end of outer (homotopy) iterations + if converged: + self.info_file.write('Optimization successful!\n') + else: + self.info_file.write('Optimization FAILED!\n') + self.info_file.write( + 'Total number of nonlinear iterations: %i\n\n'%self.iter) + +# imports here to prevent circular errors +import math +import numpy as np +from kona.options import BadKonaOption, get_opt +from kona.linalg.common import current_solution, objective_value +from kona.linalg.vectors.composite import PrimalDualVector +from kona.linalg.matrices.hessian import AndersonMultiSecant +from kona.linalg.solvers.util import EPS diff --git a/src/kona/linalg/common.py b/src/kona/linalg/common.py index 3b1e951..ea3d1c4 100644 --- a/src/kona/linalg/common.py +++ b/src/kona/linalg/common.py @@ -39,13 +39,20 @@ def current_solution(num_iter, curr_primal, curr_state=None, curr_adj=None, elif isinstance(curr_primal, DesignVector): curr_design = curr_primal out_slack = None - if isinstance(curr_dual, DualVectorEQ): + if isinstance(curr_dual, DualVectorINEQ): + out_eq = None + out_ineq = deepcopy(curr_dual.base.data) + elif isinstance(curr_dual, DualVectorEQ): out_eq = deepcopy(curr_dual.base.data) + out_ineq = None + elif isinstance(curr_dual, CompositeDualVector): + out_eq = deepcopy(curr_dual.eq.base.data) + out_ineq = deepcopy(curr_dual.ineq.base.data) elif curr_dual is None: out_eq = None + out_ineq = None else: - raise TypeError("Invalid dual vector type: must be DualVectorEQ!") - out_ineq = None + raise TypeError("Invalid dual vector type") else: raise TypeError("Invalid primal vector type: " + @@ -193,4 +200,4 @@ def factor_linear_system(at_primal, at_state): solver = at_design._memory.solver - solver.factor_linear_system(at_design.base.data, at_state.base) \ No newline at end of file + solver.factor_linear_system(at_design.base.data, at_state.base) diff --git a/src/kona/linalg/vectors/common.py b/src/kona/linalg/vectors/common.py index 3e330d3..3a9d026 100644 --- a/src/kona/linalg/vectors/common.py +++ b/src/kona/linalg/vectors/common.py @@ -580,6 +580,49 @@ def equals_lagrangian_adjoint(self, at_kkt, at_state, state_work, # solve the adjoint system dRdU(at_primal, at_state).T.solve(state_work, self, rel_tol=1e-6) + def equals_homotopy_adjoint(self, at_pd, at_state, state_work, + obj_scale=1.0, cnstr_scale=1.0, rel_tol=1e-6): + """ + Computes in-place the adjoint variables for the augmented Lagrangian, + :math:`L = f - h^T \lambda_h - g^T \lambda_g` where :math:`h` denotes the equality constraints (if any) and :math:`g` denotes the inequality constraints (if any). Note the + negative signs, which are opposite to how the Lagrangian is defined elsewhere in Kona. + + Parameters + ---------- + at_pd : PrimalDualVector + Current primal-dual solution. + at_state : StateVector + Current state point. + state_work : StateVector + Temporary work vector of State type. + obj_scale : float, optional + Scaling for the objective function. + cnstr_scale : float, optional + Scaling for the constraints. + rel_tol : float, optional + relative tolerance provided to adjoint solver + """ + # assemble the right hand side for the Lagrangian adjoint solve + assert isinstance(at_pd, PrimalDualVector), \ + "Invalid KKT vector: must be PrimalDualVector!" + at_primal = at_pd.primal + + # add the dual contribution to the Adjoint RHS, if necessary + at_dual = at_pd.get_dual() + if at_dual is not None: + dCdU(at_primal, at_state).T.product(at_dual, self, state_work) + self.times(cnstr_scale) + else: + self.equals(0.) + # add objective partial + state_work.equals_objective_partial(at_primal, at_state) + state_work.times(obj_scale) + # form the adjoint RHS + state_work.minus(self) # note the minus here !!! + state_work.times(-1.) + # solve the adjoint system + dRdU(at_primal, at_state).T.solve(state_work, self, rel_tol=rel_tol) + class DualVector(KonaVector): pass @@ -713,8 +756,33 @@ def equals_mangasarian(self, ineq, mult): self.base.data[:] = -np.fabs(ineq.base.data[:] - mult.base.data[:])**3 \ + ineq.base.data[:]**3 + mult.base.data[:]**3 + def deriv_mangasarian(self, ineq, mult): + """ + Evaluate the derivative of the Mangasarian function with respect to the first argument, + :math:`\\partial G/\\partial g = -3\\mathsf{sgn}(g-\\lambda_g)|g - \\lambda_g|^2 + 3(g)^2`, + where :math:`\\mathsf{sgn}(\\cdot)` is the sign of the argument. Stores the result + in-place. + + Parameters + ---------- + ineq : DualVectorINEQ + Inequality constraint at the current design and state. + mult : DualVectorINEQ + The multiplier corresponding to the inequality constraints. + """ + assert isinstance(ineq, DualVectorINEQ), \ + 'DualVectorINEQ >> Invalid type for ineq. ' + \ + 'Must be DualVecINEQ!' + assert isinstance(mult, DualVectorINEQ), \ + 'DualVectorINEQ >> Invalid type for mult. ' + \ + 'Must be DualVecINEQ!' + self.base.data[:] = np.multiply(-3*np.sign(ineq.base.data[:] - mult.base.data[:]), + (ineq.base.data[:] - mult.base.data[:])**2) \ + + 3*ineq.base.data[:]**2 + # package imports at the bottom to prevent import errors import numpy as np from kona.linalg.vectors.composite import ReducedKKTVector +from kona.linalg.vectors.composite import PrimalDualVector from kona.linalg.vectors.composite import CompositePrimalVector from kona.linalg.matrices.common import * diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index 7ff3098..05234db 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -252,6 +252,20 @@ def get_num_var(self): nvar += self.ineq._memory.nineq return nvar + def get_dual(self): + """ + Returns 1) eq or ineq if only one is None, 2) a CompositeDualVector if neither is none or + 3) None if both are none + """ + dual = None + if self.eq is not None and self.ineq is not None: + dual = CompositeDualVector(self.eq, self.ineq) + elif self.eq is not None: + dual = self.eq + elif self.ineq is not None: + dual = self.ineq + return dual + def equals_init_guess(self): """ Sets the primal-dual vector to the initial guess, using the initial design. @@ -317,19 +331,33 @@ def equals_KKT_conditions(self, x, state, adjoint, obj_scale=1.0, cnstr_scale=1. if cineq is not None: cineq.equals_constraints(design, state ,cnstr_scale) + def get_optimality_and_feasiblity(self): + """ + Returns the norm of the primal (opt) the dual parts of the vector (feas). If the dual + parts of the vector are both None, then feas is returned as zero. + """ + opt = self.primal.norm2 + feas = 0.0 + if self.eq is not None: + feas += self.eq.norm2**2 + if self.ineq is not None: + feas += self.ineq.norm2**2 + feas = np.sqrt(feas) + return opt, feas + def equals_homotopy_residual(self, dLdx, x, init, mu=1.0): """ Using dLdx=:math:`\\begin{pmatrix} \\nabla_x L && h && g \\end{pmatrix}`, which can be obtained from the method equals_KKT_conditions, as well as the initial values init=:math:`\\begin{pmatrix} x_0 && h(x_0,u_0) && g(x_0,u_0) \\end{pmatrix}` and the - current point :math:`\\begin{pmatrix} x && \lambda_h && \lambda_g \end{pmatrix}`, this + current point x=:math:`\\begin{pmatrix} x && \lambda_h && \lambda_g \end{pmatrix}`, this method computes the following nonlinear vector function: .. math:: r(x,\\lambda_h,\\lambda_g;\\mu) = \\begin{bmatrix} \\mu\\left[\\nabla_x f(x, u) - \\nabla_x h(x, u)^T \\lambda_{h} - \\nabla_x g(x, u)^T \\lambda_{g}\\right] + (1 - \\mu)(x - x_0) \\\\ - h(x,u) - (1-\mu)h(x_0,u_0) \\\\ + -\\mu h(x,u) - (1-\mu)\\lambda_h \\\\ -|g(x,u) - (1-\\mu)*g_0 - \\lambda_g|^3 + (g(x,u) - (1-\\mu)g_0)^3 + \\lambda_g^3 - (1-\\mu)\\hat{g} \\end{bmatrix} @@ -366,14 +394,73 @@ def equals_homotopy_residual(self, dLdx, x, init, mu=1.0): self.primal.equals_ax_p_by(1.0, x.primal, -1.0, init.primal) self.primal.equals_ax_p_by(mu, dLdx.primal, (1.0 - mu), self.primal) if dLdx.eq is not None: - # include the equality constraint part: h - (1-mu)h_0 - self.eq.equals_ax_p_by(1.0, dLdx.eq, -(1.0 - mu), init.eq) + # include the equality constraint part: -mu*h - (1-mu)*lambda_h + self.eq.equals_ax_p_by(-mu, dLdx.eq, -(1.0 - mu), x.eq) + #self.eq.equals_ax_p_by(mu, dLdx.eq, -(1.0 - mu), x.eq) if dLdx.ineq is not None: # include the inequality constraint part self.ineq.equals_ax_p_by(1.0, dLdx.ineq, -(1.0 - mu), init.ineq) self.ineq.equals_mangasarian(self.ineq, x.ineq) self.ineq.minus((1.0 - mu)*0.1) + def equals_predictor_rhs(self, dLdx, x, init, mu=1.0): + """ + Using dLdx=:math:`\\begin{pmatrix} \\nabla_x L && h && g \\end{pmatrix}`, which can be + obtained from the method equals_KKT_conditions, as well as the initial values + init=:math:`\\begin{pmatrix} x_0 && h(x_0,u_0) && g(x_0,u_0) \\end{pmatrix}` and the + current point x=:math:`\\begin{pmatrix} x && \lambda_h && \lambda_g \end{pmatrix}`, this + method computes the right-hand-side for the homotopy-predictor step, that is + + .. math:: + \partial r/\partial \\mu = + \\begin{bmatrix} + \\left[\\nabla_x f(x, u) - \\nabla_x h(x, u)^T \\lambda_{h} - \\nabla_x g(x, u)^T \\lambda_{g}\\right] - (x - x_0) \\\\ + -h(x,u) + \\lambda_h \\\\ + -3*(g_0)|g(x,u) - (1-\\mu)*g_0 - \\lambda_g|^2 + 3*g_0*(g(x,u) - (1-\\mu)g_0)^2 + \hat{g} + \\end{bmatrix} + + where :math:`h(x,u)` are the equality constraints, and :math:`g(x,u)` are the + inequality constraints. The vectors :math:`\\lambda_h` and :math:`\\lambda_g` + are the associated Lagrange multipliers. + + Parameters + ---------- + dLdx : PrimalDualVector + The total derivative of the Lagranginan with respect to the primal and dual variables. + x : PrimalDualVector + The current solution vector value corresponding to dLdx. + init: PrimalDualVector + The initial primal variable, as well as the initial constraint values. + mu : float + Homotopy parameter; must be between 0 and 1. + """ + assert isinstance(dLdx, PrimalDualVector), \ + "PrimalDualVector() >> dLdx must be a PrimalDualVector!" + assert isinstance(init, PrimalDualVector), \ + "PrimalDualVector() >> init must be a PrimalDualVector!" + if dLdx.eq is None or self.eq is None or x.eq is None or init.eq is None: + assert dLdx.eq is None and self.eq is None and x.eq is None and init.eq is None, \ + "PrimalDualVector() >> inconsistent eq component in self, dLdx, x, and/or init!" + if dLdx.ineq is None or self.ineq is None or x.ineq is None or init.ineq is None: + assert dLdx.ineq is None and self.ineq is None and x.ineq is None \ + and init.ineq is None, \ + "PrimalDualVector() >> inconsistent ineq component in self, dLdx, x, and/or init!" + if dLdx == self or x == self or init == self: + "PrimalDualVector() >> equals_predictor_rhs is not in-place!" + # construct the primal part of the rhs: dLdx - (x - x_0) + self.primal.equals_ax_p_by(1.0, dLdx.primal, -1.0, x.primal) + self.primal.plus(init.primal) + if dLdx.eq is not None: + # include the equality constraint part: -h + lambda_h + self.eq.equals_ax_p_by(-1., dLdx.eq, 1.0, x.eq) + #self.eq.equals_ax_p_by(1., dLdx.eq, 1.0, x.eq) + if dLdx.ineq is not None: + # include the inequality constraint part + self.ineq.equals_ax_p_by(1.0, dLdx.ineq, -(1.0 - mu), init.ineq) + self.ineq.deriv_mangasarian(self.ineq, x.ineq) + self.ineq.times(init.ineq) + self.ineq.plus(0.1) + def get_base_data(self, A): """ Inserts the PrimalDualVector's underlying data into the given array diff --git a/src/kona/test/test_primal_dual_vector.py b/src/kona/test/test_primal_dual_vector.py index 33b1be9..fe364ab 100644 --- a/src/kona/test/test_primal_dual_vector.py +++ b/src/kona/test/test_primal_dual_vector.py @@ -4,6 +4,7 @@ from kona.linalg.memory import KonaMemory from dummy_solver import DummySolver from kona.linalg.matrices.common import * +from kona.linalg.vectors.common import DualVectorEQ, DualVectorINEQ from kona.linalg.vectors.composite import CompositePrimalVector from kona.linalg.vectors.composite import CompositeDualVector from kona.linalg.vectors.composite import PrimalDualVector @@ -14,21 +15,24 @@ def setUp(self): solver = DummySolver(10, 5, 5) # 10 DVs, 5 equality, 5 inequality self.km = km = KonaMemory(solver) - km.primal_factory.request_num_vectors(10) + km.primal_factory.request_num_vectors(11) km.state_factory.request_num_vectors(3) - km.eq_factory.request_num_vectors(6) - km.ineq_factory.request_num_vectors(10) # doe we need this many? + km.eq_factory.request_num_vectors(7) + km.ineq_factory.request_num_vectors(11) # doe we need this many? km.allocate_memory() # get some work vectors self.primal_work1 = km.primal_factory.generate() self.primal_work2 = km.primal_factory.generate() + self.primal_work3 = km.primal_factory.generate() self.slack_work = km.ineq_factory.generate() # only used for error testing self.state_work = km.state_factory.generate() self.eq_work1 = km.eq_factory.generate() self.eq_work2 = km.eq_factory.generate() + self.eq_work3 = km.eq_factory.generate() self.ineq_work1 = km.ineq_factory.generate() self.ineq_work2 = km.ineq_factory.generate() + self.ineq_work3 = km.ineq_factory.generate() # get static eval point vectors self.design = km.primal_factory.generate() @@ -63,6 +67,8 @@ def setUp(self): ineq_vec=self.ineq_work1) self.pd_work12 = PrimalDualVector(self.primal_work2, eq_vec=self.eq_work2, ineq_vec=self.ineq_work2) + self.pd_work13 = PrimalDualVector(self.primal_work3, eq_vec=self.eq_work3, + ineq_vec=self.ineq_work3) # case 2: primal-dual vector has design and ineq self.pv2 = km.primal_factory.generate() @@ -186,6 +192,22 @@ def test_get_num_var_case3(self): def test_get_num_var_case4(self): self.assertEqual(self.pd_vec4.get_num_var(), 10) + def test_get_dual_case1(self): + dual = self.pd_vec1.get_dual() + self.assertTrue(isinstance(dual, CompositeDualVector)) + + def test_get_dual_case2(self): + dual = self.pd_vec2.get_dual() + self.assertTrue(isinstance(dual, DualVectorINEQ)) + + def test_get_dual_case3(self): + dual = self.pd_vec3.get_dual() + self.assertTrue(isinstance(dual, DualVectorEQ)) + + def test_get_dual_case4(self): + dual = self.pd_vec4.get_dual() + self.assertTrue(dual is None) + def test_kkt_conditions_case1(self): # case 1 has both equality and inequality constraints # recall: self.dual = 1, so dCeqdU^T*dual.eq + dCineqdU^Tdual.ineq = 5 + 5 = 10 @@ -281,11 +303,16 @@ def test_homotopy_residual_case1(self): # check results exp_dLdX_norm = np.sqrt(10. * 10.**2) self.assertAlmostEqual(self.pd_vec1.primal.norm2, exp_dLdX_norm, places=10) - exp_dLdEq_norm = np.sqrt(5. * 100.**2) + exp_dLdEq_norm = np.sqrt(5. * 100.5**2) self.assertAlmostEqual(self.pd_vec1.eq.norm2, exp_dLdEq_norm, places=10) exp_dLdIn_norm = np.sqrt(5. * 29701.95**2) self.assertAlmostEqual(self.pd_vec1.ineq.norm2, exp_dLdIn_norm, places=10) + # test get_optimality_and_feasiblity while we are at it + opt, feas = self.pd_vec1.get_optimality_and_feasiblity() + self.assertAlmostEqual(opt, exp_dLdX_norm, places=10) + self.assertAlmostEqual(feas, np.sqrt(exp_dLdEq_norm**2 + exp_dLdIn_norm**2), places=10) + def test_homotopy_residual_case2(self): # case 2 has inequality constraints only # recall: self.dual = 1, so dCineqdU^Tdual.ineq = 5 @@ -311,6 +338,11 @@ def test_homotopy_residual_case2(self): exp_dLdIn_norm = np.sqrt(5. * 29701.95**2) self.assertAlmostEqual(self.pd_vec2.ineq.norm2, exp_dLdIn_norm, places=10) + # test get_optimality_and_feasiblity while we are at it + opt, feas = self.pd_vec2.get_optimality_and_feasiblity() + self.assertAlmostEqual(opt, exp_dLdX_norm, places=10) + self.assertAlmostEqual(feas, exp_dLdIn_norm, places=10) + def test_homotopy_residual_case3(self): # case 3 has equality constraints only # recall: self.dual = 1, so dCeqdU^Tdual.ineq = 5 @@ -333,9 +365,14 @@ def test_homotopy_residual_case3(self): # check results exp_dLdX_norm = np.sqrt(10. * 5**2) self.assertAlmostEqual(self.pd_vec3.primal.norm2, exp_dLdX_norm, places=10) - exp_dLdEq_norm = np.sqrt(5. * 100.**2) + exp_dLdEq_norm = np.sqrt(5. * 100.5**2) self.assertAlmostEqual(self.pd_vec3.eq.norm2, exp_dLdEq_norm, places=10) + # test get_optimality_and_feasiblity while we are at it + opt, feas = self.pd_vec3.get_optimality_and_feasiblity() + self.assertAlmostEqual(opt, exp_dLdX_norm, places=10) + self.assertAlmostEqual(feas, exp_dLdEq_norm, places=10) + def test_homotopy_residual_case4(self): # case 4 has no constraints # recall: dFdU = -1 @@ -355,6 +392,43 @@ def test_homotopy_residual_case4(self): exp_dLdX_norm = np.sqrt(10. * 0**2) self.assertAlmostEqual(self.pd_vec3.primal.norm2, exp_dLdX_norm, places=10) + # test get_optimality_and_feasiblity while we are at it + opt, feas = self.pd_vec4.get_optimality_and_feasiblity() + self.assertAlmostEqual(opt, exp_dLdX_norm, places=10) + self.assertAlmostEqual(feas, 0.0, places=10) + + def test_equals_predictor_rhs_case1(self): + # case 1 has both equality and inequality constraints + # recall: self.dual = 1, so dCeqdU^T*dual.eq + dCineqdU^Tdual.ineq = 5 + 5 = 10 + dCdU(self.design, self.state).T.product(CompositeDualVector(self.dual_eq, self.dual_ineq), + self.adjoint, self.state_work) + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.minus(self.adjoint) # L = f - ceq^T*lambda_eq - cineq^T*lambda_ineq + self.state_work.times(-1.) + # We get adjoint = -11 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + dLdx = self.pd_work11 + dLdx.equals_KKT_conditions(self.at_pd1, self.state, self.adjoint) + + init = self.pd_work12 + init.equals_init_guess() + init.eq.equals_constraints(init.primal, self.state) + init.ineq.equals_constraints(init.primal, self.state) + + x = self.at_pd1 + self.pd_vec1.equals_homotopy_residual(dLdx, x, init, mu=0.5) + dmu = 1e-7 + self.pd_work13.equals_homotopy_residual(dLdx, x, init, mu=0.5+dmu) + self.pd_work13.minus(self.pd_vec1) + self.pd_work13.divide_by(dmu) + self.pd_vec1.equals_predictor_rhs(dLdx, x, init, mu=0.5) + + # check results + self.assertAlmostEqual(self.pd_vec1.primal.norm2, self.pd_work13.primal.norm2, places=5) + self.assertAlmostEqual(self.pd_vec1.eq.norm2/self.pd_work13.eq.norm2, 1.0, places=5) + self.assertAlmostEqual(self.pd_vec1.ineq.norm2/self.pd_work13.ineq.norm2, 1.0, places=5) + def test_get_base_data_case1(self): # case 1 has both equality and inequality constraints A = np.zeros((10+5+5,1)) From 1aa0e9027919dd4184cbb610d817bc77a62c6ee9 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Fri, 14 Jul 2017 09:21:34 -0400 Subject: [PATCH 09/18] introduced ReducedSpaceMultiSecant algorithm --- src/kona/algorithms/__init__.py | 1 + src/kona/algorithms/homotopy_multi_secant.py | 8 +- .../algorithms/reduced_space_multi_secant.py | 261 ++++++++++++++++++ src/kona/examples/sellar.py | 6 + .../matrices/hessian/anderson_multisecant.py | 45 ++- src/kona/linalg/matrices/hessian/basic.py | 7 +- src/kona/linalg/vectors/common.py | 16 +- src/kona/linalg/vectors/composite.py | 61 ++++ src/kona/test/test_dual_vectors.py | 6 +- src/kona/test/test_multi_secant.py | 26 +- src/kona/test/test_primal_dual_vector.py | 105 ++++++- 11 files changed, 517 insertions(+), 25 deletions(-) create mode 100644 src/kona/algorithms/reduced_space_multi_secant.py diff --git a/src/kona/algorithms/__init__.py b/src/kona/algorithms/__init__.py index b4def12..fc1a9c6 100644 --- a/src/kona/algorithms/__init__.py +++ b/src/kona/algorithms/__init__.py @@ -6,5 +6,6 @@ from predictor_corrector import PredictorCorrector from predictor_corrector_cnstr import PredictorCorrectorCnstr from homotopy_multi_secant import HomotopyMultiSecant +from reduced_space_multi_secant import ReducedSpaceMultiSecant from verifier import Verifier import util diff --git a/src/kona/algorithms/homotopy_multi_secant.py b/src/kona/algorithms/homotopy_multi_secant.py index fcb489e..1921a80 100644 --- a/src/kona/algorithms/homotopy_multi_secant.py +++ b/src/kona/algorithms/homotopy_multi_secant.py @@ -92,7 +92,7 @@ def _generate_vector(self): Create appropriate vector based on vector factory """ assert self.primal_factory is not None, \ - 'AndersonMultiSecant() >> primal_factory is not defined!' + 'HomotopyMultiSecant() >> primal_factory is not defined!' dual_eq = None dual_ineq = None primal = self.primal_factory.generate() @@ -107,7 +107,7 @@ def solve(self): self.info_file.write( '\n' + '**************************************************\n' + - '*** Using Multi-Secant Algorithm ***\n' + + '*** Using Multi-Secant Homotopy ***\n' + '**************************************************\n' + '\n') @@ -188,7 +188,7 @@ def solve(self): # find the predictor step R.equals_predictor_rhs(dLdX, X, X_init, mu=self.mu) - self.multisecant.build_difference_matrices(mu=self.mu) + self.multisecant.build_difference_matrices_for_homotopy(mu=self.mu) dXdmu_old.equals(dXdmu) self.multisecant.solve(R, dXdmu) # adjust the step size @@ -242,7 +242,7 @@ def solve(self): break # solve for corrector step - self.multisecant.build_difference_matrices(mu=self.mu) + self.multisecant.build_difference_matrices_for_homotopy(mu=self.mu) self.multisecant.solve(R, dX) # safe-guard against large steps if dX.norm2 > self.radius_max: diff --git a/src/kona/algorithms/reduced_space_multi_secant.py b/src/kona/algorithms/reduced_space_multi_secant.py new file mode 100644 index 0000000..7020c88 --- /dev/null +++ b/src/kona/algorithms/reduced_space_multi_secant.py @@ -0,0 +1,261 @@ +from kona.algorithms.base_algorithm import OptimizationAlgorithm + + +class ReducedSpaceMultiSecant(OptimizationAlgorithm): + """ + Algorithm for generic nonlinearly constrained optiimzation problems. + + Parameters + ---------- + primal_factory : VectorFactory + state_factory : VectorFactory + eq_factory : VectorFactory + ineq_factory: VectorFactory + optns : dict, optional + """ + + def __init__(self, primal_factory, state_factory, + eq_factory, ineq_factory, optns=None): + # trigger base class initialization + super(ReducedSpaceMultiSecant, self).__init__( + primal_factory, state_factory, eq_factory, ineq_factory, optns + ) + + # number of vectors required in solve() method + # PrimalDualVectors = X, dX, R, dLdX, obj_grad, dX_safe + num_pd = 6 + self.primal_factory.request_num_vectors(num_pd) + if self.eq_factory is not None: + self.eq_factory.request_num_vectors(num_pd) + if self.ineq_factory is not None: + self.ineq_factory.request_num_vectors(num_pd) + self.state_factory.request_num_vectors(3) + + # iteration counter + self.iter = 0 + + # set the type of multi-secant method + try: + # the multisecant member is for the "second"-order scheme + multisecant = get_opt( + self.optns, AndersonMultiSecant, 'multi_secant', 'type') + hessian_optns = get_opt(self.optns, {}, 'multi_secant') + hessian_optns['out_file'] = self.info_file + self.multisecant = multisecant([self.primal_factory, self.eq_factory, + self.ineq_factory], hessian_optns) + # the multisecant_safe member is for the first-order scheme + self.multisecant_safe = multisecant([self.primal_factory, self.eq_factory, + self.ineq_factory], hessian_optns) + + except Exception: + raise BadKonaOption(self.optns, 'multi_secant','type') + + self.primal_tol_abs = get_opt(self.optns, 1e-6, 'opt_tol_abs') + self.cnstr_tol_abs = get_opt(self.optns, 1e-6, 'feas_tol_abs') + self.beta = get_opt(self.optns, 1.0, 'multi_secant', 'beta') + + # TEMPORARY!!! + self.radius_max = get_opt(self.optns, 1.0, 'homotopy', 'radius_max') + + # The following data members are set by super class + # self.primal_tol + # self.cnstr_tol + # self.max_iter + # self.info_file + # self.hist_file + + def _write_header(self): + self.hist_file.write( + '# Kona homotopy-globalized multi-secant convergence history file\n' + + '# iter' + ' '*5 + + ' cost' + ' '*5 + + 'optimality ' + ' '*5 + + 'feasibility ' + ' '*5 + + 'objective ' + '\n' + ) + + def _write_history(self, opt, feas, obj, info): + self.hist_file.write( + '%7i'%self.iter + ' '*5 + + '%7i'%self.primal_factory._memory.cost + ' '*5 + + '%11e'%opt + ' '*5 + + '%11e'%feas + ' '*5 + + '%11e'%obj + ' '*5 + + info + '\n' + ) + + def _generate_vector(self): + """ + Create appropriate vector based on vector factory + """ + assert self.primal_factory is not None, \ + 'ReducedSpaceMultiSecant() >> primal_factory is not defined!' + dual_eq = None + dual_ineq = None + primal = self.primal_factory.generate() + if self.eq_factory is not None: + dual_eq = self.eq_factory.generate() + if self.ineq_factory is not None: + dual_ineq = self.ineq_factory.generate() + return PrimalDualVector(primal, eq_vec=dual_eq, ineq_vec=dual_ineq) + + def solve(self): + self._write_header() + self.info_file.write( + '\n' + + '**************************************************\n' + + '*** Using Multi-Secant Algorithm ***\n' + + '**************************************************\n' + + '\n') + + # generate primal-dual vectors, and other vectors + X = self._generate_vector() + dX = self._generate_vector() + dX_safe = self._generate_vector() + R = self._generate_vector() + dLdX = self._generate_vector() + obj_grad = self._generate_vector() + + # generate state vectors + state = self.state_factory.generate() + state_work = self.state_factory.generate() + adjoint = self.state_factory.generate() + + # evaluate the initial design, state, and adjoint before starting outer iterations + X.equals_init_guess() + state.equals_primal_solution(X.primal) + obj_val = objective_value(X.primal, state) + + # initialize the multi-secant methods + self.multisecant.set_initial_data(X) + self.multisecant_safe.set_initial_data(X) + + # get the adjoint for -h^T\lambda_h - g^T\lambda_g + adjoint.equals_homotopy_adjoint(X, state, state_work, obj_scale=0.0) + R.equals_KKT_conditions(X, state, adjoint, obj_scale=0.0) + dLdX.equals(R) + R.primal.plus(X.primal) + self.multisecant_safe.add_to_history(X, R) + + # get the adjoint for df/dx, then add it to finish dLdX computation + adjoint.equals_objective_adjoint(X.primal, state, state_work) + obj_grad.equals(0.) + obj_grad.primal.equals_total_gradient(X.primal, state, adjoint) + dLdX.primal.plus(obj_grad.primal) + self.multisecant.add_to_history(X, dLdX) + + # send initial point info to the user + solver_info = current_solution(self.iter, X.primal, state, adjoint, X.get_dual()) + if isinstance(solver_info, str): + self.info_file.write('\n' + solver_info + '\n') + + # get initial optimality and feasibility + R.equals_primaldual_residual(dLdX, X.ineq) + opt, feas = R.get_optimality_and_feasiblity() + grad_norm0 = max(opt, EPS) + feas_norm0 = max(feas, EPS) + self.info_file.write( + 'grad_norm0 = %e\n'%grad_norm0 + + 'feas_norm0 = %e\n'%feas_norm0) + + # Begin iterations + converged = False + info = ' ' + for i in xrange(self.max_iter): + self.iter += 1 + + # evaluate optimality and feasibility, and output as necessary + # Note: dLdX should be up-to-date at this point + R.equals_primaldual_residual(dLdX, X.ineq) + grad_norm, feas_norm = R.get_optimality_and_feasiblity() + self.info_file.write( + '==========================================================\n' + + 'Iteration %i\n'%self.iter) + self.info_file.write( + 'grad_norm = %e (%e <-- tolerance)\n'%( + grad_norm, self.primal_tol) + + 'feas_norm = %e (%e <-- tolerance)\n'%( + feas_norm, self.cnstr_tol)) + # write convergence history + self._write_history(grad_norm, feas_norm, obj_val, info) + + # check for convergence + if grad_norm < self.primal_tol*grad_norm0 + self.primal_tol_abs and \ + feas_norm < self.cnstr_tol*feas_norm0 + self.cnstr_tol_abs: + converged = True + break + + # get full multi-secant step and projected-gradient + self.multisecant.build_difference_matrices() + self.multisecant.solve(R, dX, self.beta) + self.multisecant_safe.build_difference_matrices() + self.multisecant_safe.solve(obj_grad, dX_safe) + info = ' ' + if dX.primal.inner(dX_safe.primal) < 0.: + # the full step is considered unsafe: compute first-order step + self.multisecant_safe.solve(R, dX, self.beta) + info += 'safe-step' + # remove history from multisecant + #self.multisecant.clear_history() + + # safe-guard against large steps + if dX.norm2 > self.radius_max: #/self.iter: + dX.times(self.radius_max/(dX.norm2)) #*self.iter)) + info += ' length restricted' + X.plus(dX) + + # j = 0 + # while True: + # j += 1 + # X.plus(dX) + # # evaluate at new X, and check if optimality improved + # state.equals_primal_solution(X.primal) + # adjoint.equals_homotopy_adjoint(X, state, state_work) + # dLdX.equals_KKT_conditions(X, state, adjoint) + # R.equals_primaldual_residual(dLdX, X.ineq) + # self.info_file.write( + # 'line search iteration %i: |R(x_k)| = %e: |R(x_k + p_k)| = %e\n'%( + # j, grad_norm**2 + feas_norm**2, R.norm2**2)) + # if R.norm2**2 < grad_norm**2 + feas_norm**2 or j >= 4: + # break + # X.minus(dX) + # dX.times(0.1) + + # evaluate at new X, construct first-order optimality conditions, and store + state.equals_primal_solution(X.primal) + obj_val = objective_value(X.primal, state) + + # get the adjoint for -h^T\lambda_h - g^T\lambda_g + adjoint.equals_homotopy_adjoint(X, state, state_work, obj_scale=0.0) + R.equals_KKT_conditions(X, state, adjoint, obj_scale=0.0) + dLdX.equals(R) + R.primal.plus(X.primal) + self.multisecant_safe.add_to_history(X, R) + + # get the adjoint for df/dx, then add it to finish dLdX computation + adjoint.equals_objective_adjoint(X.primal, state, state_work) + obj_grad.primal.equals_total_gradient(X.primal, state, adjoint) + dLdX.primal.plus(obj_grad.primal) + self.multisecant.add_to_history(X, dLdX) + + # output current solution info to the user + solver_info = current_solution(self.iter, X.primal, state, adjoint, X.get_dual()) + if isinstance(solver_info, str): + self.info_file.write('\n' + solver_info + '\n') + + # end of "Newton" iterations + if converged: + self.info_file.write('Optimization successful!\n') + else: + self.info_file.write('Optimization FAILED!\n') + self.info_file.write( + 'Total number of nonlinear iterations: %i\n\n'%self.iter) + +# imports here to prevent circular errors +import math +import numpy as np +from kona.options import BadKonaOption, get_opt +from kona.linalg.common import current_solution, objective_value +from kona.linalg.vectors.composite import PrimalDualVector +from kona.linalg.matrices.hessian import AndersonMultiSecant +from kona.linalg.solvers.util import EPS diff --git a/src/kona/examples/sellar.py b/src/kona/examples/sellar.py index d56a564..2c185bd 100644 --- a/src/kona/examples/sellar.py +++ b/src/kona/examples/sellar.py @@ -147,6 +147,8 @@ def solve_nonlinear(self, at_design, result): abs_tol = 1e-6 du = 0 while iters < max_iter: + iters += 1 + # print "Seller nonlinear solve iteration",iters y1 = u[0] y2 = u[1] @@ -157,6 +159,7 @@ def solve_nonlinear(self, at_design, result): R[0] = y1 - x1**2 - x3 - x2 + 0.2*y2 R[1] = y2 - np.sqrt(y1) - x1 - x2 + # print "\tnorm(R) = ",np.linalg.norm(R) if np.linalg.norm(R) < abs_tol: result.data[0] = u[0] result.data[1] = u[1] @@ -168,8 +171,11 @@ def solve_nonlinear(self, at_design, result): du = np.linalg.solve(dRdU, -R) + # safe-guard against y1 dropping below 3.16 u += du + u[0] = max(0.01, u[0]) + print 'Sellar nonlinear problem failed to converge' result.data[0] = u[0] result.data[1] = u[1] return -1 diff --git a/src/kona/linalg/matrices/hessian/anderson_multisecant.py b/src/kona/linalg/matrices/hessian/anderson_multisecant.py index 5b14fe3..55de260 100644 --- a/src/kona/linalg/matrices/hessian/anderson_multisecant.py +++ b/src/kona/linalg/matrices/hessian/anderson_multisecant.py @@ -1,3 +1,4 @@ +from collections import deque import numpy from kona.linalg.matrices.hessian.basic import MultiSecantApprox from kona.linalg.vectors.composite import CompositeDualVector, PrimalDualVector @@ -66,7 +67,41 @@ def add_to_history(self, x_new, r_new): self.x_hist.append(x) self.r_hist.append(r) - def build_difference_matrices(self, mu=0.0): + def clear_history(self): + """ + Clears data from x_hist and r_hist deques + """ + self.x_hist = deque() + self.r_hist = deque() + + def build_difference_matrices(self): + assert len(self.x_hist) == len(self.r_hist), \ + 'AndersonMultiSecant() >> inconsistent list sizes!' + # clear the previous x_diff and r_diff lists + del self.x_diff[:], self.r_diff[:] + self.work.equals_primaldual_residual(self.r_hist[0], self.x_hist[0].ineq) + for k in range(len(self.x_hist)-1): + # generate new vectors for x_diff and r_diff lists + dx = self._generate_vector() + dr = self._generate_vector() + self.x_diff.append(dx) + self.r_diff.append(dr) + # now get the differences + self.x_diff[k].equals_ax_p_by(1.0, self.x_hist[k+1], -1.0, self.x_hist[k]) + self.r_diff[k].equals(self.work) + self.work.equals_primaldual_residual(self.r_hist[k+1], self.x_hist[k+1].ineq) + self.r_diff[k].minus(self.work) + self.r_diff[k].times(-1.0) + + def build_difference_matrices_for_homotopy(self, mu=0.0): + """ + Constructs the solution and homotopy residual differences from the history + + Parameters + ---------- + mu : float + Homotopy continuation parameter + """ assert len(self.x_hist) == len(self.r_hist), \ 'AndersonMultiSecant() >> inconsistent list sizes!' # clear the previous x_diff and r_diff lists @@ -85,7 +120,7 @@ def build_difference_matrices(self, mu=0.0): self.r_diff[k].minus(self.work) self.r_diff[k].times(-1.0) - def solve(self, in_vec, out_vec, rel_tol=1e-15): + def solve(self, in_vec, out_vec, beta=1.0, rel_tol=1e-15): # store the difference matrices and rhs in numpy array format nvar = in_vec.get_num_var() dR = numpy.empty((nvar, len(self.x_diff))) @@ -96,9 +131,9 @@ def solve(self, in_vec, out_vec, rel_tol=1e-15): vec.get_base_data(dX[:,k]) rhs = numpy.empty((nvar)) in_vec.get_base_data(rhs) - dRinv = numpy.linalg.pinv(dR, rcond=1e-14) - sol = numpy.empty_like(rhs) - sol[:] = -rhs - numpy.matmul(dX - dR, numpy.matmul(dRinv,rhs)) + dRinv = numpy.linalg.pinv(dR, rcond=1e-3) + sol = numpy.zeros_like(rhs) + sol[:] = -beta*rhs - numpy.matmul(dX - beta*dR, numpy.matmul(dRinv,rhs)) out_vec.set_base_data(sol) # imports at the bottom to prevent circular import errors diff --git a/src/kona/linalg/matrices/hessian/basic.py b/src/kona/linalg/matrices/hessian/basic.py index 2711d57..41b49ed 100644 --- a/src/kona/linalg/matrices/hessian/basic.py +++ b/src/kona/linalg/matrices/hessian/basic.py @@ -166,14 +166,9 @@ def add_to_history(self, x_new, r_new): """ raise NotImplementedError # pragma: no cover - def build_difference_matrices(self, mu=0.0): + def build_difference_matrices(self): """ Constructs the solution and residual differences from the history - - Parameters - ---------- - mu : float - Homotopy continuation parameter """ raise NotImplementedError # pragma: no cover diff --git a/src/kona/linalg/vectors/common.py b/src/kona/linalg/vectors/common.py index 3a9d026..dc704f6 100644 --- a/src/kona/linalg/vectors/common.py +++ b/src/kona/linalg/vectors/common.py @@ -753,8 +753,11 @@ def equals_mangasarian(self, ineq, mult): assert isinstance(mult, DualVectorINEQ), \ 'DualVectorINEQ >> Invalid type for mult. ' + \ 'Must be DualVecINEQ!' - self.base.data[:] = -np.fabs(ineq.base.data[:] - mult.base.data[:])**3 \ - + ineq.base.data[:]**3 + mult.base.data[:]**3 + #self.base.data[:] = -np.fabs(ineq.base.data[:] - mult.base.data[:])**3 \ + # + ineq.base.data[:]**3 + mult.base.data[:]**3 + self.base.data[:] = -np.fabs(ineq.base.data[:] - mult.base.data[:]) \ + + ineq.base.data[:] + mult.base.data[:] + self.times(0.5) def deriv_mangasarian(self, ineq, mult): """ @@ -776,9 +779,12 @@ def deriv_mangasarian(self, ineq, mult): assert isinstance(mult, DualVectorINEQ), \ 'DualVectorINEQ >> Invalid type for mult. ' + \ 'Must be DualVecINEQ!' - self.base.data[:] = np.multiply(-3*np.sign(ineq.base.data[:] - mult.base.data[:]), - (ineq.base.data[:] - mult.base.data[:])**2) \ - + 3*ineq.base.data[:]**2 + #self.base.data[:] = np.multiply(-3*np.sign(ineq.base.data[:] - mult.base.data[:]), + # (ineq.base.data[:] - mult.base.data[:])**2) \ + # + 3*ineq.base.data[:]**2 + self.base.data[:] = -np.sign(ineq.base.data[:] - mult.base.data[:]) \ + + np.ones_like(ineq.base.data) + self.times(0.5) # package imports at the bottom to prevent import errors import numpy as np diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index 05234db..648403a 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -345,6 +345,51 @@ def get_optimality_and_feasiblity(self): feas = np.sqrt(feas) return opt, feas + def equals_primaldual_residual(self, dLdx, ineq_mult=None): + """ + Using dLdx=:math:`\\begin{pmatrix} \\nabla_x L && h && g \\end{pmatrix}`, which can be + obtained from the method equals_KKT_conditions (as well as the inequality multiplier when + inequality constraints are present) this method computes the following nonlinear vector + function: + + .. math:: + r(x,\\lambda_h,\\lambda_g;\\mu) = + \\begin{bmatrix} + \\left[\\nabla_x f(x, u) - \\nabla_x h(x, u)^T \\lambda_{h} - \\nabla_x g(x, u)^T \\lambda_{g}\\right] \\\\ + -h(x,u) \\\\ + -|g(x,u) - \\lambda_g| + g(x,u) + \\lambda_g + \\end{bmatrix} + + where :math:`h(x,u)` are the equality constraints, and :math:`g(x,u)` are the + inequality constraints. The vectors :math:`\\lambda_h` and :math:`\\lambda_g` + are the associated Lagrange multipliers. + + Parameters + ---------- + dLdx : PrimalDualVector + The total derivative of the Lagranginan with respect to the primal and dual variables. + ineq_mult : DualVectorINEQ, optional + The current multiplier associated with the inequality constraints. + """ + assert isinstance(dLdx, PrimalDualVector), \ + "PrimalDualVector() >> dLdx must be a PrimalDualVector!" + if dLdx.ineq is not None or ineq_mult is not None: + assert isinstance(ineq_mult, DualVectorINEQ), \ + "PrimalDualVector() >> ineq_mult, if present, must be a DualVectorINEQ!" + assert dLdx.ineq is not None and ineq_mult is not None, \ + "PrimalDualVector() >> dLdx.ineq and ineq_mult are inconsistent!" + if dLdx == self: + "PrimalDualVector() >> equals_primaldual_residual is not in-place!" + # construct the primal part of the residual: dLdx + self.primal.equals(dLdx.primal) + if dLdx.eq is not None: + # include the equality constraint part: -h + self.eq.equals(dLdx.eq) + self.eq.times(-1.) + if dLdx.ineq is not None: + # include the inequality constraint part + self.ineq.equals_mangasarian(dLdx.ineq, ineq_mult) + def equals_homotopy_residual(self, dLdx, x, init, mu=1.0): """ Using dLdx=:math:`\\begin{pmatrix} \\nabla_x L && h && g \\end{pmatrix}`, which can be @@ -379,6 +424,8 @@ def equals_homotopy_residual(self, dLdx, x, init, mu=1.0): """ assert isinstance(dLdx, PrimalDualVector), \ "PrimalDualVector() >> dLdx must be a PrimalDualVector!" + assert isinstance(x, PrimalDualVector), \ + "PrimalDualVector() >> x must be a PrimalDualVector!" assert isinstance(init, PrimalDualVector), \ "PrimalDualVector() >> init must be a PrimalDualVector!" if dLdx.eq is None or self.eq is None or x.eq is None or init.eq is None: @@ -403,6 +450,20 @@ def equals_homotopy_residual(self, dLdx, x, init, mu=1.0): self.ineq.equals_mangasarian(self.ineq, x.ineq) self.ineq.minus((1.0 - mu)*0.1) + def equals_projgrad_residual(self, dLdx, x, init, mu=1.0): + # construct the primal part of the residual: mu*dLdx + (1-mu)(x - x_0) + self.primal.equals_ax_p_by(1.0, x.primal, -1.0, init.primal) + self.primal.equals_ax_p_by(mu, dLdx.primal, (1.0 - mu), self.primal) + if dLdx.eq is not None: + # include the equality constraint part: -mu*h - (1-mu)*lambda_h + self.eq.equals_ax_p_by(-mu, dLdx.eq, -(1.0 - mu), x.eq) + #self.eq.equals_ax_p_by(mu, dLdx.eq, -(1.0 - mu), x.eq) + if dLdx.ineq is not None: + # include the inequality constraint part + self.ineq.equals_ax_p_by(1.0, dLdx.ineq, -(1.0 - mu), init.ineq) + self.ineq.equals_mangasarian(self.ineq, x.ineq) + self.ineq.minus((1.0 - mu)*0.1) + def equals_predictor_rhs(self, dLdx, x, init, mu=1.0): """ Using dLdx=:math:`\\begin{pmatrix} \\nabla_x L && h && g \\end{pmatrix}`, which can be diff --git a/src/kona/test/test_dual_vectors.py b/src/kona/test/test_dual_vectors.py index 65ad066..b07c706 100644 --- a/src/kona/test/test_dual_vectors.py +++ b/src/kona/test/test_dual_vectors.py @@ -61,11 +61,13 @@ def test_equals_mangasarian(self): self.mult.equals(-9) self.mang.equals_mangasarian(self.dv, self.mult) for i in xrange(len(self.mang.base.data)): - self.assertEqual(self.mang.base.data[i], -1730) + self.assertEqual(self.mang.base.data[i], -10) # linear Mangasarian + # self.assertEqual(self.mang.base.data[i], -1730) # cubic Mangasarian self.mult.equals(-11) self.mang.equals_mangasarian(self.dv, self.mult) for i in xrange(len(self.mang.base.data)): - self.assertEqual(self.mang.base.data[i], -2332) + self.assertEqual(self.mang.base.data[i], -11) # linear Mangasarian + # self.assertEqual(self.mang.base.data[i], -2332) # cubic Mangasarian class DualVectorEQTestCaseIDF(unittest.TestCase): diff --git a/src/kona/test/test_multi_secant.py b/src/kona/test/test_multi_secant.py index 18b8d4d..101ed33 100644 --- a/src/kona/test/test_multi_secant.py +++ b/src/kona/test/test_multi_secant.py @@ -367,6 +367,12 @@ def test_build_difference_matrices_case1(self): aa.x_diff[k].get_base_data(A[:]) for i in range(A.shape[0]): self.assertAlmostEqual(A[i], 1.0, places=14) + aa.build_difference_matrices_for_homotopy() + A = np.zeros(3) + for k in range(max_stored-1): + aa.x_diff[k].get_base_data(A[:]) + for i in range(A.shape[0]): + self.assertAlmostEqual(A[i], 1.0, places=14) def test_build_difference_matrices_case2(self): # case 2: constrained with equality constraints only @@ -404,6 +410,12 @@ def test_build_difference_matrices_case2(self): aa.x_diff[k].get_base_data(A[:]) for i in range(A.shape[0]): self.assertAlmostEqual(A[i], 1.0, places=14) + aa.build_difference_matrices_for_homotopy() + A = np.zeros(13) + for k in range(max_stored-1): + aa.x_diff[k].get_base_data(A[:]) + for i in range(A.shape[0]): + self.assertAlmostEqual(A[i], 1.0, places=14) def test_build_difference_matrices_case3(self): # case 3: constrained with inequality constraints only @@ -441,6 +453,12 @@ def test_build_difference_matrices_case3(self): aa.x_diff[k].get_base_data(A[:]) for i in range(A.shape[0]): self.assertAlmostEqual(A[i], 1.0, places=14) + aa.build_difference_matrices_for_homotopy() + A = np.zeros(15) + for k in range(max_stored-1): + aa.x_diff[k].get_base_data(A[:]) + for i in range(A.shape[0]): + self.assertAlmostEqual(A[i], 1.0, places=14) def test_build_difference_matrices_case4(self): # case 4: constrained with both equality and inequality constraints @@ -482,6 +500,12 @@ def test_build_difference_matrices_case4(self): aa.x_diff[k].get_base_data(A[:]) for i in range(A.shape[0]): self.assertAlmostEqual(A[i], 1.0, places=14) + aa.build_difference_matrices_for_homotopy() + A = np.zeros(18) + for k in range(max_stored-1): + aa.x_diff[k].get_base_data(A[:]) + for i in range(A.shape[0]): + self.assertAlmostEqual(A[i], 1.0, places=14) def test_solve_case1(self): # case 1: unconstrained @@ -517,7 +541,7 @@ def test_solve_case1(self): x.set_base_data(np.array([0.,0.,0.])) r.set_base_data(np.array([0.,0.,0.])) aa.add_to_history(x, r) - aa.build_difference_matrices(mu=1.0) + aa.build_difference_matrices_for_homotopy(mu=1.0) # AA should be able to recover exact inverse x.set_base_data(np.array([1.,0.,0.])) diff --git a/src/kona/test/test_primal_dual_vector.py b/src/kona/test/test_primal_dual_vector.py index fe364ab..a4339db 100644 --- a/src/kona/test/test_primal_dual_vector.py +++ b/src/kona/test/test_primal_dual_vector.py @@ -278,6 +278,105 @@ def test_kkt_conditions_case4(self): exp_dLdX_norm = np.sqrt(10. * 0.**2) self.assertEqual(self.pd_vec4.primal.norm2, exp_dLdX_norm) + def test_equals_primaldual_residual_case1(self): + # case 1 has both equality and inequality constraints + # recall: self.dual = 1, so dCeqdU^T*dual.eq + dCineqdU^Tdual.ineq = 5 + 5 = 10 + dCdU(self.design, self.state).T.product(CompositeDualVector(self.dual_eq, self.dual_ineq), + self.adjoint, self.state_work) + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.minus(self.adjoint) # L = f - ceq^T*lambda_eq - cineq^T*lambda_ineq + self.state_work.times(-1.) + # We get adjoint = -11 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + dLdx = self.pd_work11 + dLdx.equals_KKT_conditions(self.at_pd1, self.state, self.adjoint) + + init = self.pd_work12 + init.equals_init_guess() + init.eq.equals_constraints(init.primal, self.state) + init.ineq.equals_constraints(init.primal, self.state) + + ineq_mult = self.at_pd1.ineq + self.pd_vec1.equals_primaldual_residual(dLdx, ineq_mult) + + # check results + exp_dLdX_norm = np.sqrt(10. * 20.**2) + self.assertAlmostEqual(self.pd_vec1.primal.norm2, exp_dLdX_norm, places=10) + exp_dLdEq_norm = np.sqrt(5. * 200**2) + self.assertAlmostEqual(self.pd_vec1.eq.norm2, exp_dLdEq_norm, places=10) + # exp_dLdIn_norm = np.sqrt(5. * 29701.95**2) # for cubic Mangasarian + exp_dLdIn_norm = np.sqrt(5. * 1**2) # for linear Mangasarian + self.assertAlmostEqual(self.pd_vec1.ineq.norm2, exp_dLdIn_norm, places=10) + + def test_equals_primaldual_residual_case2(self): + # case 2 has inequality constraints only + # recall: self.dual = 1, so dCineqdU^Tdual.ineq = 5 + dCdU(self.design, self.state).T.product(self.dual_ineq, self.adjoint, self.state_work) + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.minus(self.adjoint) # L = f - cineq^T*lambda_ineq + self.state_work.times(-1.) + # We get adjoint = -6 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + dLdx = self.pd_work21 + dLdx.equals_KKT_conditions(self.at_pd2, self.state, self.adjoint) + init = self.pd_work22 + init.equals_init_guess() + init.ineq.equals_constraints(init.primal, self.state) + + ineq_mult = self.at_pd2.ineq + self.pd_vec2.equals_primaldual_residual(dLdx, ineq_mult) + + # check results + exp_dLdX_norm = np.sqrt(10. * 10**2) + self.assertAlmostEqual(self.pd_vec2.primal.norm2, exp_dLdX_norm, places=10) + # exp_dLdIn_norm = np.sqrt(5. * 29701.95**2) # for cubic Mangasarian + exp_dLdIn_norm = np.sqrt(5. * 1**2) # for linear Mangasarian + self.assertAlmostEqual(self.pd_vec2.ineq.norm2, exp_dLdIn_norm, places=10) + + def test_equals_primaldual_residual_case3(self): + # case 3 has equality constraints only + # recall: self.dual = 1, so dCeqdU^Tdual.ineq = 5 + dCdU(self.design, self.state).T.product(self.dual_eq, self.adjoint, self.state_work) + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.minus(self.adjoint) # L = f - ceq^T*lambda_eq + self.state_work.times(-1.) + # We get adjoint = -6 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + dLdx = self.pd_work31 + dLdx.equals_KKT_conditions(self.at_pd3, self.state, self.adjoint) + init = self.pd_work32 + init.equals_init_guess() + init.eq.equals_constraints(init.primal, self.state) + + self.pd_vec3.equals_primaldual_residual(dLdx) + + # check results + exp_dLdX_norm = np.sqrt(10. * 10**2) + self.assertAlmostEqual(self.pd_vec3.primal.norm2, exp_dLdX_norm, places=10) + exp_dLdEq_norm = np.sqrt(5. * 200**2) + self.assertAlmostEqual(self.pd_vec3.eq.norm2, exp_dLdEq_norm, places=10) + + def test_equals_primaldual_residual_case4(self): + # case 4 has no constraints + # recall: dFdU = -1 + self.state_work.equals_objective_partial(self.design, self.state) + self.state_work.times(-1.) + # We get adjoint = -1 + dRdU(self.design, self.state).T.solve(self.state_work, self.adjoint) + dLdx = self.pd_work41 + dLdx.equals_KKT_conditions(self.at_pd4, self.state, self.adjoint) + init = self.pd_work42 + init.equals_init_guess() + + self.pd_vec4.equals_primaldual_residual(dLdx) + + # check results + exp_dLdX_norm = np.sqrt(10. * 0**2) + self.assertAlmostEqual(self.pd_vec3.primal.norm2, exp_dLdX_norm, places=10) + def test_homotopy_residual_case1(self): # case 1 has both equality and inequality constraints # recall: self.dual = 1, so dCeqdU^T*dual.eq + dCineqdU^Tdual.ineq = 5 + 5 = 10 @@ -305,7 +404,8 @@ def test_homotopy_residual_case1(self): self.assertAlmostEqual(self.pd_vec1.primal.norm2, exp_dLdX_norm, places=10) exp_dLdEq_norm = np.sqrt(5. * 100.5**2) self.assertAlmostEqual(self.pd_vec1.eq.norm2, exp_dLdEq_norm, places=10) - exp_dLdIn_norm = np.sqrt(5. * 29701.95**2) + # exp_dLdIn_norm = np.sqrt(5. * 29701.95**2) # for cubic Mangasarian + exp_dLdIn_norm = np.sqrt(5. * 0.95**2) # for linear Mangasarian self.assertAlmostEqual(self.pd_vec1.ineq.norm2, exp_dLdIn_norm, places=10) # test get_optimality_and_feasiblity while we are at it @@ -335,7 +435,8 @@ def test_homotopy_residual_case2(self): # check results exp_dLdX_norm = np.sqrt(10. * 5**2) self.assertAlmostEqual(self.pd_vec2.primal.norm2, exp_dLdX_norm, places=10) - exp_dLdIn_norm = np.sqrt(5. * 29701.95**2) + # exp_dLdIn_norm = np.sqrt(5. * 29701.95**2) # for cubic Mangasarian + exp_dLdIn_norm = np.sqrt(5. * 0.95**2) # for linear Mangasarian self.assertAlmostEqual(self.pd_vec2.ineq.norm2, exp_dLdIn_norm, places=10) # test get_optimality_and_feasiblity while we are at it From 150064e3ceb9a7ca2a88aabab79e4a326717c653 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Fri, 14 Jul 2017 09:22:32 -0400 Subject: [PATCH 10/18] beginning Schur-based preconditioner for multi-secant algorithm --- src/kona/linalg/matrices/preconds/schur.py | 70 ++++++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100644 src/kona/linalg/matrices/preconds/schur.py diff --git a/src/kona/linalg/matrices/preconds/schur.py b/src/kona/linalg/matrices/preconds/schur.py new file mode 100644 index 0000000..a9c4460 --- /dev/null +++ b/src/kona/linalg/matrices/preconds/schur.py @@ -0,0 +1,70 @@ +import numpy as np +from kona.linalg.matrices.hessian.basic import BaseHessian + +class ApproxSchur(BaseHessian): + """ + Uses an approximate SVD to approximate the Schur complement of the primal-dual matrix. + + Parameters + ---------- + + + Attributes + ---------- + subspace_size : int + Number of iterations in the Lanczos algorithm, and the number of singular values + approximated by the decomposition. + """ + def __init__(self, vector_factories, optns=None): + super(ApproxSchur, self).__init__(vector_factories, optns) + + # request vectors needed (work) + dual_work = None + if self.eq_factory is not None: + self.eq_factory.request_num_vectors(1) + if self.ineq_factory is not None: + self.ineq_factory.request_num_vectors(1) + + # initialize the total constraint jacobian block + self.dCdX = TotalConstraintJacobian(vector_factories) + + # Questions for Alp: + # 1) LowRankSVD needs a rev_factory for rectangular matrices; what if CompositeDualVector? + # 2) Would it be safe to implement a CompositeDualVector factory? + # 3) Looks like LowRankSVD assumes BaseVectors, correct? + + # initialize the low-rank SVD + # Will these definitions persist!!!! + def fwd_mat_vec(in_vec, out_vec): + self.dCdX.product(in_vec, out_vec) + def rev_mat_vec(in_vec, out_vec): + self.dCdX.T.product(in_vec, out_vec) + self.Schur_inv = LowRankSVD(fwd_mat_vec, self.primal_factory, + rev_mat_vec, [self.eq_factory, self.ineq_factory]) + + def linearize(self, at_primal, at_state, scale=1.0): + # need to generate dual_work here + + def product(self, in_vec, out_vec): + assert isinstance(in_vec, PrimalDualVector), \ + "ApproxSchur() >> in_vec must be a PrimalDualVector!" + assert isinstance(out_vec, PrimalDualVector), \ + "ApproxSchur() >> out_vec must be a PrimalDualVector!" + # set some aliases + dual_work = self.dual_work + in_primal = in_vec.primal + in_dual = in_vec.get_dual() + out_primal = out_vec.primal + out_dual = out_vec.get_dual() + if dual_work is not None: + # Step 1: w <-- A*u_p + u_d + self.dCdX.product(in_primal, dual_work) + dual_work.plus(in_dual) + # Step 2: v_d <-- (A*A^T)^-1*w + self.Schur_inv.inv_schur_product(dual_work, out_dual) + # Step 3: v_p <-- u_p - A^T*v_d + self.dCdX.T.product(out_dual, out_primal) + else: + # no preconditioning for unconstrained case + out_primal.equals(0.) + out_primal.equals_ax_p_by(1., in_primal, -1., out_primal) From 30fbfab96738c87af8836b391042082fad1edee0 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Wed, 19 Jul 2017 13:31:51 -0400 Subject: [PATCH 11/18] added an approximate inverse for A*A^T based on low-rank SVD --- .../linalg/matrices/preconds/low_rank_svd.py | 16 +- src/kona/test/test_low_rank_svd.py | 147 +++++++++++++++++- 2 files changed, 160 insertions(+), 3 deletions(-) diff --git a/src/kona/linalg/matrices/preconds/low_rank_svd.py b/src/kona/linalg/matrices/preconds/low_rank_svd.py index d52a5cb..accbe99 100644 --- a/src/kona/linalg/matrices/preconds/low_rank_svd.py +++ b/src/kona/linalg/matrices/preconds/low_rank_svd.py @@ -56,8 +56,8 @@ def __init__(self, fwd_mat_vec, fwd_factory, self._allocated = False # request vector memory for future allocation - self.fwd_factory.request_num_vectors(2*self.subspace_size + 1) - self.rev_factory.request_num_vectors(2*self.subspace_size) + self.fwd_factory.request_num_vectors(2*self.subspace_size + 2) + self.rev_factory.request_num_vectors(2*self.subspace_size + 1) def linearize(self): if not self._allocated: @@ -129,3 +129,15 @@ def approx_rev_prod(self, in_vec, out_vec): out_vec.equals(0.0) for i in xrange(len(self.V)): out_vec.equals_ax_p_by(1., out_vec, SUT_vec[i], self.V[i]) + + def inv_schur_prod(self, in_vec, out_vec): + UT_vec = np.zeros(len(self.U)-1) + for i in xrange(len(self.U)-1): + UT_vec[i] = self.U[i].inner(in_vec) + S2invUT_vec = np.zeros(len(self.U)-1) + for i in xrange(len(self.U)-1): + S2invUT_vec[i] = UT_vec[i]*(1./self.S[i,i]**2 - 1./self.S[-1,-1]**2) + out_vec.equals(in_vec) + out_vec.times(1./self.S[-1,-1]**2) + for i in xrange(len(self.U)-1): + out_vec.equals_ax_p_by(1.0, out_vec, S2invUT_vec[i], self.U[i]) diff --git a/src/kona/test/test_low_rank_svd.py b/src/kona/test/test_low_rank_svd.py index a8ef3b8..4855db1 100644 --- a/src/kona/test/test_low_rank_svd.py +++ b/src/kona/test/test_low_rank_svd.py @@ -1,13 +1,18 @@ import unittest +import inspect, os +import numpy as np + from kona.examples import Sellar +from kona.user import UserSolver from kona.linalg.memory import KonaMemory from kona.linalg.matrices.common import dCdU, dRdU from kona.linalg.matrices.hessian import TotalConstraintJacobian from kona.linalg.matrices.hessian import LagrangianHessian from kona.linalg.matrices.preconds import LowRankSVD from kona.linalg.vectors.composite import ReducedKKTVector -from kona.linalg.vectors.composite import CompositePrimalVector +from kona.linalg.vectors.composite import CompositeFactory, CompositePrimalVector, \ + CompositeDualVector class LowRankSVDTestCase(unittest.TestCase): @@ -148,5 +153,145 @@ def mat_vec(in_vec, out_vec): self.assertTrue(rel_error <= 1e-8) + def test_inv_schur_prod_eq(self): + '''LowRankSVD approximation of inverse of Schur complement (equality only)''' + max_lanczos = 10 + solver = UserSolver(20, num_eq=10) + km = KonaMemory(solver) + pf = km.primal_factory + eqf = km.eq_factory + + path = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) + A = np.loadtxt(path+'/synthetic_jac.dat') + Sinv = np.loadtxt(path+'/synthetic_schur.dat') + + # define the Low-rank SVD + def fwd_mat_vec(in_vec, out_vec): + out_vec.base.data[:] = A.dot(in_vec.base.data) + + def rev_mat_vec(in_vec, out_vec): + out_vec.base.data[:] = A.T.dot(in_vec.base.data) + + svd_optns = {'lanczos_size': max_lanczos} + self.svd = LowRankSVD( + fwd_mat_vec, pf, rev_mat_vec, eqf, svd_optns) + + # add request that will be used to create in_vec and out_vec + km.eq_factory.request_num_vectors(2) + km.allocate_memory() + + in_vec = km.eq_factory.generate() + out_vec = km.eq_factory.generate() + + # loop over and check each column in the approximate Schur + self.svd.linearize() + for i in xrange(km.neq): + in_vec.base.data = np.zeros_like(in_vec.base.data) + in_vec.base.data[i] = 1.0 + out_vec.base.data = np.zeros_like(out_vec.base.data) + self.svd.inv_schur_prod(in_vec, out_vec) + for j in xrange(km.neq): + self.assertAlmostEqual(Sinv[j,i]/out_vec.base.data[j], 1.0, places=9) + + def test_inv_schur_prod_ineq(self): + '''LowRankSVD approximation of inverse of Schur complement (inequality only)''' + max_lanczos = 10 + solver = UserSolver(20, num_ineq=10) + km = KonaMemory(solver) + pf = km.primal_factory + ineqf = km.ineq_factory + + path = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) + A = np.loadtxt(path+'/synthetic_jac.dat') + Sinv = np.loadtxt(path+'/synthetic_schur.dat') + + # define the Low-rank SVD + def fwd_mat_vec(in_vec, out_vec): + out_vec.base.data[:] = A.dot(in_vec.base.data) + + def rev_mat_vec(in_vec, out_vec): + out_vec.base.data[:] = A.T.dot(in_vec.base.data) + + svd_optns = {'lanczos_size': max_lanczos} + self.svd = LowRankSVD( + fwd_mat_vec, pf, rev_mat_vec, ineqf, svd_optns) + + # add request that will be used to create in_vec and out_vec + km.ineq_factory.request_num_vectors(2) + km.allocate_memory() + + in_vec = km.ineq_factory.generate() + out_vec = km.ineq_factory.generate() + + # loop over and check each column in the approximate Schur + self.svd.linearize() + for i in xrange(km.nineq): + in_vec.base.data = np.zeros_like(in_vec.base.data) + in_vec.base.data[i] = 1.0 + out_vec.base.data = np.zeros_like(out_vec.base.data) + self.svd.inv_schur_prod(in_vec, out_vec) + for j in xrange(km.nineq): + self.assertAlmostEqual(Sinv[j,i]/out_vec.base.data[j], 1.0, places=9) + + def test_inv_schur_prod_comp(self): + '''LowRankSVD approximation of inverse of Schur complement (CompositeDualVector)''' + max_lanczos = 10 + solver = UserSolver(20, num_eq=6, num_ineq=4) + km = KonaMemory(solver) + pf = km.primal_factory + eqf = km.eq_factory + ineqf = km.ineq_factory + dualf = CompositeFactory(km, CompositeDualVector) + + path = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) + A = np.loadtxt(path+'/synthetic_jac.dat') + Sinv = np.loadtxt(path+'/synthetic_schur.dat') + + # define the Low-rank SVD + def fwd_mat_vec(in_vec, out_vec): + out_copy = np.zeros(10) + out_copy[:] = A.dot(in_vec.base.data) + out_vec.eq.base.data[:] = out_copy[0:6] + out_vec.ineq.base.data[:] = out_copy[6:] + + def rev_mat_vec(in_vec, out_vec): + in_copy = np.zeros(10) + in_copy[0:6] = in_vec.eq.base.data[:] + in_copy[6:] = in_vec.ineq.base.data[:] + out_vec.base.data[:] = A.T.dot(in_copy) + + svd_optns = {'lanczos_size': max_lanczos} + self.svd = LowRankSVD( + fwd_mat_vec, pf, rev_mat_vec, dualf, svd_optns) + + # add request that will be used to create in_vec and out_vec + dualf.request_num_vectors(2) + km.allocate_memory() + + in_vec = dualf.generate() + out_vec = dualf.generate() + + # loop over and check each column in the approximate Schur + self.svd.linearize() + for i in xrange(km.neq): + in_vec.equals(0.) + in_vec.eq.base.data[i] = 1.0 + out_vec.equals(0.) + self.svd.inv_schur_prod(in_vec, out_vec) + for j in xrange(km.neq): + self.assertAlmostEqual(Sinv[j,i]/out_vec.eq.base.data[j], 1.0, places=9) + for j in xrange(km.nineq): + self.assertAlmostEqual(Sinv[km.neq+j,i]/out_vec.ineq.base.data[j], 1.0, places=9) + for i in xrange(km.nineq): + in_vec.equals(0.) + in_vec.ineq.base.data[i] = 1.0 + out_vec.equals(0.) + self.svd.inv_schur_prod(in_vec, out_vec) + for j in xrange(km.neq): + self.assertAlmostEqual(Sinv[j,km.neq+i]/out_vec.eq.base.data[j], 1.0, places=9) + for j in xrange(km.nineq): + self.assertAlmostEqual(Sinv[km.neq+j,km.neq+i]/out_vec.ineq.base.data[j], 1.0, + places=9) + if __name__ == "__main__": unittest.main() From 1d02d8544477e2212ee0e4c8e9b56b32fb580d5e Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Wed, 19 Jul 2017 22:04:31 -0400 Subject: [PATCH 12/18] finished and tested draft version of ApproxSchur preconditioner --- src/kona/linalg/matrices/preconds/schur.py | 91 +++++++++----- src/kona/optimizer.py | 2 +- src/kona/test/test_approx_schur.py | 132 +++++++++++++++++++++ 3 files changed, 196 insertions(+), 29 deletions(-) create mode 100644 src/kona/test/test_approx_schur.py diff --git a/src/kona/linalg/matrices/preconds/schur.py b/src/kona/linalg/matrices/preconds/schur.py index a9c4460..0d6d466 100644 --- a/src/kona/linalg/matrices/preconds/schur.py +++ b/src/kona/linalg/matrices/preconds/schur.py @@ -1,49 +1,82 @@ import numpy as np +from kona.linalg.vectors.composite import CompositeFactory, PrimalDualVector, CompositeDualVector from kona.linalg.matrices.hessian.basic import BaseHessian +from kona.linalg.matrices.hessian import TotalConstraintJacobian +from kona.linalg.matrices.preconds.low_rank_svd import LowRankSVD class ApproxSchur(BaseHessian): """ - Uses an approximate SVD to approximate the Schur complement of the primal-dual matrix. + Uses an approximate SVD to estimate the inverse Schur complement of the primal-dual matrix. + Specifically, this preconditioner forms an approximate inverse of the augmented matrix - Parameters - ---------- + .. math:: + \\begin{bmatrix} I && -A^T \\\\ -A && 0 \\end{bmatrix} + + where :math:`A` is the total constraint Jacobian. The preconditioner is based on an LU + decomposition with the Schur complement approximated by a low-rank update to a diagonal matrix. + The action of the preconditioner on a vector :math:`[u_p^T u_v^T]^T` is given by + + .. math:: + \\begin{bmatrix} u_p - \\tilde{A}^T (\\tilde{A}\\tilde{A}^T)^{-1}(\\tilde{A}u_p - u_d) \\\\ + - (\\tilde{A}\\tilde{A}^T)^{-1}(\\tilde{A}u_p - u_d) \\end{bmatrix} + where :math:`\\tilde{A}` is the low-rank approximation to the constraint Jacobian, and + :math:`(\\tilde{A}\\tilde{A}^{T})^{-1}` is found using a pseudo-inverse-like approximation: see + LowRankSVD.inv_schur_prod for further details on this. Attributes ---------- - subspace_size : int - Number of iterations in the Lanczos algorithm, and the number of singular values - approximated by the decomposition. + dCdX : TotalConstraintJacobian + Used in Lanczos bi-diagonalization to find low-rank SVD + dCdX_approx : LowRankSVD + Used to approximate the inverse of (dCdX)*(dCdX^T) """ - def __init__(self, vector_factories, optns=None): + def __init__(self, vector_factories, optns=None): super(ApproxSchur, self).__init__(vector_factories, optns) - # request vectors needed (work) - dual_work = None - if self.eq_factory is not None: - self.eq_factory.request_num_vectors(1) - if self.ineq_factory is not None: - self.ineq_factory.request_num_vectors(1) + # set-up the dual factory + self.dual_factory = None + if self.eq_factory is not None and self.ineq_factory is not None: + self.dual_factory = CompositeFactory(self.eq_factory._memory, CompositeDualVector) + elif self.eq_factory is not None: + self.dual_factory = self.eq_factory + else: + self.dual_factory = self.ineq_factory + + # request vectors needed (dual_work) + self.dual_work = None + if self.dual_factory is not None: + self.dual_factory.request_num_vectors(1) # initialize the total constraint jacobian block self.dCdX = TotalConstraintJacobian(vector_factories) - # Questions for Alp: - # 1) LowRankSVD needs a rev_factory for rectangular matrices; what if CompositeDualVector? - # 2) Would it be safe to implement a CompositeDualVector factory? - # 3) Looks like LowRankSVD assumes BaseVectors, correct? - # initialize the low-rank SVD - # Will these definitions persist!!!! def fwd_mat_vec(in_vec, out_vec): self.dCdX.product(in_vec, out_vec) + def rev_mat_vec(in_vec, out_vec): self.dCdX.T.product(in_vec, out_vec) - self.Schur_inv = LowRankSVD(fwd_mat_vec, self.primal_factory, - rev_mat_vec, [self.eq_factory, self.ineq_factory]) + + self.dCdX_approx = LowRankSVD(fwd_mat_vec, self.primal_factory, + rev_mat_vec, self.dual_factory, optns) + + self._allocated = False def linearize(self, at_primal, at_state, scale=1.0): - # need to generate dual_work here + # store references to the evaluation point, and save the scaling of constraint terms + self.at_design = at_primal + self.at_state = at_state + self.scale = scale + + # get the total constraint Jacobian ready, and then use it in Lanczos bi-diagonalization + self.dCdX.linearize(self.at_design, self.at_state, scale=self.scale) + self.dCdX_approx.linearize() + + # check if dual_work needs to be allocated + if not self._allocated: + if self.dual_factory is not None: + self.dual_work = self.dual_factory.generate() def product(self, in_vec, out_vec): assert isinstance(in_vec, PrimalDualVector), \ @@ -58,13 +91,15 @@ def product(self, in_vec, out_vec): out_dual = out_vec.get_dual() if dual_work is not None: # Step 1: w <-- A*u_p + u_d - self.dCdX.product(in_primal, dual_work) + self.dCdX_approx.approx_fwd_prod(in_primal, dual_work) dual_work.plus(in_dual) - # Step 2: v_d <-- (A*A^T)^-1*w - self.Schur_inv.inv_schur_product(dual_work, out_dual) - # Step 3: v_p <-- u_p - A^T*v_d - self.dCdX.T.product(out_dual, out_primal) + # Step 2: v_d <-- -(A*A^T)^-1*w + self.dCdX_approx.inv_schur_prod(dual_work, out_dual) + out_dual.times(-1.) + # Step 3: v_p <-- A^T*v_d + self.dCdX_approx.approx_rev_prod(out_dual, out_primal) else: # no preconditioning for unconstrained case out_primal.equals(0.) - out_primal.equals_ax_p_by(1., in_primal, -1., out_primal) + # v_p <-- v_p + u_p + out_primal.equals_ax_p_by(1., in_primal, 1., out_primal) diff --git a/src/kona/optimizer.py b/src/kona/optimizer.py index fd37006..4d39c25 100644 --- a/src/kona/optimizer.py +++ b/src/kona/optimizer.py @@ -119,4 +119,4 @@ def solve(self, print_opts=False): import numpy as np from kona.options import print_dict from kona.user import UserSolver -from kona.linalg.memory import KonaMemory \ No newline at end of file +from kona.linalg.memory import KonaMemory diff --git a/src/kona/test/test_approx_schur.py b/src/kona/test/test_approx_schur.py new file mode 100644 index 0000000..d71e95a --- /dev/null +++ b/src/kona/test/test_approx_schur.py @@ -0,0 +1,132 @@ +import numpy as np +import unittest +import inspect +import os + +from kona import Optimizer +from kona.algorithms.base_algorithm import OptimizationAlgorithm +from kona.linalg.matrices.preconds.schur import ApproxSchur +from kona.user import UserSolver +from kona.linalg.vectors.composite import PrimalDualVector +# from kona.linalg.matrices.hessian.constraint_jacobian import TotalConstraintJacobian +from kona.linalg.matrices.hessian import TotalConstraintJacobian + +class SolverForApproxSchur(UserSolver): + """A minimal class used to provided synthetic total-constraint Jacobians""" + + def __init__(self, ndv=20, neq=10, nineq=0): + super(SolverForApproxSchur, self).__init__( + num_design=ndv, + num_state=0, + num_eq=neq, + num_ineq=nineq) + self.num_design = ndv + self.num_eq = neq + self.num_ineq = nineq + path = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) + self.dCdX = np.loadtxt(path+'/synthetic_jac.dat') + + def multiply_dCEQdX(self, at_design, at_state, in_vec): + if self.num_eq > 0: + return np.dot(self.dCdX[0:self.num_eq,:], in_vec) + else: + return 0.0 + + def multiply_dCEQdX_T(self, at_design, at_state, in_vec): + if self.num_eq > 0: + return np.dot(in_vec.T, self.dCdX[0:self.num_eq,:]).T + else: + return 0.0 + +class AlgorithmForApproxSchur(OptimizationAlgorithm,unittest.TestCase): + + def __init__(self, primal_factory, state_factory, + eq_factory, ineq_factory, optns=None): + # trigger base class initialization + super(AlgorithmForApproxSchur, self).__init__( + primal_factory, state_factory, eq_factory, ineq_factory, optns + ) + # number of vectors required in solve() method + self.primal_factory.request_num_vectors(3) + self.state_factory.request_num_vectors(1) + if self.eq_factory is not None: + self.eq_factory.request_num_vectors(2) + if self.ineq_factory is not None: + self.eq_factory.request_num_vectors(2) + max_lanczos = 10 + optns = {'lanczos_size': max_lanczos} + self.precond = ApproxSchur([primal_factory, state_factory, eq_factory, ineq_factory], optns) + km = self.primal_factory + self.num_design = km._memory.ndv + self.num_eq = km._memory.neq + self.num_ineq = km._memory.nineq + # algorithm also needs direct access to the Jacobian + self.dCdX = TotalConstraintJacobian([primal_factory, state_factory, eq_factory, + ineq_factory]) + + def _generate_vector(self): + dual_eq = None + dual_ineq = None + primal = self.primal_factory.generate() + if self.eq_factory is not None: + dual_eq = self.eq_factory.generate() + if self.ineq_factory is not None: + dual_ineq = self.ineq_factory.generate() + return PrimalDualVector(primal, eq_vec=dual_eq, ineq_vec=dual_ineq) + + def exact_product(self, in_vec, out_vec): + """Apply the exact augmented matrix using constraint Jacobian from file""" + self.dCdX.T.product(in_vec.get_dual(), out_vec.primal) + out_vec.primal.equals_ax_p_by(1.0, in_vec.primal, -1.0, out_vec.primal) + self.dCdX.product(in_vec.primal, out_vec.get_dual()) + out_vec.get_dual().times(-1.) + + def solve(self): + # at_state and at_design are not used, but they need to be defined + at_state = self.state_factory.generate() + at_design = self.primal_factory.generate() + u_vec = self._generate_vector() + v_vec = self._generate_vector() + self.dCdX.linearize(at_design, at_state) + self.precond.linearize(at_design, at_state) + for i in xrange(20): + u_vec.equals(0.) + v_vec.equals(0.) + u_vec.primal.base.data[i] = 1.0 + self.exact_product(u_vec, v_vec) + u_vec.equals(0.) + self.precond.product(v_vec, u_vec) + for j in xrange(20): + if j == i: + self.assertAlmostEqual(u_vec.primal.base.data[j], 1., places=10) + else: + self.assertAlmostEqual(u_vec.primal.base.data[j], 0., places=10) + for j in xrange(10): + self.assertAlmostEqual(u_vec.eq.base.data[j], 0., places=10) + + for i in xrange(10): + u_vec.equals(0.) + v_vec.equals(0.) + u_vec.eq.base.data[i] = 1.0 + self.exact_product(u_vec, v_vec) + u_vec.equals(0.) + self.precond.product(v_vec, u_vec) + for j in xrange(20): + self.assertAlmostEqual(u_vec.primal.base.data[j], 0., places=10) + for j in xrange(10): + if j == i: + self.assertAlmostEqual(u_vec.eq.base.data[j], 1., places=10) + else: + self.assertAlmostEqual(u_vec.eq.base.data[j], 0., places=10) + +class ApproxSchurTestCase(unittest.TestCase): + + def test_ApproxSchur_eq(self): + '''test Approximate Schur precondition (equality only)''' + solver = SolverForApproxSchur(ndv=20, neq=10) + algorithm = AlgorithmForApproxSchur + optimizer = Optimizer(solver, algorithm) + optimizer.solve() + +if __name__ == "__main__": + unittest.main() From bc43b07b6ec7ceb9c64bf43ed541d4201c4b464a Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Fri, 21 Jul 2017 10:38:51 -0400 Subject: [PATCH 13/18] introduced regularization parameter into AndersonMultiSecant to deal with nonconvexity --- .../algorithms/reduced_space_multi_secant.py | 39 ++++++++++++++--- .../matrices/hessian/anderson_multisecant.py | 43 ++++++++++++++++--- src/kona/linalg/matrices/preconds/schur.py | 29 ++++++++----- src/kona/linalg/vectors/composite.py | 5 +++ 4 files changed, 92 insertions(+), 24 deletions(-) diff --git a/src/kona/algorithms/reduced_space_multi_secant.py b/src/kona/algorithms/reduced_space_multi_secant.py index 7020c88..34b5674 100644 --- a/src/kona/algorithms/reduced_space_multi_secant.py +++ b/src/kona/algorithms/reduced_space_multi_secant.py @@ -1,6 +1,5 @@ from kona.algorithms.base_algorithm import OptimizationAlgorithm - class ReducedSpaceMultiSecant(OptimizationAlgorithm): """ Algorithm for generic nonlinearly constrained optiimzation problems. @@ -34,6 +33,22 @@ def __init__(self, primal_factory, state_factory, # iteration counter self.iter = 0 + # set the preconditioner for the multi-secant methods + self.precond = get_opt(self.optns, None, 'multi_secant', 'precond') + if self.precond is None: + # use the identity preconditioner + self.precond = IdentityMatrix() + elif self.precond is 'approx_schur': + # use the SVD-based approximate Schur preconditioner + precond_optns = get_opt(self.optns, {}, 'multi_secant') + self.precond = ApproxSchur([self.primal_factory, self.state_factory, self.eq_factory, + self.ineq_factory], precond_optns) + elif self.precond is 'idf_schur': + self.precond = ReducedSchurPreconditioner( + [primal_factory, state_factory, eq_factory, ineq_factory]) + else: + raise BadKonaOption(self.optns, 'precond') + # set the type of multi-secant method try: # the multisecant member is for the "second"-order scheme @@ -41,6 +56,7 @@ def __init__(self, primal_factory, state_factory, self.optns, AndersonMultiSecant, 'multi_secant', 'type') hessian_optns = get_opt(self.optns, {}, 'multi_secant') hessian_optns['out_file'] = self.info_file + self.hess_reg = get_opt(self.optns, 0.0, 'multi_secant', 'hess_reg') self.multisecant = multisecant([self.primal_factory, self.eq_factory, self.ineq_factory], hessian_optns) # the multisecant_safe member is for the first-order scheme @@ -50,6 +66,7 @@ def __init__(self, primal_factory, state_factory, except Exception: raise BadKonaOption(self.optns, 'multi_secant','type') + # set remaining options self.primal_tol_abs = get_opt(self.optns, 1e-6, 'opt_tol_abs') self.cnstr_tol_abs = get_opt(self.optns, 1e-6, 'feas_tol_abs') self.beta = get_opt(self.optns, 1.0, 'multi_secant', 'beta') @@ -130,6 +147,9 @@ def solve(self): self.multisecant.set_initial_data(X) self.multisecant_safe.set_initial_data(X) + # initialize the preconditioner + self.precond.linearize(X.primal, state) + # get the adjoint for -h^T\lambda_h - g^T\lambda_g adjoint.equals_homotopy_adjoint(X, state, state_work, obj_scale=0.0) R.equals_KKT_conditions(X, state, adjoint, obj_scale=0.0) @@ -186,14 +206,14 @@ def solve(self): break # get full multi-secant step and projected-gradient - self.multisecant.build_difference_matrices() - self.multisecant.solve(R, dX, self.beta) - self.multisecant_safe.build_difference_matrices() - self.multisecant_safe.solve(obj_grad, dX_safe) + self.multisecant.build_difference_matrices(self.hess_reg) + self.multisecant.solve(R, dX, self.beta, self.precond.product) + self.multisecant_safe.build_difference_matrices(self.hess_reg) + self.multisecant_safe.solve(obj_grad, dX_safe, precond=self.precond.product) info = ' ' - if dX.primal.inner(dX_safe.primal) < 0.: + if False: #dX.primal.inner(dX_safe.primal) < 0.: # the full step is considered unsafe: compute first-order step - self.multisecant_safe.solve(R, dX, self.beta) + self.multisecant_safe.solve(R, dX, self.beta, self.precond.product) info += 'safe-step' # remove history from multisecant #self.multisecant.clear_history() @@ -225,6 +245,8 @@ def solve(self): state.equals_primal_solution(X.primal) obj_val = objective_value(X.primal, state) + self.precond.linearize(X.primal, state) + # get the adjoint for -h^T\lambda_h - g^T\lambda_g adjoint.equals_homotopy_adjoint(X, state, state_work, obj_scale=0.0) R.equals_KKT_conditions(X, state, adjoint, obj_scale=0.0) @@ -258,4 +280,7 @@ def solve(self): from kona.linalg.common import current_solution, objective_value from kona.linalg.vectors.composite import PrimalDualVector from kona.linalg.matrices.hessian import AndersonMultiSecant +from kona.linalg.matrices.common import IdentityMatrix +from kona.linalg.matrices.preconds.schur import ApproxSchur +from kona.linalg.matrices.preconds import ReducedSchurPreconditioner from kona.linalg.solvers.util import EPS diff --git a/src/kona/linalg/matrices/hessian/anderson_multisecant.py b/src/kona/linalg/matrices/hessian/anderson_multisecant.py index 55de260..1d88a87 100644 --- a/src/kona/linalg/matrices/hessian/anderson_multisecant.py +++ b/src/kona/linalg/matrices/hessian/anderson_multisecant.py @@ -1,7 +1,8 @@ from collections import deque import numpy from kona.linalg.matrices.hessian.basic import MultiSecantApprox -from kona.linalg.vectors.composite import CompositeDualVector, PrimalDualVector +from kona.linalg.vectors.composite import PrimalDualVector +from kona.linalg.solvers.util import EPS class AndersonMultiSecant(MultiSecantApprox): """ @@ -74,12 +75,25 @@ def clear_history(self): self.x_hist = deque() self.r_hist = deque() - def build_difference_matrices(self): + def build_difference_matrices(self, alpha=0.0): + """ + Constructs the difference matrices based on first-order optimality + + Parameters + ---------- + alpha : Float + controls the amount of Hessian regularization + """ assert len(self.x_hist) == len(self.r_hist), \ 'AndersonMultiSecant() >> inconsistent list sizes!' + assert alpha >= 0.0, \ + 'AndersonMultiSecant() >> alpha must be non-negative in build_difference_matrices!' # clear the previous x_diff and r_diff lists del self.x_diff[:], self.r_diff[:] self.work.equals_primaldual_residual(self.r_hist[0], self.x_hist[0].ineq) + if alpha > EPS: + # include Hessian regularization + self.work.primal.equals_ax_p_by(1.0, self.work.primal, alpha, self.x_hist[0].primal) for k in range(len(self.x_hist)-1): # generate new vectors for x_diff and r_diff lists dx = self._generate_vector() @@ -90,6 +104,10 @@ def build_difference_matrices(self): self.x_diff[k].equals_ax_p_by(1.0, self.x_hist[k+1], -1.0, self.x_hist[k]) self.r_diff[k].equals(self.work) self.work.equals_primaldual_residual(self.r_hist[k+1], self.x_hist[k+1].ineq) + if alpha > EPS: + # include Hessian regularization + self.work.primal.equals_ax_p_by(1.0, self.work.primal, alpha, + self.x_hist[k+1].primal) self.r_diff[k].minus(self.work) self.r_diff[k].times(-1.0) @@ -120,7 +138,7 @@ def build_difference_matrices_for_homotopy(self, mu=0.0): self.r_diff[k].minus(self.work) self.r_diff[k].times(-1.0) - def solve(self, in_vec, out_vec, beta=1.0, rel_tol=1e-15): + def solve(self, in_vec, out_vec, beta=1.0, precond=None, rel_tol=1e-15): # store the difference matrices and rhs in numpy array format nvar = in_vec.get_num_var() dR = numpy.empty((nvar, len(self.x_diff))) @@ -131,9 +149,22 @@ def solve(self, in_vec, out_vec, beta=1.0, rel_tol=1e-15): vec.get_base_data(dX[:,k]) rhs = numpy.empty((nvar)) in_vec.get_base_data(rhs) - dRinv = numpy.linalg.pinv(dR, rcond=1e-3) + dRinv = numpy.linalg.pinv(dR, rcond=1e-6) sol = numpy.zeros_like(rhs) - sol[:] = -beta*rhs - numpy.matmul(dX - beta*dR, numpy.matmul(dRinv,rhs)) - out_vec.set_base_data(sol) + # sol[:] = -beta*rhs - numpy.matmul(dX - beta*dR, numpy.matmul(dRinv,rhs)) + # out_vec.set_base_data(sol) + dRinv_r = numpy.zeros_like(rhs) + dRinv_r = numpy.matmul(dRinv,rhs) + sol = numpy.matmul(dR, dRinv_r) - rhs + # move the base data into work so it can be preconditioned + self.work.set_base_data(sol) + if precond is None: + out_vec.equals(self.work) + else: + precond(self.work, out_vec) + out_vec.times(beta) + sol = numpy.matmul(dX, dRinv_r) + self.work.set_base_data(sol) + out_vec.equals_ax_p_by(1., out_vec, -1., self.work) # imports at the bottom to prevent circular import errors diff --git a/src/kona/linalg/matrices/preconds/schur.py b/src/kona/linalg/matrices/preconds/schur.py index 0d6d466..4bfadde 100644 --- a/src/kona/linalg/matrices/preconds/schur.py +++ b/src/kona/linalg/matrices/preconds/schur.py @@ -48,18 +48,23 @@ def __init__(self, vector_factories, optns=None): if self.dual_factory is not None: self.dual_factory.request_num_vectors(1) - # initialize the total constraint jacobian block - self.dCdX = TotalConstraintJacobian(vector_factories) + # initialize the total constraint jacobian block + self.dCdX = TotalConstraintJacobian(vector_factories) - # initialize the low-rank SVD - def fwd_mat_vec(in_vec, out_vec): - self.dCdX.product(in_vec, out_vec) + # initialize the low-rank SVD + def fwd_mat_vec(in_vec, out_vec): + self.dCdX.product(in_vec, out_vec) - def rev_mat_vec(in_vec, out_vec): - self.dCdX.T.product(in_vec, out_vec) + def rev_mat_vec(in_vec, out_vec): + self.dCdX.T.product(in_vec, out_vec) - self.dCdX_approx = LowRankSVD(fwd_mat_vec, self.primal_factory, - rev_mat_vec, self.dual_factory, optns) + self.dCdX_approx = LowRankSVD(fwd_mat_vec, self.primal_factory, + rev_mat_vec, self.dual_factory, optns) + + else: + # if self.dual_factory is None, there are no constraints + self.dCdX = None + self.dCdX_approx = None self._allocated = False @@ -70,13 +75,15 @@ def linearize(self, at_primal, at_state, scale=1.0): self.scale = scale # get the total constraint Jacobian ready, and then use it in Lanczos bi-diagonalization - self.dCdX.linearize(self.at_design, self.at_state, scale=self.scale) - self.dCdX_approx.linearize() + if self.dual_factory is not None: + self.dCdX.linearize(self.at_design, self.at_state, scale=self.scale) + self.dCdX_approx.linearize() # check if dual_work needs to be allocated if not self._allocated: if self.dual_factory is not None: self.dual_work = self.dual_factory.generate() + self._allocated = True def product(self, in_vec, out_vec): assert isinstance(in_vec, PrimalDualVector), \ diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index 9cb1b14..bf138f0 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -363,6 +363,11 @@ def get_dual(self): dual = self.ineq return dual + @property + def dual(self): + # added this property so we could use IDF preconditioner with PrimalDualVectors + return self.get_dual() + def equals_init_guess(self): """ Sets the primal-dual vector to the initial guess, using the initial design. From ab06c732c895c338d6d4311f0820a69d7d3adfa4 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Fri, 21 Jul 2017 11:31:27 -0400 Subject: [PATCH 14/18] removed multisecant_safe from ReducedSpaceMultiSecant --- .../algorithms/reduced_space_multi_secant.py | 73 +++---------------- src/kona/linalg/vectors/composite.py | 2 +- 2 files changed, 13 insertions(+), 62 deletions(-) diff --git a/src/kona/algorithms/reduced_space_multi_secant.py b/src/kona/algorithms/reduced_space_multi_secant.py index 34b5674..afefe2d 100644 --- a/src/kona/algorithms/reduced_space_multi_secant.py +++ b/src/kona/algorithms/reduced_space_multi_secant.py @@ -21,8 +21,8 @@ def __init__(self, primal_factory, state_factory, ) # number of vectors required in solve() method - # PrimalDualVectors = X, dX, R, dLdX, obj_grad, dX_safe - num_pd = 6 + # PrimalDualVectors = X, dX, R, dLdX + num_pd = 4 self.primal_factory.request_num_vectors(num_pd) if self.eq_factory is not None: self.eq_factory.request_num_vectors(num_pd) @@ -59,10 +59,6 @@ def __init__(self, primal_factory, state_factory, self.hess_reg = get_opt(self.optns, 0.0, 'multi_secant', 'hess_reg') self.multisecant = multisecant([self.primal_factory, self.eq_factory, self.ineq_factory], hessian_optns) - # the multisecant_safe member is for the first-order scheme - self.multisecant_safe = multisecant([self.primal_factory, self.eq_factory, - self.ineq_factory], hessian_optns) - except Exception: raise BadKonaOption(self.optns, 'multi_secant','type') @@ -128,10 +124,8 @@ def solve(self): # generate primal-dual vectors, and other vectors X = self._generate_vector() dX = self._generate_vector() - dX_safe = self._generate_vector() R = self._generate_vector() dLdX = self._generate_vector() - obj_grad = self._generate_vector() # generate state vectors state = self.state_factory.generate() @@ -143,25 +137,14 @@ def solve(self): state.equals_primal_solution(X.primal) obj_val = objective_value(X.primal, state) - # initialize the multi-secant methods + # initialize the multi-secant method and preconditioner self.multisecant.set_initial_data(X) - self.multisecant_safe.set_initial_data(X) - - # initialize the preconditioner self.precond.linearize(X.primal, state) - # get the adjoint for -h^T\lambda_h - g^T\lambda_g - adjoint.equals_homotopy_adjoint(X, state, state_work, obj_scale=0.0) - R.equals_KKT_conditions(X, state, adjoint, obj_scale=0.0) + # get the adjoint for dLdX, compute KKT conditions, and add to history + adjoint.equals_homotopy_adjoint(X, state, state_work) + R.equals_KKT_conditions(X, state, adjoint) dLdX.equals(R) - R.primal.plus(X.primal) - self.multisecant_safe.add_to_history(X, R) - - # get the adjoint for df/dx, then add it to finish dLdX computation - adjoint.equals_objective_adjoint(X.primal, state, state_work) - obj_grad.equals(0.) - obj_grad.primal.equals_total_gradient(X.primal, state, adjoint) - dLdX.primal.plus(obj_grad.primal) self.multisecant.add_to_history(X, dLdX) # send initial point info to the user @@ -208,56 +191,24 @@ def solve(self): # get full multi-secant step and projected-gradient self.multisecant.build_difference_matrices(self.hess_reg) self.multisecant.solve(R, dX, self.beta, self.precond.product) - self.multisecant_safe.build_difference_matrices(self.hess_reg) - self.multisecant_safe.solve(obj_grad, dX_safe, precond=self.precond.product) - info = ' ' - if False: #dX.primal.inner(dX_safe.primal) < 0.: - # the full step is considered unsafe: compute first-order step - self.multisecant_safe.solve(R, dX, self.beta, self.precond.product) - info += 'safe-step' - # remove history from multisecant - #self.multisecant.clear_history() # safe-guard against large steps - if dX.norm2 > self.radius_max: #/self.iter: - dX.times(self.radius_max/(dX.norm2)) #*self.iter)) + info = ' ' + if dX.norm2 > self.radius_max: + dX.times(self.radius_max/(dX.norm2)) info += ' length restricted' X.plus(dX) - # j = 0 - # while True: - # j += 1 - # X.plus(dX) - # # evaluate at new X, and check if optimality improved - # state.equals_primal_solution(X.primal) - # adjoint.equals_homotopy_adjoint(X, state, state_work) - # dLdX.equals_KKT_conditions(X, state, adjoint) - # R.equals_primaldual_residual(dLdX, X.ineq) - # self.info_file.write( - # 'line search iteration %i: |R(x_k)| = %e: |R(x_k + p_k)| = %e\n'%( - # j, grad_norm**2 + feas_norm**2, R.norm2**2)) - # if R.norm2**2 < grad_norm**2 + feas_norm**2 or j >= 4: - # break - # X.minus(dX) - # dX.times(0.1) - # evaluate at new X, construct first-order optimality conditions, and store state.equals_primal_solution(X.primal) obj_val = objective_value(X.primal, state) self.precond.linearize(X.primal, state) - # get the adjoint for -h^T\lambda_h - g^T\lambda_g - adjoint.equals_homotopy_adjoint(X, state, state_work, obj_scale=0.0) - R.equals_KKT_conditions(X, state, adjoint, obj_scale=0.0) + # get the adjoint for dLdX + adjoint.equals_homotopy_adjoint(X, state, state_work) + R.equals_KKT_conditions(X, state, adjoint) dLdX.equals(R) - R.primal.plus(X.primal) - self.multisecant_safe.add_to_history(X, R) - - # get the adjoint for df/dx, then add it to finish dLdX computation - adjoint.equals_objective_adjoint(X.primal, state, state_work) - obj_grad.primal.equals_total_gradient(X.primal, state, adjoint) - dLdX.primal.plus(obj_grad.primal) self.multisecant.add_to_history(X, dLdX) # output current solution info to the user diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index bf138f0..ed40237 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -381,7 +381,7 @@ def equals_init_guess(self): def equals_KKT_conditions(self, x, state, adjoint, obj_scale=1.0, cnstr_scale=1.0): """ Calculates the total derivative of the Lagrangian - :math:`\\mathcal{L}(x, u) = f(x, u)+ \\lambda_{h}^T h(x, u) + \\lambda_{g}^T g(x, u)` with + :math:`\\mathcal{L}(x, u) = f(x, u) - \\lambda_{h}^T h(x, u) - \\lambda_{g}^T g(x, u)` with respect to :math:`\\begin{pmatrix}x && \\lambda_{h} && \\lambda_{g}\\end{pmatrix}^T`, where :math:`h` denotes the equality constraints (if any) and :math:`g` denotes the inequality constraints (if any). Note that these (total) derivatives do not From 2932af680c39fcfdb2e71ccacb8e894e411fc6d4 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Mon, 24 Jul 2017 09:00:28 -0400 Subject: [PATCH 15/18] renamed beta and alpha in multi-secant to be consistent with paper --- .../algorithms/reduced_space_multi_secant.py | 10 ++++------ .../matrices/hessian/anderson_multisecant.py | 20 +++++++++---------- 2 files changed, 14 insertions(+), 16 deletions(-) diff --git a/src/kona/algorithms/reduced_space_multi_secant.py b/src/kona/algorithms/reduced_space_multi_secant.py index afefe2d..43c9e67 100644 --- a/src/kona/algorithms/reduced_space_multi_secant.py +++ b/src/kona/algorithms/reduced_space_multi_secant.py @@ -65,10 +65,8 @@ def __init__(self, primal_factory, state_factory, # set remaining options self.primal_tol_abs = get_opt(self.optns, 1e-6, 'opt_tol_abs') self.cnstr_tol_abs = get_opt(self.optns, 1e-6, 'feas_tol_abs') - self.beta = get_opt(self.optns, 1.0, 'multi_secant', 'beta') - - # TEMPORARY!!! - self.radius_max = get_opt(self.optns, 1.0, 'homotopy', 'radius_max') + self.alpha = get_opt(self.optns, 1.0, 'multi_secant', 'alpha') + self.radius_max = get_opt(self.optns, 1.0, 'multi_secant', 'radius_max') # The following data members are set by super class # self.primal_tol @@ -188,9 +186,9 @@ def solve(self): converged = True break - # get full multi-secant step and projected-gradient + # get full multi-secant step self.multisecant.build_difference_matrices(self.hess_reg) - self.multisecant.solve(R, dX, self.beta, self.precond.product) + self.multisecant.solve(R, dX, self.alpha, self.precond.product) # safe-guard against large steps info = ' ' diff --git a/src/kona/linalg/matrices/hessian/anderson_multisecant.py b/src/kona/linalg/matrices/hessian/anderson_multisecant.py index 1d88a87..cb8860e 100644 --- a/src/kona/linalg/matrices/hessian/anderson_multisecant.py +++ b/src/kona/linalg/matrices/hessian/anderson_multisecant.py @@ -75,25 +75,25 @@ def clear_history(self): self.x_hist = deque() self.r_hist = deque() - def build_difference_matrices(self, alpha=0.0): + def build_difference_matrices(self, beta=0.0): """ Constructs the difference matrices based on first-order optimality Parameters ---------- - alpha : Float + beta : Float controls the amount of Hessian regularization """ assert len(self.x_hist) == len(self.r_hist), \ 'AndersonMultiSecant() >> inconsistent list sizes!' - assert alpha >= 0.0, \ - 'AndersonMultiSecant() >> alpha must be non-negative in build_difference_matrices!' + assert beta >= 0.0, \ + 'AndersonMultiSecant() >> beta must be non-negative in build_difference_matrices!' # clear the previous x_diff and r_diff lists del self.x_diff[:], self.r_diff[:] self.work.equals_primaldual_residual(self.r_hist[0], self.x_hist[0].ineq) - if alpha > EPS: + if beta > EPS: # include Hessian regularization - self.work.primal.equals_ax_p_by(1.0, self.work.primal, alpha, self.x_hist[0].primal) + self.work.primal.equals_ax_p_by(1.0, self.work.primal, beta, self.x_hist[0].primal) for k in range(len(self.x_hist)-1): # generate new vectors for x_diff and r_diff lists dx = self._generate_vector() @@ -104,9 +104,9 @@ def build_difference_matrices(self, alpha=0.0): self.x_diff[k].equals_ax_p_by(1.0, self.x_hist[k+1], -1.0, self.x_hist[k]) self.r_diff[k].equals(self.work) self.work.equals_primaldual_residual(self.r_hist[k+1], self.x_hist[k+1].ineq) - if alpha > EPS: + if beta > EPS: # include Hessian regularization - self.work.primal.equals_ax_p_by(1.0, self.work.primal, alpha, + self.work.primal.equals_ax_p_by(1.0, self.work.primal, beta, self.x_hist[k+1].primal) self.r_diff[k].minus(self.work) self.r_diff[k].times(-1.0) @@ -138,7 +138,7 @@ def build_difference_matrices_for_homotopy(self, mu=0.0): self.r_diff[k].minus(self.work) self.r_diff[k].times(-1.0) - def solve(self, in_vec, out_vec, beta=1.0, precond=None, rel_tol=1e-15): + def solve(self, in_vec, out_vec, alpha=1.0, precond=None, rel_tol=1e-15): # store the difference matrices and rhs in numpy array format nvar = in_vec.get_num_var() dR = numpy.empty((nvar, len(self.x_diff))) @@ -162,7 +162,7 @@ def solve(self, in_vec, out_vec, beta=1.0, precond=None, rel_tol=1e-15): out_vec.equals(self.work) else: precond(self.work, out_vec) - out_vec.times(beta) + out_vec.times(alpha) sol = numpy.matmul(dX, dRinv_r) self.work.set_base_data(sol) out_vec.equals_ax_p_by(1., out_vec, -1., self.work) From 7ca799ac40dd343be255ca56b962f8ef2759b40d Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Mon, 24 Jul 2017 11:50:32 -0400 Subject: [PATCH 16/18] forgot to add synthetic data to repo for tests --- src/kona/test/synthetic_jac.dat | 10 ++++++++++ src/kona/test/synthetic_schur.dat | 10 ++++++++++ 2 files changed, 20 insertions(+) create mode 100644 src/kona/test/synthetic_jac.dat create mode 100644 src/kona/test/synthetic_schur.dat diff --git a/src/kona/test/synthetic_jac.dat b/src/kona/test/synthetic_jac.dat new file mode 100644 index 0000000..b73c0b3 --- /dev/null +++ b/src/kona/test/synthetic_jac.dat @@ -0,0 +1,10 @@ +0.040615645110386886 -0.018035062548443507 0.06670424593420433 0.008723832446152015 -0.02819360370572436 0.18772926255280378 -0.06143062176741613 0.1010770311582592 0.17617264785870174 0.00885157740710856 0.009058130571168926 0.0773294968325857 -0.033136178944166456 -0.11075142237230237 -0.008911072058798934 -0.039380262884348045 -0.044160399230537696 0.08485045215404045 0.018619585485323595 -0.05588504515549134 +-0.029006184463245016 -0.040733688102892705 0.030020751110519137 0.11698630021308515 -0.012150163414000792 0.011687132700784936 -0.07546047839992885 0.09942158986360695 0.15776193923961876 0.04319502178890986 -0.05926697822892368 0.018766069745538 -0.007261570759872578 -0.0837723778804382 -0.003862791344994239 -0.036809925947878325 0.029512260572063086 0.0312218007647646 -0.002855369369552209 -0.08575510772417029 +-0.028727107266665358 0.0629430137655778 0.09055118178804553 0.022795047107499786 -0.053848027112525884 -0.05149622266417116 -0.03912947801262293 0.0979022418013978 0.007785013759290341 -0.02586003398214672 0.018921878595235603 0.11279846810276382 -0.06832584837338347 0.027860046767303005 -0.0233588700899392 -0.038170700956331446 0.037762397683918286 0.010714347902609346 -0.01257548859955052 0.035846991619844566 +0.15858190954895449 -0.07267651833735865 0.08825508876014221 0.13002015021430527 -0.015417357822478679 0.06849523636182098 -0.07754625896787269 0.0528666253287894 0.12496655952664254 0.01925604002178876 0.030122530685916613 -0.05788294113936446 -0.06192041114461269 -0.020839697413008913 0.001753135341243263 0.010246179911848098 0.07030315570451685 0.05524738310848036 0.009475290135544446 -0.13515186213285227 +0.09568063550104519 -0.0820010273782307 -0.13236473654907954 0.10713807140096361 0.07346972333938351 0.2547406882937501 -0.08848496814454729 -0.028675488121314173 0.10096295982113138 0.12815871193505698 -0.02202229893262441 -0.19910482764783247 0.026015253966127666 -0.1921816591523116 -0.17275661814440574 0.15122323440697713 0.125575359658206 -0.07323158494383558 -0.1358131537832159 -0.13884650726601627 +-0.09479509641539914 0.08603056342109677 0.007410056767788465 -0.07035519114235869 -0.03019408175276288 -0.08626137674628144 -0.012350518455337158 -0.0476009730545892 -0.03245260095893065 -0.15128068265798233 0.02743977452179608 0.1428836057763587 -0.07872017566976491 0.04565476219212711 0.09146607967936145 -0.06678520062039717 -0.04164901658805981 0.03995394904597895 0.031697052935775266 0.004392056679850483 +-0.10790111070060804 -0.020416297655665456 -0.012938138351845214 0.04588492972601843 -0.03817998007752414 -0.006744090406322487 -0.10031624226697332 0.07749671960571432 0.043028289564393904 -0.0625128131975816 -0.09372311771059952 0.13623228326952805 -0.026778094243513694 -0.0627040486316276 -0.07932920726255133 0.0681126915101031 -0.008775070974662963 0.07943282178082772 0.05033230211871626 -0.0761127376154925 +0.00013589859315313833 0.004246629972954941 0.06818773687469637 0.045716786191511274 -0.02320482119548351 0.07473637542732657 -0.026308717529531626 0.020599624384222408 0.011586048514356836 -0.01854286007464473 -0.037025112583256595 0.10178247983966691 -0.04491359739348439 -0.060899423859043234 -0.017662741483740224 -0.027293549841328806 -0.019580640067815583 0.05860212574726484 -0.051451292520953734 -0.03377908248735524 +-0.055629091801367936 0.026410323939507947 0.02382033958662485 -0.020179267330509815 -0.03584045902081505 -0.027304646121125146 -0.030380547277023173 -0.062016978492557835 0.04392921669074589 -0.003746544912717557 -0.005585117134228435 0.056990482176091456 0.06593788257526367 0.022576033350715466 -0.002430986946676293 -0.06398251979751039 -0.0034325587534481047 0.06904848799957829 0.05261872588651662 0.10008162378048605 +0.05752757684558572 -0.1794890554222471 -0.1326445902284921 0.1963791526030564 0.025019462681374203 0.09820587238098732 -0.1710357492018347 -0.006384578570330581 0.18950154295958538 0.1535088007968005 -0.07178323017635989 -0.19337832183050174 -0.01370226430513569 -0.28087554809536946 -0.33619486922112335 0.11391364789341585 0.13821184786488902 -0.1081723497025987 -0.04130455272219124 -0.28359968708348826 diff --git a/src/kona/test/synthetic_schur.dat b/src/kona/test/synthetic_schur.dat new file mode 100644 index 0000000..81bcc1d --- /dev/null +++ b/src/kona/test/synthetic_schur.dat @@ -0,0 +1,10 @@ +30.621421089526827 -15.059149993981238 4.859071039315542 -9.594928821257 -10.287125532283007 -0.6776139778566428 -3.9225744574761516 -20.905186454857116 -8.826727341779302 8.893266875714493 +-15.059149993981238 52.94317806685921 -7.623402639356191 -10.172396647665181 12.715266075344365 8.990566411853324 -9.376391256568649 -0.19527876849702921 -5.949705672266701 -13.436158907114763 +4.859071039315541 -7.623402639356191 41.852248575297274 -3.885163138593024 11.039616970497311 0.9761850650195578 -3.8401387776511307 -19.72053935260709 2.038767420361335 0.17087463378318812 +-9.594928821257 -10.172396647665181 -3.885163138593023 24.47830855390713 4.38437181828273 1.6937816157826973 10.571964369883789 -2.138217515801466 10.405068434535501 -5.462921779384152 +-10.287125532283007 12.71526607534436 11.039616970497311 4.384371818282732 32.78154854612652 13.659862967577105 7.319751591328366 -13.770870214973954 9.070543076354886 -19.196158539757636 +-0.6776139778566446 8.990566411853324 0.9761850650195623 1.6937816157826955 13.659862967577105 36.160184371193274 -10.856016987562517 -9.58008681659966 1.8621629565448865 0.9405699790242421 +-3.9225744574761507 -9.376391256568649 -3.840138777651129 10.57196436988379 7.319751591328367 -10.856016987562514 32.18012884498713 -14.262773223119027 0.42858509499868885 -12.03696964508661 +-20.90518645485712 -0.19527876849702921 -19.72053935260709 -2.138217515801466 -13.770870214973954 -9.580086816599664 -14.262773223119027 75.41947328783678 4.1184667365477585 8.067559572926918 +-8.826727341779302 -5.949705672266704 2.0387674203613346 10.4050684345355 9.070543076354886 1.8621629565448847 0.42858509499868885 4.1184667365477585 39.07939349318166 -0.9509403284467568 +8.893266875714492 -13.43615890711476 0.1708746337831899 -5.4629217793841525 -19.196158539757636 0.940569979024243 -12.03696964508661 8.067559572926918 -0.9509403284467546 19.48411517108413 From dd17cf051f52443e8dc9f9206fb3583101241747 Mon Sep 17 00:00:00 2001 From: Jason Hicken Date: Sat, 29 Jul 2017 08:56:08 -0400 Subject: [PATCH 17/18] tried a filter --- .../algorithms/reduced_space_multi_secant.py | 41 +++++++++++++++++-- .../matrices/hessian/anderson_multisecant.py | 3 ++ src/kona/linalg/vectors/composite.py | 3 +- 3 files changed, 43 insertions(+), 4 deletions(-) diff --git a/src/kona/algorithms/reduced_space_multi_secant.py b/src/kona/algorithms/reduced_space_multi_secant.py index 43c9e67..66ae618 100644 --- a/src/kona/algorithms/reduced_space_multi_secant.py +++ b/src/kona/algorithms/reduced_space_multi_secant.py @@ -28,7 +28,7 @@ def __init__(self, primal_factory, state_factory, self.eq_factory.request_num_vectors(num_pd) if self.ineq_factory is not None: self.ineq_factory.request_num_vectors(num_pd) - self.state_factory.request_num_vectors(3) + self.state_factory.request_num_vectors(5) # iteration counter self.iter = 0 @@ -67,6 +67,7 @@ def __init__(self, primal_factory, state_factory, self.cnstr_tol_abs = get_opt(self.optns, 1e-6, 'feas_tol_abs') self.alpha = get_opt(self.optns, 1.0, 'multi_secant', 'alpha') self.radius_max = get_opt(self.optns, 1.0, 'multi_secant', 'radius_max') + self.filter = SimpleFilter() # The following data members are set by super class # self.primal_tol @@ -129,6 +130,8 @@ def solve(self): state = self.state_factory.generate() state_work = self.state_factory.generate() adjoint = self.state_factory.generate() + state_save = self.state_factory.generate() + adjoint_save = self.state_factory.generate() # evaluate the initial design, state, and adjoint before starting outer iterations X.equals_init_guess() @@ -192,11 +195,42 @@ def solve(self): # safe-guard against large steps info = ' ' - if dX.norm2 > self.radius_max: - dX.times(self.radius_max/(dX.norm2)) + if dX.primal.norm2 > self.radius_max: + dX.times(self.radius_max/(dX.primal.norm2)) info += ' length restricted' X.plus(dX) + # # evaluate at new X, construct first-order optimality conditions, and store + # state_save.equals(state) + # state.equals_primal_solution(X.primal) + # obj_val = objective_value(X.primal, state) + # # get the adjoint for dLdX + # adjoint_save.equals(adjoint) + # adjoint.equals_homotopy_adjoint(X, state, state_work) + # R.equals_KKT_conditions(X, state, adjoint) + # dLdX.equals(R) + # self.multisecant.add_to_history(X, dLdX) + # + # R.equals_primaldual_residual(dLdX, X.ineq) + # grad_norm, feas_norm = R.get_optimality_and_feasiblity() + # if self.filter.dominates(obj_val, feas_norm): + # # point is acceptable + # self.precond.linearize(X.primal, state) + # + # else: # point is not acceptable + # info += ', rejected' + # X.minus(dX) + # state.equals(state_save) + # obj_val = objective_value(X.primal, state) + # adjoint.equals(adjoint_save) + # R.equals_KKT_conditions(X, state, adjoint) + # dLdX.equals(R) + # + # # output current solution info to the user + # solver_info = current_solution(self.iter, X.primal, state, adjoint, X.get_dual()) + # if isinstance(solver_info, str): + # self.info_file.write('\n' + solver_info + '\n') + # evaluate at new X, construct first-order optimality conditions, and store state.equals_primal_solution(X.primal) obj_val = objective_value(X.primal, state) @@ -233,3 +267,4 @@ def solve(self): from kona.linalg.matrices.preconds.schur import ApproxSchur from kona.linalg.matrices.preconds import ReducedSchurPreconditioner from kona.linalg.solvers.util import EPS +from kona.algorithms.util.filter import SimpleFilter diff --git a/src/kona/linalg/matrices/hessian/anderson_multisecant.py b/src/kona/linalg/matrices/hessian/anderson_multisecant.py index cb8860e..fc17c7d 100644 --- a/src/kona/linalg/matrices/hessian/anderson_multisecant.py +++ b/src/kona/linalg/matrices/hessian/anderson_multisecant.py @@ -150,6 +150,9 @@ def solve(self, in_vec, out_vec, alpha=1.0, precond=None, rel_tol=1e-15): rhs = numpy.empty((nvar)) in_vec.get_base_data(rhs) dRinv = numpy.linalg.pinv(dR, rcond=1e-6) + if len(self.x_diff) > 0: + U, s, V = numpy.linalg.svd(dR) + print s sol = numpy.zeros_like(rhs) # sol[:] = -beta*rhs - numpy.matmul(dX - beta*dR, numpy.matmul(dRinv,rhs)) # out_vec.set_base_data(sol) diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index ed40237..e894e4b 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -459,7 +459,7 @@ def equals_primaldual_residual(self, dLdx, ineq_mult=None): \\begin{bmatrix} \\left[\\nabla_x f(x, u) - \\nabla_x h(x, u)^T \\lambda_{h} - \\nabla_x g(x, u)^T \\lambda_{g}\\right] \\\\ -h(x,u) \\\\ - -|g(x,u) - \\lambda_g| + g(x,u) + \\lambda_g + \\frac{1}{2}(|g(x,u) - \\lambda_g| - g(x,u) - \\lambda_g) \\end{bmatrix} where :math:`h(x,u)` are the equality constraints, and :math:`g(x,u)` are the @@ -491,6 +491,7 @@ def equals_primaldual_residual(self, dLdx, ineq_mult=None): if dLdx.ineq is not None: # include the inequality constraint part self.ineq.equals_mangasarian(dLdx.ineq, ineq_mult) + self.ineq.times(-1.) def equals_homotopy_residual(self, dLdx, x, init, mu=1.0): """ From c7ee37ef3cc609dd0e7c68ad4297a37d3ba3bd92 Mon Sep 17 00:00:00 2001 From: valeradonchev Date: Wed, 4 Oct 2023 13:37:19 -0400 Subject: [PATCH 18/18] Ran 2to3 on all files Used 2to3 from official Python library to change code from Python 2 to Python 3 syntax. --- src/kona/__init__.py | 10 ++-- src/kona/algorithms/__init__.py | 22 ++++----- src/kona/algorithms/composite_step_rsnk.py | 2 +- src/kona/algorithms/constrained_rsnk.py | 4 +- src/kona/algorithms/homotopy_multi_secant.py | 4 +- src/kona/algorithms/predictor_corrector.py | 2 +- .../algorithms/predictor_corrector_cnstr.py | 2 +- .../algorithms/reduced_space_multi_secant.py | 2 +- .../algorithms/reduced_space_quasi_newton.py | 2 +- src/kona/algorithms/unconstrained_rsnk.py | 2 +- src/kona/algorithms/util/__init__.py | 6 +-- src/kona/algorithms/util/filter.py | 2 +- src/kona/algorithms/util/linesearch.py | 4 +- src/kona/algorithms/util/merit.py | 2 +- src/kona/examples/__init__.py | 16 +++---- src/kona/examples/sellar.py | 2 +- src/kona/examples/spiral.py | 2 +- src/kona/linalg/__init__.py | 8 ++-- src/kona/linalg/matrices/hessian/__init__.py | 18 +++---- .../matrices/hessian/anderson_multisecant.py | 2 +- src/kona/linalg/matrices/hessian/lbfgs.py | 6 +-- src/kona/linalg/matrices/hessian/lsr1.py | 14 +++--- .../matrices/hessian/reduced_hessian.py | 4 +- .../linalg/matrices/hessian/reduced_kkt.py | 2 +- src/kona/linalg/matrices/preconds/__init__.py | 4 +- .../linalg/matrices/preconds/idf_schur.py | 2 +- .../linalg/matrices/preconds/low_rank_svd.py | 24 +++++----- src/kona/linalg/memory.py | 6 +-- src/kona/linalg/solvers/__init__.py | 4 +- src/kona/linalg/solvers/krylov/__init__.py | 12 ++--- src/kona/linalg/solvers/krylov/fgmres.py | 8 ++-- src/kona/linalg/solvers/krylov/flecs.py | 20 ++++---- src/kona/linalg/solvers/krylov/gcrot.py | 16 +++---- .../linalg/solvers/krylov/line_search_cg.py | 2 +- src/kona/linalg/solvers/krylov/stcg.py | 2 +- src/kona/linalg/solvers/util.py | 18 +++---- src/kona/linalg/vectors/common.py | 4 +- src/kona/linalg/vectors/composite.py | 28 +++++------ src/kona/optimizer.py | 2 +- src/kona/options.py | 2 +- src/kona/test/test_FLECS_solver.py | 10 ++-- src/kona/test/test_GCROT_solver.py | 2 +- src/kona/test/test_approx_schur.py | 12 ++--- src/kona/test/test_aug_lag_merit.py | 2 +- src/kona/test/test_constraint_jacobian.py | 16 +++---- src/kona/test/test_dual_vectors.py | 4 +- src/kona/test/test_equality_composite_step.py | 2 +- .../test/test_equality_predictor_corrector.py | 2 +- src/kona/test/test_krylov_util.py | 2 +- src/kona/test/test_lagrangian_hessian.py | 16 +++---- src/kona/test/test_low_rank_svd.py | 48 +++++++++---------- src/kona/test/test_primal_dual_vector.py | 2 +- src/kona/test/test_reduced_kkt_matrix.py | 20 ++++---- src/kona/test/test_verifier.py | 2 +- src/kona/user/__init__.py | 6 +-- 55 files changed, 220 insertions(+), 220 deletions(-) diff --git a/src/kona/__init__.py b/src/kona/__init__.py index 0c53ad6..3aa025d 100644 --- a/src/kona/__init__.py +++ b/src/kona/__init__.py @@ -1,5 +1,5 @@ -from optimizer import Optimizer -import examples -import algorithms -import linalg -import user +from .optimizer import Optimizer +from . import examples +from . import algorithms +from . import linalg +from . import user diff --git a/src/kona/algorithms/__init__.py b/src/kona/algorithms/__init__.py index fc1a9c6..cd8d998 100644 --- a/src/kona/algorithms/__init__.py +++ b/src/kona/algorithms/__init__.py @@ -1,11 +1,11 @@ -import base_algorithm -from reduced_space_quasi_newton import ReducedSpaceQuasiNewton -from unconstrained_rsnk import UnconstrainedRSNK -from constrained_rsnk import ConstrainedRSNK -from composite_step_rsnk import CompositeStepRSNK -from predictor_corrector import PredictorCorrector -from predictor_corrector_cnstr import PredictorCorrectorCnstr -from homotopy_multi_secant import HomotopyMultiSecant -from reduced_space_multi_secant import ReducedSpaceMultiSecant -from verifier import Verifier -import util +from . import base_algorithm +from .reduced_space_quasi_newton import ReducedSpaceQuasiNewton +from .unconstrained_rsnk import UnconstrainedRSNK +from .constrained_rsnk import ConstrainedRSNK +from .composite_step_rsnk import CompositeStepRSNK +from .predictor_corrector import PredictorCorrector +from .predictor_corrector_cnstr import PredictorCorrectorCnstr +from .homotopy_multi_secant import HomotopyMultiSecant +from .reduced_space_multi_secant import ReducedSpaceMultiSecant +from .verifier import Verifier +from . import util diff --git a/src/kona/algorithms/composite_step_rsnk.py b/src/kona/algorithms/composite_step_rsnk.py index 264c170..6325a5a 100644 --- a/src/kona/algorithms/composite_step_rsnk.py +++ b/src/kona/algorithms/composite_step_rsnk.py @@ -180,7 +180,7 @@ def solve(self): ############################### min_radius_active = False krylov_tol = 0.00095 - for i in xrange(self.max_iter): + for i in range(self.max_iter): # advance iteration counter self.iter += 1 diff --git a/src/kona/algorithms/constrained_rsnk.py b/src/kona/algorithms/constrained_rsnk.py index bfd0bed..60c4bfb 100644 --- a/src/kona/algorithms/constrained_rsnk.py +++ b/src/kona/algorithms/constrained_rsnk.py @@ -194,7 +194,7 @@ def solve(self): # BEGIN NEWTON LOOP HERE ############################### min_radius_active = False - for i in xrange(self.max_iter): + for i in range(self.max_iter): # advance iteration counter self.iter += 1 @@ -354,7 +354,7 @@ def filter_step(self, X, state, P, kkt_rhs, kkt_work, state_work, dual_work): filter_success = False min_radius_active = False max_filter_iter = 3 - for i in xrange(max_filter_iter): + for i in range(max_filter_iter): self.info_file.write('\nFilter Step : iter %i\n'%(i + 1)) X.plus(P) state_work.equals(state) # save state for reverting later diff --git a/src/kona/algorithms/homotopy_multi_secant.py b/src/kona/algorithms/homotopy_multi_secant.py index 1921a80..c7fc1fe 100644 --- a/src/kona/algorithms/homotopy_multi_secant.py +++ b/src/kona/algorithms/homotopy_multi_secant.py @@ -158,7 +158,7 @@ def solve(self): # Begin outer (homotopy) loop here converged = False - for i in xrange(self.max_iter): + for i in range(self.max_iter): self.iter += 1 self.inner = 0 @@ -222,7 +222,7 @@ def solve(self): ': inner_feas_norm0 = %e\n'%inner_feas_norm0) # inner (corrector) iterations - for k in xrange(self.inner_max_iter): + for k in range(self.inner_max_iter): self.inner += 1 self.info_file.write( 'Corr. Iter. %i'%self.inner + diff --git a/src/kona/algorithms/predictor_corrector.py b/src/kona/algorithms/predictor_corrector.py index 23c5e5f..93be034 100644 --- a/src/kona/algorithms/predictor_corrector.py +++ b/src/kona/algorithms/predictor_corrector.py @@ -231,7 +231,7 @@ def solve(self): max_newton = self.inner_maxiter inner_iters = 0 dx_newt.equals(0.0) - for i in xrange(max_newton): + for i in range(max_newton): self.info_file.write('\n') self.info_file.write(' Inner Newton iteration %i\n'%(i+1)) diff --git a/src/kona/algorithms/predictor_corrector_cnstr.py b/src/kona/algorithms/predictor_corrector_cnstr.py index 1a047c4..06b52b7 100644 --- a/src/kona/algorithms/predictor_corrector_cnstr.py +++ b/src/kona/algorithms/predictor_corrector_cnstr.py @@ -378,7 +378,7 @@ def solve(self): max_newton = self.inner_maxiter inner_iters = 0 dx_newt.equals(0.0) - for i in xrange(max_newton): + for i in range(max_newton): self.info_file.write('\n') self.info_file.write(' Inner Newton iteration %i\n'%(i+1)) diff --git a/src/kona/algorithms/reduced_space_multi_secant.py b/src/kona/algorithms/reduced_space_multi_secant.py index 66ae618..f47bd19 100644 --- a/src/kona/algorithms/reduced_space_multi_secant.py +++ b/src/kona/algorithms/reduced_space_multi_secant.py @@ -165,7 +165,7 @@ def solve(self): # Begin iterations converged = False info = ' ' - for i in xrange(self.max_iter): + for i in range(self.max_iter): self.iter += 1 # evaluate optimality and feasibility, and output as necessary diff --git a/src/kona/algorithms/reduced_space_quasi_newton.py b/src/kona/algorithms/reduced_space_quasi_newton.py index 8a879d2..c92a385 100644 --- a/src/kona/algorithms/reduced_space_quasi_newton.py +++ b/src/kona/algorithms/reduced_space_quasi_newton.py @@ -92,7 +92,7 @@ def solve(self): self.iter = 0 converged = False self._write_header() - for i in xrange(self.max_iter): + for i in range(self.max_iter): info.write('========== Outer Iteration %i ==========\n'%(i+1)) if not state.equals_primal_solution(x): info.write('WARNING: Nonlinear solution failed to converge!\n') diff --git a/src/kona/algorithms/unconstrained_rsnk.py b/src/kona/algorithms/unconstrained_rsnk.py index 2d34530..fbeb88a 100644 --- a/src/kona/algorithms/unconstrained_rsnk.py +++ b/src/kona/algorithms/unconstrained_rsnk.py @@ -162,7 +162,7 @@ def solve(self): self._write_header(obj_scale) converged = False grad_tol = self.primal_tol*grad_norm0 - for i in xrange(self.max_iter): + for i in range(self.max_iter): self.info_file.write( '==================================================\n') diff --git a/src/kona/algorithms/util/__init__.py b/src/kona/algorithms/util/__init__.py index ae5b50c..d2f4c88 100644 --- a/src/kona/algorithms/util/__init__.py +++ b/src/kona/algorithms/util/__init__.py @@ -1,3 +1,3 @@ -import merit -import linesearch -import filter +from . import merit +from . import linesearch +from . import filter diff --git a/src/kona/algorithms/util/filter.py b/src/kona/algorithms/util/filter.py index b0fb12f..51d14ef 100644 --- a/src/kona/algorithms/util/filter.py +++ b/src/kona/algorithms/util/filter.py @@ -34,7 +34,7 @@ def dominates(self, obj, cnstr_norm): # check if new point dominates any filter point dominated_by = [] - for i in xrange(len(self.points)): + for i in range(len(self.points)): if new_point >= self.points[i]: # new point is not acceptable return False diff --git a/src/kona/algorithms/util/linesearch.py b/src/kona/algorithms/util/linesearch.py index 399248e..05b3f2a 100644 --- a/src/kona/algorithms/util/linesearch.py +++ b/src/kona/algorithms/util/linesearch.py @@ -228,7 +228,7 @@ def _zoom(self, alpha_low, alpha_hi, phi_low, phi_hi, dphi_low, dphi_hi): self.out_file.write('\n') # START OF BIG FOR-LOOP - for i in xrange(self.max_iter): + for i in range(self.max_iter): self.out_file.write(' Strong-Wolfe Zoom : iter %i\n'%(i+1)) @@ -305,7 +305,7 @@ def find_step_length(self, merit): self.out_file.write('\n') # START OF BIG FOR-LOOP - for i in xrange(self.max_iter): + for i in range(self.max_iter): self.out_file.write(' Strong-Wolfe Line Search : iter %i\n'%(i+1)) self.out_file.write(' ----------------------------------\n') diff --git a/src/kona/algorithms/util/merit.py b/src/kona/algorithms/util/merit.py index 7627b29..5b2a6e7 100644 --- a/src/kona/algorithms/util/merit.py +++ b/src/kona/algorithms/util/merit.py @@ -409,7 +409,7 @@ def eval_func(self, alpha): else: self.mult_eq.equals(self.mult_eq_start) self.func_val += self.mult_eq.inner(self.cnstr_eq) - print self.func_val + print(self.func_val) if self.cnstr_ineq is not None: if not self.freeze_mults: self.mult_ineq.equals_ax_p_by( diff --git a/src/kona/examples/__init__.py b/src/kona/examples/__init__.py index a133645..503e6ce 100644 --- a/src/kona/examples/__init__.py +++ b/src/kona/examples/__init__.py @@ -1,8 +1,8 @@ -from rosenbrock import Rosenbrock -from sphere_constrained import SphereConstrained -from exponential_constrained import ExponentialConstrained -from simple_2by2 import Simple2x2 -from constrained_2by2 import Constrained2x2 -from spiral import Spiral -from sellar import Sellar -from simple_mdo import SimpleMDF, SimpleIDF +from .rosenbrock import Rosenbrock +from .sphere_constrained import SphereConstrained +from .exponential_constrained import ExponentialConstrained +from .simple_2by2 import Simple2x2 +from .constrained_2by2 import Constrained2x2 +from .spiral import Spiral +from .sellar import Sellar +from .simple_mdo import SimpleMDF, SimpleIDF diff --git a/src/kona/examples/sellar.py b/src/kona/examples/sellar.py index 1cc58a7..f44b556 100644 --- a/src/kona/examples/sellar.py +++ b/src/kona/examples/sellar.py @@ -175,7 +175,7 @@ def solve_nonlinear(self, at_design, result): u += du u[0] = max(0.01, u[0]) - print 'Sellar nonlinear problem failed to converge' + print('Sellar nonlinear problem failed to converge') result.data[0] = u[0] result.data[1] = u[1] return -1 diff --git a/src/kona/examples/spiral.py b/src/kona/examples/spiral.py index ff52e0f..3ac937a 100644 --- a/src/kona/examples/spiral.py +++ b/src/kona/examples/spiral.py @@ -114,7 +114,7 @@ def solve_nonlinear(self, at_design, result): max_iter = 100 converged = False # start Newton loop - for iters in xrange(max_iter): + for iters in range(max_iter): # check convergence with new residual if np.linalg.norm(self.PDE.R) <= rel_tol*norm0: converged = True diff --git a/src/kona/linalg/__init__.py b/src/kona/linalg/__init__.py index d95fe3d..128190c 100644 --- a/src/kona/linalg/__init__.py +++ b/src/kona/linalg/__init__.py @@ -1,4 +1,4 @@ -import vectors -import matrices -import solvers -import common +from . import vectors +from . import matrices +from . import solvers +from . import common diff --git a/src/kona/linalg/matrices/hessian/__init__.py b/src/kona/linalg/matrices/hessian/__init__.py index 17c6898..78793ac 100644 --- a/src/kona/linalg/matrices/hessian/__init__.py +++ b/src/kona/linalg/matrices/hessian/__init__.py @@ -1,12 +1,12 @@ -import basic +from . import basic -from lbfgs import LimitedMemoryBFGS -from lsr1 import LimitedMemorySR1 -from anderson_multisecant import AndersonMultiSecant +from .lbfgs import LimitedMemoryBFGS +from .lsr1 import LimitedMemorySR1 +from .anderson_multisecant import AndersonMultiSecant -from reduced_hessian import ReducedHessian -from reduced_kkt import ReducedKKTMatrix +from .reduced_hessian import ReducedHessian +from .reduced_kkt import ReducedKKTMatrix -from constraint_jacobian import TotalConstraintJacobian -from augmented_kkt_matrix import AugmentedKKTMatrix -from lagrangian_hessian import LagrangianHessian +from .constraint_jacobian import TotalConstraintJacobian +from .augmented_kkt_matrix import AugmentedKKTMatrix +from .lagrangian_hessian import LagrangianHessian diff --git a/src/kona/linalg/matrices/hessian/anderson_multisecant.py b/src/kona/linalg/matrices/hessian/anderson_multisecant.py index fc17c7d..12d6dc1 100644 --- a/src/kona/linalg/matrices/hessian/anderson_multisecant.py +++ b/src/kona/linalg/matrices/hessian/anderson_multisecant.py @@ -152,7 +152,7 @@ def solve(self, in_vec, out_vec, alpha=1.0, precond=None, rel_tol=1e-15): dRinv = numpy.linalg.pinv(dR, rcond=1e-6) if len(self.x_diff) > 0: U, s, V = numpy.linalg.svd(dR) - print s + print(s) sol = numpy.zeros_like(rhs) # sol[:] = -beta*rhs - numpy.matmul(dX - beta*dR, numpy.matmul(dRinv,rhs)) # out_vec.set_base_data(sol) diff --git a/src/kona/linalg/matrices/hessian/lbfgs.py b/src/kona/linalg/matrices/hessian/lbfgs.py index 3ade7b7..232adee 100644 --- a/src/kona/linalg/matrices/hessian/lbfgs.py +++ b/src/kona/linalg/matrices/hessian/lbfgs.py @@ -68,12 +68,12 @@ def solve(self, u_vec, v_vec, rel_tol=1e-15): rho = numpy.zeros(num_stored) alpha = numpy.zeros(num_stored) - for k in xrange(num_stored): + for k in range(num_stored): rho[k] = 1.0 / (s_dot_s_list[k] * lambda0 + s_dot_y_list[k]) v_vec.equals(u_vec) - for k in xrange(num_stored-1, -1, -1): + for k in range(num_stored-1, -1, -1): alpha[k] = rho[k] * s_list[k].inner(v_vec) if lambda0 > 0.0: v_vec.equals_ax_p_by(1.0, v_vec, @@ -91,7 +91,7 @@ def solve(self, u_vec, v_vec, rel_tol=1e-15): else: v_vec.divide_by(self.norm_init) - for k in xrange(num_stored): + for k in range(num_stored): beta = rho[k] * y_list[k].inner(v_vec) if lambda0 > 0.0: beta += rho[k] * lambda0 * s_list[k].inner(v_vec) diff --git a/src/kona/linalg/matrices/hessian/lsr1.py b/src/kona/linalg/matrices/hessian/lsr1.py index 7f3ce83..8bab84d 100644 --- a/src/kona/linalg/matrices/hessian/lsr1.py +++ b/src/kona/linalg/matrices/hessian/lsr1.py @@ -100,18 +100,18 @@ def product(self, u_vec, v_vec): num_stored = len(s_list) Bs = [] - for k in xrange(num_stored): + for k in range(num_stored): Bs.append(self.vec_fac.generate()) Bs[k].equals(s_list[k]) v_vec.equals(u_vec) - for i in xrange(num_stored): + for i in range(num_stored): denom = 1.0 / (y_list[i].inner(s_list[i]) - Bs[i].inner(s_list[i])) fac = (y_list[i].inner(u_vec) - Bs[i].inner(u_vec))*denom v_vec.equals_ax_p_by(1.0, v_vec, fac, y_list[i]) v_vec.equals_ax_p_by(1.0, v_vec, -fac, Bs[i]) - for j in xrange(i+1, num_stored): + for j in range(i+1, num_stored): fac = (y_list[i].inner(s_list[j]) - Bs[i].inner(s_list[j]))*denom Bs[j].equals_ax_p_by(1.0, Bs[j], fac, y_list[i]) @@ -133,7 +133,7 @@ def solve(self, u_vec, v_vec, rel_tol=1e-15): alpha = numpy.zeros(num_stored) z_list = [] - for k in xrange(num_stored): + for k in range(num_stored): z_list.append(self.vec_fac.generate()) z_list[k].equals_ax_p_by(1.0 - lambda0, s_list[k], (lambda0 - 1.0)/norm_init, @@ -141,8 +141,8 @@ def solve(self, u_vec, v_vec, rel_tol=1e-15): alpha = self._check_threshold(z_list, 0, alpha) - for i in xrange(1, num_stored): - for j in xrange(1, num_stored): + for i in range(1, num_stored): + for j in range(1, num_stored): prod = (1.0 - lambda0) * z_list[i-1].inner(y_list[j]) + \ lambda0 * norm_init * z_list[i-1].inner(s_list[j]) z_list[j].equals_ax_p_by(1.0, z_list[j], @@ -151,6 +151,6 @@ def solve(self, u_vec, v_vec, rel_tol=1e-15): alpha = self._check_threshold(z_list, i, alpha) v_vec.equals(u_vec) - for k in xrange(num_stored): + for k in range(num_stored): v_vec.equals_ax_p_by(1.0, v_vec, alpha[k] * z_list[k].inner(u_vec), z_list[k]) diff --git a/src/kona/linalg/matrices/hessian/reduced_hessian.py b/src/kona/linalg/matrices/hessian/reduced_hessian.py index 3fb032e..4982d62 100644 --- a/src/kona/linalg/matrices/hessian/reduced_hessian.py +++ b/src/kona/linalg/matrices/hessian/reduced_hessian.py @@ -94,14 +94,14 @@ def linearize(self, at_design, at_state, at_adjoint, scale=1.0): self.w_adj = self.state_factory.generate() self.lambda_adj = self.state_factory.generate() self.state_work = [] - for i in xrange(4): + for i in range(4): self.state_work.append(self.state_factory.generate()) # generate primal vectors self.pert_design = self.primal_factory.generate() self.reduced_grad = self.primal_factory.generate() self.primal_work = [] - for i in xrange(2): + for i in range(2): self.primal_work.append(self.primal_factory.generate()) self._allocated = True diff --git a/src/kona/linalg/matrices/hessian/reduced_kkt.py b/src/kona/linalg/matrices/hessian/reduced_kkt.py index af2222d..8edeed1 100644 --- a/src/kona/linalg/matrices/hessian/reduced_kkt.py +++ b/src/kona/linalg/matrices/hessian/reduced_kkt.py @@ -126,7 +126,7 @@ def linearize(self, at_kkt, at_state, at_adjoint, self.w_adj = self.state_factory.generate() self.lambda_adj = self.state_factory.generate() self.state_work = [] - for i in xrange(3): + for i in range(3): self.state_work.append(self.state_factory.generate()) # generate primal vectors diff --git a/src/kona/linalg/matrices/preconds/__init__.py b/src/kona/linalg/matrices/preconds/__init__.py index 5b1bbd1..4065357 100644 --- a/src/kona/linalg/matrices/preconds/__init__.py +++ b/src/kona/linalg/matrices/preconds/__init__.py @@ -1,2 +1,2 @@ -from idf_schur import ReducedSchurPreconditioner -from low_rank_svd import LowRankSVD +from .idf_schur import ReducedSchurPreconditioner +from .low_rank_svd import LowRankSVD diff --git a/src/kona/linalg/matrices/preconds/idf_schur.py b/src/kona/linalg/matrices/preconds/idf_schur.py index 19e96f1..90ae5e1 100644 --- a/src/kona/linalg/matrices/preconds/idf_schur.py +++ b/src/kona/linalg/matrices/preconds/idf_schur.py @@ -88,7 +88,7 @@ def linearize(self, at_primal, at_state, scale=1.0): # design vectors self.design_prod = self.primal_factory.generate() self.design_work = [] - for i in xrange(2): + for i in range(2): self.design_work.append(self.primal_factory.generate()) # dual vectors self.dual_prod = None diff --git a/src/kona/linalg/matrices/preconds/low_rank_svd.py b/src/kona/linalg/matrices/preconds/low_rank_svd.py index accbe99..8e052da 100644 --- a/src/kona/linalg/matrices/preconds/low_rank_svd.py +++ b/src/kona/linalg/matrices/preconds/low_rank_svd.py @@ -69,7 +69,7 @@ def linearize(self): self.V = [] self.P = [] self.U = [] - for i in xrange(self.subspace_size): + for i in range(self.subspace_size): self.Q.append(self.fwd_factory.generate()) self.V.append(self.fwd_factory.generate()) self.P.append(self.rev_factory.generate()) @@ -97,47 +97,47 @@ def linearize(self): self.S = np.diag(s_tmp) # calculate V = Q*v_tmp - for j in xrange(len(self.V)): + for j in range(len(self.V)): self.V[j].equals(0.0) - for i in xrange(len(v_tmp[:, j])): + for i in range(len(v_tmp[:, j])): self.q_work.equals(self.Q[i]) self.q_work.times(v_tmp[i, j]) self.V[j].plus(self.q_work) # calculate U = P*u_tmp - for j in xrange(len(self.U)): + for j in range(len(self.U)): self.U[j].equals(0.0) - for i in xrange(len(u_tmp[:, j])): + for i in range(len(u_tmp[:, j])): self.p_work.equals(self.P[i]) self.p_work.times(u_tmp[i, j]) self.U[j].plus(self.p_work) def approx_fwd_prod(self, in_vec, out_vec): VT_in = np.zeros(len(self.V)) - for i in xrange(len(self.V)): + for i in range(len(self.V)): VT_in[i] = self.V[i].inner(in_vec) SVT_in = np.dot(self.S, VT_in) out_vec.equals(0.0) - for i in xrange(len(self.U)): + for i in range(len(self.U)): out_vec.equals_ax_p_by(1., out_vec, SVT_in[i], self.U[i]) def approx_rev_prod(self, in_vec, out_vec): UT_vec = np.zeros(len(self.U)) - for i in xrange(len(self.U)): + for i in range(len(self.U)): UT_vec[i] = self.U[i].inner(in_vec) SUT_vec = np.dot(self.S, UT_vec) out_vec.equals(0.0) - for i in xrange(len(self.V)): + for i in range(len(self.V)): out_vec.equals_ax_p_by(1., out_vec, SUT_vec[i], self.V[i]) def inv_schur_prod(self, in_vec, out_vec): UT_vec = np.zeros(len(self.U)-1) - for i in xrange(len(self.U)-1): + for i in range(len(self.U)-1): UT_vec[i] = self.U[i].inner(in_vec) S2invUT_vec = np.zeros(len(self.U)-1) - for i in xrange(len(self.U)-1): + for i in range(len(self.U)-1): S2invUT_vec[i] = UT_vec[i]*(1./self.S[i,i]**2 - 1./self.S[-1,-1]**2) out_vec.equals(in_vec) out_vec.times(1./self.S[-1,-1]**2) - for i in xrange(len(self.U)-1): + for i in range(len(self.U)-1): out_vec.equals_ax_p_by(1.0, out_vec, S2invUT_vec[i], self.U[i]) diff --git a/src/kona/linalg/memory.py b/src/kona/linalg/memory.py index 6f05729..65ca029 100644 --- a/src/kona/linalg/memory.py +++ b/src/kona/linalg/memory.py @@ -24,7 +24,7 @@ class VectorFactory(object): def __init__(self, memory, vec_type=None): self.num_vecs = 0 self._memory = memory - if vec_type not in self._memory.vector_stack.keys(): + if vec_type not in list(self._memory.vector_stack.keys()): raise TypeError('VectorFactory() >> Unknown vector type!') else: self._vec_type = vec_type @@ -70,7 +70,7 @@ def __init__(self, filename, rank): # only produce a file handle for the root processor if rank == 0: if isinstance(filename, str): - self.file = open(filename, 'w', 0) + self.file = open(filename, 'w') else: self.file = filename else: @@ -181,7 +181,7 @@ def pop_vector(self, vec_type): BaseVector User-defined vector data structure. """ - if vec_type not in self.vector_stack.keys(): + if vec_type not in list(self.vector_stack.keys()): raise TypeError('KonaMemory.pop_vector() >> ' + 'Unknown vector type!') else: diff --git a/src/kona/linalg/solvers/__init__.py b/src/kona/linalg/solvers/__init__.py index b8343ce..4bab5bf 100644 --- a/src/kona/linalg/solvers/__init__.py +++ b/src/kona/linalg/solvers/__init__.py @@ -1,2 +1,2 @@ -import krylov -import util +from . import krylov +from . import util diff --git a/src/kona/linalg/solvers/krylov/__init__.py b/src/kona/linalg/solvers/krylov/__init__.py index 308b555..99392e4 100644 --- a/src/kona/linalg/solvers/krylov/__init__.py +++ b/src/kona/linalg/solvers/krylov/__init__.py @@ -1,6 +1,6 @@ -import basic -from stcg import STCG -from flecs import FLECS -from fgmres import FGMRES -from gcrot import GCROT -from line_search_cg import LineSearchCG +from . import basic +from .stcg import STCG +from .flecs import FLECS +from .fgmres import FGMRES +from .gcrot import GCROT +from .line_search_cg import LineSearchCG diff --git a/src/kona/linalg/solvers/krylov/fgmres.py b/src/kona/linalg/solvers/krylov/fgmres.py index d067b34..0bddeb8 100644 --- a/src/kona/linalg/solvers/krylov/fgmres.py +++ b/src/kona/linalg/solvers/krylov/fgmres.py @@ -93,7 +93,7 @@ def solve(self, mat_vec, b, x, precond): ################ lin_depend = False - for i in xrange(self.max_iter): + for i in range(self.max_iter): # check convergence and linear dependence if lin_depend and (beta > self.rel_tol*norm0): @@ -123,7 +123,7 @@ def solve(self, mat_vec, b, x, precond): # apply old Givens rotations to new column of the Hessenberg matrix # then generate new Givens rotation matrix and apply it to the last # two elements of H[i, :] and g - for k in xrange(i): + for k in range(i): H[k, i], H[k+1, i] = apply_givens( sn[k], cn[k], H[k, i], H[k+1, i]) @@ -135,7 +135,7 @@ def solve(self, mat_vec, b, x, precond): if self.check_LSgrad and iters > 1: # check the gradient of the least-squares problem y[:i] = numpy.zeros(i) - for k in xrange(i-1, -1, -1): + for k in range(i-1, -1, -1): y[k], y[k+1] = apply_givens(-sn[k], cn[k], y[k], y[k+1]) rLS = numpy.dot(H[:i+1,:i+1], y[:i+1]) if numpy.sqrt(rLS.dot(rLS)) < 1000*EPS: @@ -152,7 +152,7 @@ def solve(self, mat_vec, b, x, precond): # solve the least squares system y[:i] = solve_tri(H[:i, :i], g[:i], lower=False) - for k in xrange(i): + for k in range(i): x.equals_ax_p_by(1.0, x, y[k], Z[k]) if self.check_res: diff --git a/src/kona/linalg/solvers/krylov/flecs.py b/src/kona/linalg/solvers/krylov/flecs.py index 325fbc8..c17054e 100644 --- a/src/kona/linalg/solvers/krylov/flecs.py +++ b/src/kona/linalg/solvers/krylov/flecs.py @@ -187,13 +187,13 @@ def apply_correction(self, cnstr, step): VtVH = VtV_dual_r.dot(H_r) A = ZtZ_prim_r + VtVH + VtVH.T rhs = numpy.zeros(self.iters) - for k in xrange(self.iters): + for k in range(self.iters): rhs[k] = self.V[k].dual.inner(cnstr) self.y = numpy.linalg.solve(A, rhs) # construct the design update # leave the dual solution untouched - for k in xrange(self.iters): + for k in range(self.iters): step.primal.equals_ax_p_by( 1.0, step.primal, self.y[k], self.Z[k].primal) @@ -241,7 +241,7 @@ def solve_subspace_problems(self): # compute the RHS for the augmented Lagrangian problem rhs_aug = numpy.zeros(self.iters) - for k in xrange(self.iters): + for k in range(self.iters): rhs_aug[k] = -self.g[0]*(self.VtZ_prim[0, k] + self.mu*VtVH[0, k]) radius_aug = self.radius @@ -253,12 +253,12 @@ def solve_subspace_problems(self): L = numpy.linalg.cholesky(ZtZ_prim_r) rhs_aug = solve_tri(L, rhs_tmp, lower=True) - for j in xrange(self.iters): + for j in range(self.iters): rhs_tmp[:] = Hess_aug[:,j] vec_tmp = solve_tri(L, rhs_tmp, lower=True) Hess_aug[:, j] = vec_tmp[:] - for j in xrange(self.iters): + for j in range(self.iters): rhs_tmp[:] = Hess_aug[j,:] vec_tmp = solve_tri(L, rhs_tmp, lower=True) Hess_aug[j,:] = vec_tmp[:] @@ -361,7 +361,7 @@ def solve(self, mat_vec, b, x, precond): self.lin_depend = False self.neg_curv = False self.trust_active = False - for i in xrange(self.max_iter): + for i in range(self.max_iter): # advance iteration counter self.iters += 1 @@ -386,7 +386,7 @@ def solve(self, mat_vec, b, x, precond): self.lin_depend = True # compute new row and column of the VtZ matrix - for k in xrange(i+1): + for k in range(i+1): self.VtZ_prim[k, i] = self.V[k].primal.inner(self.Z[i].primal) self.VtZ_prim[i+1, k] = self.V[i+1].primal.inner( self.Z[k].primal) @@ -440,7 +440,7 @@ def solve(self, mat_vec, b, x, precond): # compute solution: augmented-Lagrangian for primal, FGMRES for dual x.equals(0.0) - for k in xrange(self.iters): + for k in range(self.iters): x.primal.equals_ax_p_by( 1.0, x.primal, self.y_aug[k], self.Z[k].primal) x.dual.equals_ax_p_by( @@ -454,7 +454,7 @@ def solve(self, mat_vec, b, x, precond): if self.check_res: # calculate true residual for the solution self.V[0].equals(0.0) - for k in xrange(self.iters): + for k in range(self.iters): self.V[0].equals_ax_p_by( 1.0, self.V[0], self.y_mult[k], self.Z[k]) mat_vec(self.V[0], res) @@ -526,7 +526,7 @@ def re_solve(self, b, x): # always use composite-step approach in re-solve x.equals(0.0) - for k in xrange(self.iters): + for k in range(self.iters): x.primal.equals_ax_p_by( 1.0, x.primal, self.y_aug[k], self.Z[k].primal) x.dual.equals_ax_p_by( diff --git a/src/kona/linalg/solvers/krylov/gcrot.py b/src/kona/linalg/solvers/krylov/gcrot.py index edac27b..8f9ae94 100644 --- a/src/kona/linalg/solvers/krylov/gcrot.py +++ b/src/kona/linalg/solvers/krylov/gcrot.py @@ -127,7 +127,7 @@ def solve(self, mat_vec, b, x, precond): norm0 = res.norm2 # find initial guess from recycled subspace - for k in xrange(self.num_stored): + for k in range(self.num_stored): alpha = res.inner(self.C[k]) x.equals_ax_p_by(1.0, x, alpha, self.U[k]) res.equals_ax_p_by(1.0, res, -alpha, self.C[k]) @@ -146,7 +146,7 @@ def solve(self, mat_vec, b, x, precond): ########################## W = [] Z = [] - for j in xrange(self.max_outer): + for j in range(self.max_outer): # begin nested FGMRES(fgmres_iter) fgmres_iter = self.max_recycle - self.num_stored + self.max_iter @@ -168,7 +168,7 @@ def solve(self, mat_vec, b, x, precond): g[0] = beta inner_iters = 0 lin_depend = False - for i in xrange(fgmres_iter): + for i in range(fgmres_iter): # check convergence and linear dependence if lin_depend and (beta > self.rel_tol*norm0): @@ -208,7 +208,7 @@ def solve(self, mat_vec, b, x, precond): # apply old Givens rotations to new column of the Hessenberg # matrix then generate new Givens rotation matrix and apply it # to the last two elements of H[i, :] and g - for k in xrange(i): + for k in range(i): H[k, i], H[k+1, i] = apply_givens( sn[k], cn[k], H[k, i], H[k+1, i]) @@ -228,20 +228,20 @@ def solve(self, mat_vec, b, x, precond): # first, solve to get y = R^{-1} g y[:i] = solve_tri(H[:i, :i], g[:i], lower=False) U_new.equals(0.0) - for k in xrange(i): + for k in range(i): U_new.equals_ax_p_by(1.0, U_new, y[k], Z[k]) # update U_new -= U * B - for k in xrange(self.num_stored): + for k in range(self.num_stored): tmp = numpy.dot(B[k, :i], y[:i]) U_new.equals_ax_p_by(1.0, U_new, -tmp, self.U[k]) # finished with g, so undo rotations to find C_new y[:i] = g[:i] y[i] = 0.0 - for k in xrange(i-1, -1, -1): + for k in range(i-1, -1, -1): y[k], y[k+1] = apply_givens(-sn[k], cn[k], y[k], y[k+1]) C_new.equals(0.0) - for k in xrange(i+1): + for k in range(i+1): C_new.equals_ax_p_by(1.0, C_new, y[k], W[k]) # normalize and scale new vectors and update solution and res diff --git a/src/kona/linalg/solvers/krylov/line_search_cg.py b/src/kona/linalg/solvers/krylov/line_search_cg.py index 4e85310..0fd6c7e 100644 --- a/src/kona/linalg/solvers/krylov/line_search_cg.py +++ b/src/kona/linalg/solvers/krylov/line_search_cg.py @@ -39,7 +39,7 @@ def solve(self, mat_vec, neg_grad, p, precond=None): # START OF BIG FOR LOOP ####################### - for i in xrange(self.max_iter): + for i in range(self.max_iter): mat_vec(d, Bd) curv = d.inner(Bd) diff --git a/src/kona/linalg/solvers/krylov/stcg.py b/src/kona/linalg/solvers/krylov/stcg.py index 1455500..5d040a9 100644 --- a/src/kona/linalg/solvers/krylov/stcg.py +++ b/src/kona/linalg/solvers/krylov/stcg.py @@ -73,7 +73,7 @@ def solve(self, mat_vec, b, x, precond): # START OF BIG FOR LOOP ####################### - for i in xrange(self.max_iter): + for i in range(self.max_iter): # to be included from c++ # if (ptin.get("dynamic",false)) diff --git a/src/kona/linalg/solvers/util.py b/src/kona/linalg/solvers/util.py index 3e63306..c1e25d4 100644 --- a/src/kona/linalg/solvers/util.py +++ b/src/kona/linalg/solvers/util.py @@ -182,7 +182,7 @@ def lanczos_tridiag(mat_vec, Q, Q_init=False): Q[0].equals(1.0) Q[0].divide_by(Q[0].norm2) - for i in xrange(subspace_size-1): + for i in range(subspace_size-1): # perform the matrix vector product mat_vec(Q[i], Q[i+1]) # orthogonalize the new vector @@ -248,7 +248,7 @@ def lanczos_bidiag(fwd_mat_vec, Q, q_work, Q[0].equals(1.0) Q[0].divide_by(Q[0].norm2) - for j in xrange(subspace_size): + for j in range(subspace_size): fwd_mat_vec(Q[j], P[j]) if j > 0: p_work.equals(P[j-1]) @@ -333,8 +333,8 @@ def solve_trust_reduced(H, g, radius): if (fnc < 0.0): # i.e. norm_2(y) < raidus # compute predicted decrease in objective and return pred = 0.0 - for i in xrange(n): - for j in xrange(n): + for i in range(n): + for j in range(n): pred -= 0.5*y[i]*H[i,j]*y[j] pred -= g[i]*y[i] return y, lam, pred @@ -346,7 +346,7 @@ def solve_trust_reduced(H, g, radius): dlam = 0.1*max(-eigmin, EPS) lam_h = max(-eigmin, 0.0) + dlam y, fnc_h, dfnc = secular_function(H, g, lam_h, radius) - for k in xrange(max_brk): + for k in range(max_brk): if fnc_h > 0.0: break dlam *= 0.1 @@ -357,7 +357,7 @@ def solve_trust_reduced(H, g, radius): dlam = np.sqrt(EPS) lam_l = max(-eigmin, 0.0) + dlam y, fnc_l, dfnc = secular_function(H, g, lam_l, radius) - for k in xrange(max_brk): + for k in range(max_brk): if fnc_l < 0.0: break dlam *= 100.0 @@ -376,7 +376,7 @@ def solve_trust_reduced(H, g, radius): y, fnc, dfnc = secular_function(H, g, lam, radius) res0 = abs(fnc) - for l in xrange(max_Newt): + for l in range(max_Newt): # check if y lies on the trust region; if so, exit loop if (abs(fnc) < tol*res0) or (abs(dlam) < lam_tol): break @@ -500,7 +500,7 @@ def mod_GS_normalize(i, Hsbg, w): return # begin main Gram-Schmidt loop - for k in xrange(i+1): + for k in range(i+1): prod = w[i+1].inner(w[k]) Hsbg[k, i] = prod w[i+1].equals_ax_p_by(1.0, w[i+1], -prod, w[k]) @@ -552,7 +552,7 @@ def mod_gram_schmidt(i, B, C, w, normalize=False): raise ValueError('mod_gram_schmidt failed : w = NaN') # begin main Gram-Schmidt loop - for k in xrange(len(C)): + for k in range(len(C)): prod = w.inner(C[k]) B[k, i] = prod w.equals_ax_p_by(1.0, w, -prod, C[k]) diff --git a/src/kona/linalg/vectors/common.py b/src/kona/linalg/vectors/common.py index dc704f6..4d1b737 100644 --- a/src/kona/linalg/vectors/common.py +++ b/src/kona/linalg/vectors/common.py @@ -276,11 +276,11 @@ def enforce_bounds(self): Element-wise enforcement of design bounds. """ if self.lb is not None: - for i in xrange(len(self.base.data)): + for i in range(len(self.base.data)): if self.base.data[i] < self.lb: self.base.data[i] = self.lb if self.ub is not None: - for i in xrange(len(self.base.data)): + for i in range(len(self.base.data)): if self.base.data[i] > self.ub: self.base.data[i] = self.ub diff --git a/src/kona/linalg/vectors/composite.py b/src/kona/linalg/vectors/composite.py index e894e4b..da51fab 100644 --- a/src/kona/linalg/vectors/composite.py +++ b/src/kona/linalg/vectors/composite.py @@ -109,7 +109,7 @@ def _check_type(self, vec): raise TypeError('CompositeVector() >> ' + 'Wrong vector type. Must be %s' % type(self)) else: - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): try: self._vectors[i]._check_type(vec._vectors[i]) except TypeError: @@ -131,11 +131,11 @@ def equals(self, rhs): """ if isinstance(rhs, (float, int, np.float64, np.int64, np.float32, np.int32)): - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): self._vectors[i].equals(rhs) else: self._check_type(rhs) - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): self._vectors[i].equals(rhs._vectors[i]) def plus(self, vector): @@ -150,7 +150,7 @@ def plus(self, vector): Vector to be added. """ self._check_type(vector) - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): self._vectors[i].plus(vector._vectors[i]) def minus(self, vector): @@ -165,7 +165,7 @@ def minus(self, vector): Vector to be subtracted. """ self._check_type(vector) - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): self._vectors[i].minus(vector._vectors[i]) def times(self, factor): @@ -181,11 +181,11 @@ def times(self, factor): """ if isinstance(factor, (float, int, np.float64, np.int64, np.float32, np.int32)): - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): self._vectors[i].times(factor) else: self._check_type(factor) - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): self._vectors[i].times(factor._vectors[i]) def divide_by(self, value): @@ -201,7 +201,7 @@ def divide_by(self, value): """ if isinstance(value, (float, int, np.float64, np.int64, np.float32, np.int32)): - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): self._vectors[i].divide_by(value) else: raise TypeError( @@ -221,7 +221,7 @@ def equals_ax_p_by(self, a, x, b, y): """ self._check_type(x) self._check_type(y) - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): self._vectors[i].equals_ax_p_by(a, x._vectors[i], b, y._vectors[i]) def inner(self, vector): @@ -234,7 +234,7 @@ def inner(self, vector): """ self._check_type(vector) total_prod = 0. - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): total_prod += self._vectors[i].inner(vector._vectors[i]) return total_prod @@ -248,7 +248,7 @@ def exp(self, vector): vector : CompositeVector """ self._check_type(vector) - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): self._vectors[i].exp(vector) def log(self, vector): @@ -261,7 +261,7 @@ def log(self, vector): vector : CompositeVector """ self._check_type(vector) - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): self._vectors[i].log(vector) def pow(self, power): @@ -272,7 +272,7 @@ def pow(self, power): ---------- power : float """ - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): self._vectors[i].pow(power) @property @@ -301,7 +301,7 @@ def infty(self): float : Infinity norm. """ norms = [] - for i in xrange(len(self._vectors)): + for i in range(len(self._vectors)): norms.append(self._vectors[i].infty) return max(norms) diff --git a/src/kona/optimizer.py b/src/kona/optimizer.py index 4d39c25..e6c4943 100644 --- a/src/kona/optimizer.py +++ b/src/kona/optimizer.py @@ -61,7 +61,7 @@ def __init__(self, solver, algorithm, optns=None): def _process_options(self, optns): # this is a recursive dictionary merge function def update(d, u): - for k, v in u.iteritems(): + for k, v in u.items(): if isinstance(v, collections.Mapping): r = update(d.get(k, {}), v) d[k] = r diff --git a/src/kona/options.py b/src/kona/options.py index 17e5de4..b69ba8d 100644 --- a/src/kona/options.py +++ b/src/kona/options.py @@ -30,7 +30,7 @@ def get_opt(optns, default, *keys): return val def print_dict(obj, pre='', out_file=sys.stdout): - for k, v in obj.items(): + for k, v in list(obj.items()): if hasattr(v, '__iter__'): out_file.write('%s%s : {\n'%(pre, k)) print_dict(v, pre='%s '%pre, out_file=out_file) diff --git a/src/kona/test/test_FLECS_solver.py b/src/kona/test/test_FLECS_solver.py index dea8f90..45d0b36 100644 --- a/src/kona/test/test_FLECS_solver.py +++ b/src/kona/test/test_FLECS_solver.py @@ -95,8 +95,8 @@ def test_radius_inactive_with_large_mu(self): total_data[3] = self.x.dual.base.data[:] diff = abs(total_data - expected) diff = max(diff) - print diff - print self.krylov.trust_active + print(diff) + print(self.krylov.trust_active) self.assertTrue(diff <= 1.e-3 and not self.krylov.trust_active) def test_radius_active(self): @@ -110,9 +110,9 @@ def test_radius_active(self): # compare actual result to expected exp_norm = self.krylov.radius actual_norm = self.x.primal.norm2 - print exp_norm - print actual_norm - print self.krylov.trust_active + print(exp_norm) + print(actual_norm) + print(self.krylov.trust_active) self.assertTrue( (exp_norm - actual_norm) <= 1e-1 and self.krylov.trust_active) diff --git a/src/kona/test/test_GCROT_solver.py b/src/kona/test/test_GCROT_solver.py index e99452d..7ade5c0 100644 --- a/src/kona/test/test_GCROT_solver.py +++ b/src/kona/test/test_GCROT_solver.py @@ -159,7 +159,7 @@ def test_multiple_recycle(self): self.krylov.max_krylov = 100 # solve the same problem 20 times - for k in xrange(20): + for k in range(20): self.x.equals(0.0) self.krylov.solve(self.mat_vec, self.b, self.x, self.precond.product) if k % 5 == 0: diff --git a/src/kona/test/test_approx_schur.py b/src/kona/test/test_approx_schur.py index d71e95a..71105e6 100644 --- a/src/kona/test/test_approx_schur.py +++ b/src/kona/test/test_approx_schur.py @@ -89,31 +89,31 @@ def solve(self): v_vec = self._generate_vector() self.dCdX.linearize(at_design, at_state) self.precond.linearize(at_design, at_state) - for i in xrange(20): + for i in range(20): u_vec.equals(0.) v_vec.equals(0.) u_vec.primal.base.data[i] = 1.0 self.exact_product(u_vec, v_vec) u_vec.equals(0.) self.precond.product(v_vec, u_vec) - for j in xrange(20): + for j in range(20): if j == i: self.assertAlmostEqual(u_vec.primal.base.data[j], 1., places=10) else: self.assertAlmostEqual(u_vec.primal.base.data[j], 0., places=10) - for j in xrange(10): + for j in range(10): self.assertAlmostEqual(u_vec.eq.base.data[j], 0., places=10) - for i in xrange(10): + for i in range(10): u_vec.equals(0.) v_vec.equals(0.) u_vec.eq.base.data[i] = 1.0 self.exact_product(u_vec, v_vec) u_vec.equals(0.) self.precond.product(v_vec, u_vec) - for j in xrange(20): + for j in range(20): self.assertAlmostEqual(u_vec.primal.base.data[j], 0., places=10) - for j in xrange(10): + for j in range(10): if j == i: self.assertAlmostEqual(u_vec.eq.base.data[j], 1., places=10) else: diff --git a/src/kona/test/test_aug_lag_merit.py b/src/kona/test/test_aug_lag_merit.py index 8621cd7..de6f3e0 100644 --- a/src/kona/test/test_aug_lag_merit.py +++ b/src/kona/test/test_aug_lag_merit.py @@ -60,7 +60,7 @@ def setUp(self): def test_init_func(self): '''AugmentedLagrangian merit initialization''' self.merit.reset(self.kkt_start, self.u_start, self.search_dir, self.mu) - print abs(self.merit.func_val-self.merit_val_init) + print(abs(self.merit.func_val-self.merit_val_init)) self.assertTrue(abs(self.merit.func_val-self.merit_val_init) <= 1e-10) def test_eval_func(self): diff --git a/src/kona/test/test_constraint_jacobian.py b/src/kona/test/test_constraint_jacobian.py index 2b2756f..155d003 100644 --- a/src/kona/test/test_constraint_jacobian.py +++ b/src/kona/test/test_constraint_jacobian.py @@ -81,14 +81,14 @@ def test_constrained_product(self): self.A.linearize(X.primal.design, state) self.A.product(in_vec.primal.design, out_vec.dual) - print '-----------------------------' - print 'Constraint Hessian' - print '-----------------------------' - print 'FD product:' - print dLdX_pert.dual.base.data - print 'Analytical product:' - print out_vec.dual.base.data - print '-----------------------------' + print('-----------------------------') + print('Constraint Hessian') + print('-----------------------------') + print('FD product:') + print(dLdX_pert.dual.base.data) + print('Analytical product:') + print(out_vec.dual.base.data) + print('-----------------------------') dLdX.equals_ax_p_by(1.0, dLdX_pert, -1.0, out_vec) diff_norm = dLdX.dual.norm2 diff --git a/src/kona/test/test_dual_vectors.py b/src/kona/test/test_dual_vectors.py index 3c1b446..b813ebf 100644 --- a/src/kona/test/test_dual_vectors.py +++ b/src/kona/test/test_dual_vectors.py @@ -62,12 +62,12 @@ def test_equals_mangasarian(self): self.dv.equals_constraints(at_design, at_state) self.mult.equals(-9) self.mang.equals_mangasarian(self.dv, self.mult) - for i in xrange(len(self.mang.base.data)): + for i in range(len(self.mang.base.data)): self.assertEqual(self.mang.base.data[i], -10) # linear Mangasarian # self.assertEqual(self.mang.base.data[i], -1730) # cubic Mangasarian self.mult.equals(-11) self.mang.equals_mangasarian(self.dv, self.mult) - for i in xrange(len(self.mang.base.data)): + for i in range(len(self.mang.base.data)): self.assertEqual(self.mang.base.data[i], -11) # linear Mangasarian # self.assertEqual(self.mang.base.data[i], -2332) # cubic Mangasarian diff --git a/src/kona/test/test_equality_composite_step.py b/src/kona/test/test_equality_composite_step.py index 5bcdb8c..e76528a 100644 --- a/src/kona/test/test_equality_composite_step.py +++ b/src/kona/test/test_equality_composite_step.py @@ -58,7 +58,7 @@ def test_with_simple_constrained(self): optimizer = Optimizer(solver, algorithm, optns) optimizer.solve() - print solver.curr_design + print(solver.curr_design) expected = -1.*numpy.ones(solver.num_design) diff = abs(solver.curr_design - expected) diff --git a/src/kona/test/test_equality_predictor_corrector.py b/src/kona/test/test_equality_predictor_corrector.py index d9248c7..5d730a8 100644 --- a/src/kona/test/test_equality_predictor_corrector.py +++ b/src/kona/test/test_equality_predictor_corrector.py @@ -48,7 +48,7 @@ def test_with_simple_constrained(self): optimizer = Optimizer(solver, algorithm, optns) optimizer.solve() - print solver.curr_design + print(solver.curr_design) expected = -1.*np.ones(solver.num_design) diff = abs(solver.curr_design - expected) diff --git a/src/kona/test/test_krylov_util.py b/src/kona/test/test_krylov_util.py index 702f1ea..8b1cb55 100644 --- a/src/kona/test/test_krylov_util.py +++ b/src/kona/test/test_krylov_util.py @@ -151,7 +151,7 @@ def rev_mat_vec(in_vec, out_vec): p_work = df.generate() Q = [] P = [] - for i in xrange(lanczos_size): + for i in range(lanczos_size): Q.append(pf.generate()) P.append(df.generate()) Q.append(pf.generate()) diff --git a/src/kona/test/test_lagrangian_hessian.py b/src/kona/test/test_lagrangian_hessian.py index 034b119..32787ca 100644 --- a/src/kona/test/test_lagrangian_hessian.py +++ b/src/kona/test/test_lagrangian_hessian.py @@ -78,14 +78,14 @@ def test_constrained_product(self): self.W.linearize(X, state, adjoint) self.W.multiply_W(in_vec.primal, out_vec.primal) - print '-----------------------------' - print 'Constrained Hessian' - print '-----------------------------' - print 'FD product:' - print dLdX_pert.primal.base.data - print 'Analytical product:' - print out_vec.primal.base.data - print '-----------------------------' + print('-----------------------------') + print('Constrained Hessian') + print('-----------------------------') + print('FD product:') + print(dLdX_pert.primal.base.data) + print('Analytical product:') + print(out_vec.primal.base.data) + print('-----------------------------') dLdX.equals_ax_p_by(1.0, dLdX_pert, -1.0, out_vec) diff_norm = dLdX.primal.norm2 diff --git a/src/kona/test/test_low_rank_svd.py b/src/kona/test/test_low_rank_svd.py index 4855db1..ec34b3e 100644 --- a/src/kona/test/test_low_rank_svd.py +++ b/src/kona/test/test_low_rank_svd.py @@ -65,32 +65,32 @@ def rev_mat_vec(in_vec, out_vec): self.A.product(in_vec.primal.design, out_vec_exact.dual) self.svd.approx_fwd_prod(in_vec.primal.design, out_vec_approx.dual) - print 'Constraint Jacobian test:' + print('Constraint Jacobian test:') - print 'Exact fwd product =', out_vec_exact.dual.norm2 - print 'Approx fwd product =', out_vec_approx.dual.norm2 + print('Exact fwd product =', out_vec_exact.dual.norm2) + print('Approx fwd product =', out_vec_approx.dual.norm2) out_vec_approx.dual.minus(out_vec_exact.dual) abs_error = out_vec_approx.dual.norm2 rel_error = abs_error/out_vec_exact.dual.norm2 - print 'Abs error norm =', abs_error - print 'Rel error norm =', rel_error + print('Abs error norm =', abs_error) + print('Rel error norm =', rel_error) self.assertTrue(rel_error <= 1e-8) self.A.T.product(in_vec.dual, out_vec_exact.primal.design) self.svd.approx_rev_prod(in_vec.dual, out_vec_approx.primal.design) - print 'Exact rev product =', out_vec_exact.primal.design.norm2 - print 'Approx fwd product =', out_vec_approx.primal.design.norm2 + print('Exact rev product =', out_vec_exact.primal.design.norm2) + print('Approx fwd product =', out_vec_approx.primal.design.norm2) out_vec_approx.primal.design.minus(out_vec_exact.primal.design) abs_error = out_vec_approx.primal.design.norm2 rel_error = abs_error/out_vec_exact.primal.design.norm2 - print 'Abs error norm =', abs_error - print 'Rel error norm =', rel_error + print('Abs error norm =', abs_error) + print('Rel error norm =', rel_error) self.assertTrue(rel_error <= 1e-8) @@ -140,16 +140,16 @@ def mat_vec(in_vec, out_vec): self.svd.approx_fwd_prod( in_vec.primal.design, out_vec_approx.primal.design) - print 'Hessian test:' - print 'Exact product =', out_vec_exact.primal.design.norm2 - print 'Approx product =', out_vec_approx.primal.design.norm2 + print('Hessian test:') + print('Exact product =', out_vec_exact.primal.design.norm2) + print('Approx product =', out_vec_approx.primal.design.norm2) out_vec_approx.primal.design.minus(out_vec_exact.primal.design) abs_error = out_vec_approx.primal.design.norm2 rel_error = abs_error/out_vec_exact.primal.design.norm2 - print 'Abs error norm =', abs_error - print 'Rel error norm =', rel_error + print('Abs error norm =', abs_error) + print('Rel error norm =', rel_error) self.assertTrue(rel_error <= 1e-8) @@ -185,12 +185,12 @@ def rev_mat_vec(in_vec, out_vec): # loop over and check each column in the approximate Schur self.svd.linearize() - for i in xrange(km.neq): + for i in range(km.neq): in_vec.base.data = np.zeros_like(in_vec.base.data) in_vec.base.data[i] = 1.0 out_vec.base.data = np.zeros_like(out_vec.base.data) self.svd.inv_schur_prod(in_vec, out_vec) - for j in xrange(km.neq): + for j in range(km.neq): self.assertAlmostEqual(Sinv[j,i]/out_vec.base.data[j], 1.0, places=9) def test_inv_schur_prod_ineq(self): @@ -225,12 +225,12 @@ def rev_mat_vec(in_vec, out_vec): # loop over and check each column in the approximate Schur self.svd.linearize() - for i in xrange(km.nineq): + for i in range(km.nineq): in_vec.base.data = np.zeros_like(in_vec.base.data) in_vec.base.data[i] = 1.0 out_vec.base.data = np.zeros_like(out_vec.base.data) self.svd.inv_schur_prod(in_vec, out_vec) - for j in xrange(km.nineq): + for j in range(km.nineq): self.assertAlmostEqual(Sinv[j,i]/out_vec.base.data[j], 1.0, places=9) def test_inv_schur_prod_comp(self): @@ -273,23 +273,23 @@ def rev_mat_vec(in_vec, out_vec): # loop over and check each column in the approximate Schur self.svd.linearize() - for i in xrange(km.neq): + for i in range(km.neq): in_vec.equals(0.) in_vec.eq.base.data[i] = 1.0 out_vec.equals(0.) self.svd.inv_schur_prod(in_vec, out_vec) - for j in xrange(km.neq): + for j in range(km.neq): self.assertAlmostEqual(Sinv[j,i]/out_vec.eq.base.data[j], 1.0, places=9) - for j in xrange(km.nineq): + for j in range(km.nineq): self.assertAlmostEqual(Sinv[km.neq+j,i]/out_vec.ineq.base.data[j], 1.0, places=9) - for i in xrange(km.nineq): + for i in range(km.nineq): in_vec.equals(0.) in_vec.ineq.base.data[i] = 1.0 out_vec.equals(0.) self.svd.inv_schur_prod(in_vec, out_vec) - for j in xrange(km.neq): + for j in range(km.neq): self.assertAlmostEqual(Sinv[j,km.neq+i]/out_vec.eq.base.data[j], 1.0, places=9) - for j in xrange(km.nineq): + for j in range(km.nineq): self.assertAlmostEqual(Sinv[km.neq+j,km.neq+i]/out_vec.ineq.base.data[j], 1.0, places=9) diff --git a/src/kona/test/test_primal_dual_vector.py b/src/kona/test/test_primal_dual_vector.py index a4339db..c7d275c 100644 --- a/src/kona/test/test_primal_dual_vector.py +++ b/src/kona/test/test_primal_dual_vector.py @@ -560,7 +560,7 @@ def test_get_base_data_case3(self): def test_get_base_data_case4(self): # case 4 has no constraints A = np.zeros((10,1)) - print A.shape + print(A.shape) self.at_pd4.get_base_data(A[:,0]) B = np.ones_like(A) B[0:10] *= 10 diff --git a/src/kona/test/test_reduced_kkt_matrix.py b/src/kona/test/test_reduced_kkt_matrix.py index 0005d48..7557c08 100644 --- a/src/kona/test/test_reduced_kkt_matrix.py +++ b/src/kona/test/test_reduced_kkt_matrix.py @@ -83,16 +83,16 @@ def test_equality_constrained_product(self): self.KKT_matrix.linearize(X, state, adjoint) self.KKT_matrix.product(in_vec, out_vec) - print '----------------------' - print 'Equality Constraints' - print '----------------------' - print 'FD product:' - print dLdX_pert.primal.base.data - print dLdX_pert.dual.base.data - print 'Analytical product:' - print out_vec.primal.base.data - print out_vec.dual.base.data - print '----------------------' + print('----------------------') + print('Equality Constraints') + print('----------------------') + print('FD product:') + print(dLdX_pert.primal.base.data) + print(dLdX_pert.dual.base.data) + print('Analytical product:') + print(out_vec.primal.base.data) + print(out_vec.dual.base.data) + print('----------------------') dLdX.equals_ax_p_by(1.0, dLdX_pert, -1.0, out_vec) diff_norm = dLdX.norm2 diff --git a/src/kona/test/test_verifier.py b/src/kona/test/test_verifier.py index 8ee6f39..2b3204a 100644 --- a/src/kona/test/test_verifier.py +++ b/src/kona/test/test_verifier.py @@ -29,7 +29,7 @@ def test_sellar(self): optimizer = Optimizer(solver, algorithm, optns) optimizer.solve() - self.failUnless('Output inspected by hand...') + self.assertTrue('Output inspected by hand...') if __name__ == "__main__": unittest.main() diff --git a/src/kona/user/__init__.py b/src/kona/user/__init__.py index b4f1cf8..383b9fa 100644 --- a/src/kona/user/__init__.py +++ b/src/kona/user/__init__.py @@ -1,3 +1,3 @@ -from base_vectors import BaseVector -from user_solver import UserSolver -from user_solver import UserSolverIDF +from .base_vectors import BaseVector +from .user_solver import UserSolver +from .user_solver import UserSolverIDF