Skip to content

Add active set recovery to residuum controller #109

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Add active set recovery to residuum controller #109

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e9ddf13
Add active set recovery to residuum controller
chrhansk Aug 8, 2024
db8edc7
Add recovery step test
chrhansk Aug 8, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 19 additions & 0 deletions pygradflow/implicit_func.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,25 @@ def active_set_at_point(self, p):
ub = problem.var_ub
return super().compute_active_set_box(p, lb, ub)

def flow_rhs(self, iterate: Iterate, rho: float) -> np.ndarray:

lb = self.problem.var_lb
ub = self.problem.var_ub
x = iterate.x

xflow = -iterate.aug_lag_deriv_x(rho)

# Project flow onto tnagent cone of box constraints at given x
at_lb = np.where(x < lb - 1e-8)[0]
at_ub = np.where(x > ub + 1e-8)[0]

xflow[at_lb] = np.maximum(0.0, xflow[at_lb])
xflow[at_ub] = np.minimum(0.0, xflow[at_ub])

yflow = iterate.aug_lag_deriv_y()

return np.concatenate([xflow, yflow])

def projection_initial(self, iterate: Iterate, rho: float, tau=None):
x_0 = self.orig_iterate.x
dt = self.dt
Expand Down
88 changes: 75 additions & 13 deletions pygradflow/step/residuum_ratio_control.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
from pygradflow.params import Params
from pygradflow.problem import Problem
from pygradflow.step.newton_control import NewtonController
from pygradflow.step.solver.step_solver import StepResult
from pygradflow.step.step_control import StepControlResult


Expand All @@ -29,7 +30,8 @@ def step(self, iterate, rho, dt, display, timer):
mid_step = next(next_steps)
mid_iterate = mid_step.iterate

mid_norm = np.linalg.norm(func.value_at(mid_iterate, rho))
mid_value = func.value_at(mid_iterate, rho)
mid_norm = np.linalg.norm(mid_value)

self.display_step(0, mid_step)

Expand All @@ -43,21 +45,81 @@ def step(self, iterate, rho, dt, display, timer):
theta = mid_norm / orig_norm
accepted = theta <= params.theta_max

if accepted:
lamb_mod = self.controller.update(theta)
lamb_n = max(params.lamb_min, lamb / lamb_mod)
else:
lamb_n = lamb * params.lamb_inc
if self.controller.error_sum > 0.0:
self.controller.reset()

logger.debug(
"StepController: theta: %e, accepted: %s, lamb: %e, lamb_n: %e",
"StepController: theta: %e, accepted: %s, lamb: %e",
theta,
accepted,
lamb,
lamb_n,
)

self.lamb = lamb_n
return StepControlResult.from_step_result(mid_step, lamb_n, accepted)
if accepted:
lamb_mod = self.controller.update(theta)
lamb_n = max(params.lamb_min, lamb / lamb_mod)
self.lamb = lamb_n
return StepControlResult.from_step_result(mid_step, lamb_n, accepted)

# Step would be rejected
lamb_n = lamb * params.lamb_inc

if self.controller.error_sum > 0.0:
self.controller.reset()

# Recovery starts here

active_set = func.compute_active_set(iterate, rho)
mid_active_set = func.compute_active_set(mid_iterate, rho)

if (active_set == mid_active_set).all():
return StepControlResult.from_step_result(mid_step, lamb_n, accepted)

mid_primal_value = mid_value[: self.problem.num_vars]

mid_primal_norm = np.linalg.norm(mid_primal_value)

mid_unchanged_primal_value = mid_primal_value[mid_active_set == active_set]
mid_unchanged_primal_norm = np.linalg.norm(mid_unchanged_primal_value)

norm_factor = mid_unchanged_primal_norm / mid_primal_norm

if norm_factor >= 1e-6:
self.lamb = lamb_n
return StepControlResult.from_step_result(mid_step, lamb_n, accepted)

dir = func.flow_rhs(mid_iterate, rho)

# Track back until first active step change...

primal_dir = dir[: self.problem.num_vars]
dual_dir = dir[self.problem.num_vars :]

pos_dir = primal_dir > 0.0
neg_dir = primal_dir < 0.0

lb = problem.var_lb
ub = problem.var_ub

pos_ratio = (ub[pos_dir] - mid_primal_value[pos_dir]) / primal_dir[pos_dir]
neg_ratio = (lb[neg_dir] - mid_primal_value[neg_dir]) / primal_dir[neg_dir]

recovery_dt = np.minimum(np.min(pos_ratio), np.min(neg_ratio))
recovery_dt = np.minimum(dt, recovery_dt)

recovery_dx = recovery_dt * primal_dir
recovery_dy = recovery_dt * dual_dir

recovery_step = StepResult(iterate, recovery_dx, recovery_dy, active_set)

recovery_func = ImplicitFunc(problem, iterate, recovery_dt)

recovery_value = recovery_func.value_at(recovery_step.iterate, rho)
recovery_norm = np.linalg.norm(recovery_value)

recovery_theta = recovery_norm / orig_norm
accepted = recovery_theta <= params.theta_max

if accepted:
lamb_n = lamb
else:
lamb_n = 2.0 * lamb

return StepControlResult.from_step_result(recovery_step, lamb_n, accepted=accepted)
48 changes: 48 additions & 0 deletions tests/pygradflow/linprog.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
import numpy as np
import scipy as sp

from pygradflow.problem import Problem


class LinearProgram(Problem):
"""
Linear program: minimize c^T x subject to Ax <= b
"""

def __init__(self, A, b, c):
(num_cons, num_vars) = A.shape

var_lb = np.full(shape=(num_vars,), fill_value=-np.inf)
var_ub = np.full(shape=(num_vars,), fill_value=np.inf)

cons_lb = np.full(shape=(num_cons,), fill_value=-np.inf)
cons_ub = b

self.c = c
self.b = b
self.A = sp.sparse.coo_matrix(A)

super().__init__(var_lb, var_ub, cons_lb=cons_lb, cons_ub=cons_ub)

def obj(self, x):
return np.dot(self.c, x)

def obj_grad(self, x):
return self.c

def cons(self, x):
return self.A.dot(x)

def cons_jac(self, x):
return self.A

def lag_hess(self, x, y):
return sp.sparse.coo_matrix(
np.zeros(
(
self.num_vars,
self.num_vars,
),
dtype=float,
)
)
154 changes: 154 additions & 0 deletions tests/pygradflow/test_solve_linear.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
import numpy as np
import pytest

from pygradflow.integration.integration_solver import IntegrationSolver
from pygradflow.params import Params, PenaltyUpdate, StepControlType
from pygradflow.solver import Solver

from .linprog import LinearProgram


@pytest.fixture
def linprog():
A = np.array(
[
[9.51056516e-01, 6.90983006e-01],
[-3.63271264e-01, 1.11803399e00],
[-1.17557050e00, 2.22044605e-16],
[-3.63271264e-01, -1.11803399e00],
[9.51056516e-01, -6.90983006e-01],
[-2.17234844e-01, 1.07058388e00],
[-7.11984239e-02, 1.02313377e00],
[7.48379961e-02, 9.75683661e-01],
[2.20874416e-01, 9.28233552e-01],
[3.66910836e-01, 8.80783443e-01],
[5.12947256e-01, 8.33333333e-01],
[6.58983676e-01, 7.85883224e-01],
[8.05020096e-01, 7.38433115e-01],
[-1.08531503e00, 1.24225999e-01],
[-9.95059562e-01, 2.48451997e-01],
[-9.04804091e-01, 3.72677996e-01],
[-8.14548620e-01, 4.96903995e-01],
[-7.24293149e-01, 6.21129994e-01],
[-6.34037678e-01, 7.45355992e-01],
[-5.43782206e-01, 8.69581991e-01],
[-4.53526735e-01, 9.93807990e-01],
[-4.53526735e-01, -9.93807990e-01],
[-5.43782206e-01, -8.69581991e-01],
[-6.34037678e-01, -7.45355992e-01],
[-7.24293149e-01, -6.21129994e-01],
[-8.14548620e-01, -4.96903995e-01],
[-9.04804091e-01, -3.72677996e-01],
[-9.95059562e-01, -2.48451997e-01],
[-1.08531503e00, -1.24225999e-01],
[8.05020096e-01, -7.38433115e-01],
[6.58983676e-01, -7.85883224e-01],
[5.12947256e-01, -8.33333333e-01],
[3.66910836e-01, -8.80783443e-01],
[2.20874416e-01, -9.28233552e-01],
[7.48379961e-02, -9.75683661e-01],
[-7.11984239e-02, -1.02313377e00],
[-2.17234844e-01, -1.07058388e00],
[9.51056516e-01, 5.37431227e-01],
[9.51056516e-01, 3.83879448e-01],
[9.51056516e-01, 2.30327669e-01],
[9.51056516e-01, 7.67758895e-02],
[9.51056516e-01, -7.67758895e-02],
[9.51056516e-01, -2.30327669e-01],
[9.51056516e-01, -3.83879448e-01],
[9.51056516e-01, -5.37431227e-01],
]
)

b = np.array(
[
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
0.95105652,
]
)

c = np.array([0.0, 0.1])

return LinearProgram(A, b, c)


@pytest.fixture
def x0():
return np.array([0.0, 0.5])


@pytest.fixture
def y0():
return 0.0


def test_solve_with_recovery(linprog, x0, y0):
params = Params(
collect_path=True,
rho=1e-1,
step_control_type=StepControlType.ResiduumRatio,
penalty_update=PenaltyUpdate.Constant,
)

solver = Solver(linprog, params)

result = solver.solve(x0, y0)

import matplotlib.pyplot as plt

primal_path = result.primal_path

path_x = primal_path[0, :]
path_y = primal_path[1, :]

# plt.plot(path_x, path_y, "o-")
# plt.show()

# import pdb

# pdb.set_trace()

assert result.success
Loading