From 572b2e037e3975af0e2938bf8e6ddd94829975ee Mon Sep 17 00:00:00 2001 From: Thomas Keefe Date: Thu, 30 Sep 2021 17:02:58 -0400 Subject: [PATCH 1/6] Add questions about L1 ball proj --- ya_glm/opt/constraint/convex.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ya_glm/opt/constraint/convex.py b/ya_glm/opt/constraint/convex.py index 85deb64..2228ae4 100644 --- a/ya_glm/opt/constraint/convex.py +++ b/ya_glm/opt/constraint/convex.py @@ -58,8 +58,9 @@ def is_proximable(self): # See https://gist.github.com/mblondel/6f3b7aaad90606b98f71 # for more algorithms. def project_simplex(v, z=1): - if np.sum(v) <= z: - return v + # z is what the entries need to add up to, e.g. z=1 for probability simplex + if np.sum(v) <= z: # don't we want the simplex to mean sum == z not sum <= z? + return v # also this doesn't work when v has, say, all negative entries n_features = v.shape[0] u = np.sort(v)[::-1] From b47bfeb24c9f58a251d186ae980741b2d1721546 Mon Sep 17 00:00:00 2001 From: Thomas Keefe Date: Thu, 30 Sep 2021 17:50:25 -0400 Subject: [PATCH 2/6] Add L2Ball Constraint and tests --- ya_glm/opt/constraint/convex.py | 13 ++++++++++++ ya_glm/opt/constraint/tests.py | 35 +++++++++++++++++++++++++++++++++ 2 files changed, 48 insertions(+) create mode 100644 ya_glm/opt/constraint/tests.py diff --git a/ya_glm/opt/constraint/convex.py b/ya_glm/opt/constraint/convex.py index 2228ae4..82601c0 100644 --- a/ya_glm/opt/constraint/convex.py +++ b/ya_glm/opt/constraint/convex.py @@ -75,3 +75,16 @@ def project_simplex(v, z=1): def project_l1_ball(v, z=1): return np.sign(v) * project_simplex(np.abs(v), z) + +class L2Ball(Constraint): + + def __init__(self, mult=1): + assert mult > 0 + self.mult = mult + + def _prox(self, x, step=1): + return x / np.max([np.linalg.norm(x)/self.mult, 1]) + + @property + def is_proximable(self): + return True diff --git a/ya_glm/opt/constraint/tests.py b/ya_glm/opt/constraint/tests.py new file mode 100644 index 0000000..5dbc92b --- /dev/null +++ b/ya_glm/opt/constraint/tests.py @@ -0,0 +1,35 @@ +import numpy as np +from unittest import TestCase +from ya_glm.opt.constraint import convex + +class TestProjectionsOnConstraints(TestCase): + def setUp(self): + pass + + def assert_arrays_close(self, test_array, ref_array): + "Custom assertion for arrays being almost equal" + try: + np.testing.assert_allclose(test_array, ref_array) + except AssertionError: + self.fail() + + def test_Positive(self): + cons = convex.Positive() + v = np.array([-1, 0, 2, 3, -2]) + self.assert_arrays_close(cons.prox(v), [0, 0, 2, 3, 0]) + self.assertEqual(cons.prox(-2), 0) + + + def test_L2Ball(self): + cons1 = convex.L2Ball(1) + self.assert_arrays_close(cons1.prox([0,0,0]), [0,0,0]) + self.assert_arrays_close(cons1.prox([1,0,0]), [1,0,0]) + self.assert_arrays_close(cons1.prox([0.5,0,0]), [0.5,0,0]) + self.assert_arrays_close(cons1.prox([1,1,1]), np.array([1,1,1])/np.sqrt(3)) + self.assert_arrays_close(cons1.prox([1,-1,1]), np.array([1,-1,1])/np.sqrt(3)) + + cons4 = convex.L2Ball(4) + self.assert_arrays_close(cons4.prox([0,0,0]), [0,0,0]) + self.assert_arrays_close(cons4.prox([1,0,0]), [1,0,0]) + self.assert_arrays_close(cons4.prox([0.5,0,0]), [0.5,0,0]) + self.assert_arrays_close(cons4.prox([-4,3,0]), np.array([-4,3,0])/(5/4)) From 7aad227843486b5403508295449b38f699e15b15 Mon Sep 17 00:00:00 2001 From: Thomas Keefe Date: Sat, 2 Oct 2021 14:49:45 -0400 Subject: [PATCH 3/6] Add Isotonic constraint and tests --- ya_glm/opt/constraint/convex.py | 17 ++++++++++++++++- ya_glm/opt/constraint/tests.py | 21 +++++++++++++++++++++ 2 files changed, 37 insertions(+), 1 deletion(-) diff --git a/ya_glm/opt/constraint/convex.py b/ya_glm/opt/constraint/convex.py index 82601c0..ff08dc9 100644 --- a/ya_glm/opt/constraint/convex.py +++ b/ya_glm/opt/constraint/convex.py @@ -1,5 +1,5 @@ import numpy as np - +from sklearn.isotonic import isotonic_regression from ya_glm.opt.base import Func @@ -88,3 +88,18 @@ def _prox(self, x, step=1): @property def is_proximable(self): return True + +class Isotonic(Constraint): + """Constraint for x1 <= ... <= xn or + x1 >= ... >= xn """ + def __init__(self, increasing=True): + self.increasing = increasing + + def _prox(self, x, step=1): + # computes the projection of x onto the monotone cone + # using the PAVA algorithm + return isotonic_regression(x, increasing=self.increasing) + + @property + def is_proximable(self): + return True diff --git a/ya_glm/opt/constraint/tests.py b/ya_glm/opt/constraint/tests.py index 5dbc92b..4ac39f2 100644 --- a/ya_glm/opt/constraint/tests.py +++ b/ya_glm/opt/constraint/tests.py @@ -33,3 +33,24 @@ def test_L2Ball(self): self.assert_arrays_close(cons4.prox([1,0,0]), [1,0,0]) self.assert_arrays_close(cons4.prox([0.5,0,0]), [0.5,0,0]) self.assert_arrays_close(cons4.prox([-4,3,0]), np.array([-4,3,0])/(5/4)) + + def test_Isotonic(self): + cons = convex.Isotonic(increasing=True) + for v in [ + np.arange(5), + np.array([-1, 0, 2, 3, -2]), + np.array([-1, 3, 0, 3, 2]) + ]: + result = cons.prox(v) + lags = result[1:] - result[:-1] + self.assertTrue((lags >= 0).all()) + + cons = convex.Isotonic(increasing=False) + for v in [ + np.arange(5), + np.array([-1, 0, 2, 3, -2]), + np.array([-1, 3, 0, 3, 2]) + ]: + result = cons.prox(v) + lags = result[1:] - result[:-1] + self.assertTrue((lags <= 0).all()) From 308ff285b6a59b8c20664909583915169682749e Mon Sep 17 00:00:00 2001 From: Thomas Keefe Date: Sat, 2 Oct 2021 14:53:46 -0400 Subject: [PATCH 4/6] Add comment on isotonic regression --- ya_glm/opt/constraint/convex.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/ya_glm/opt/constraint/convex.py b/ya_glm/opt/constraint/convex.py index ff08dc9..587fb11 100644 --- a/ya_glm/opt/constraint/convex.py +++ b/ya_glm/opt/constraint/convex.py @@ -92,6 +92,10 @@ def is_proximable(self): class Isotonic(Constraint): """Constraint for x1 <= ... <= xn or x1 >= ... >= xn """ + # TODO: allow for general isotonic regression + # where the order relations are a simple directed + # graph. For an algorithm see Nemeth and Nemeth, "How to project onto an + # isotone projection cone", JLLA 2010 def __init__(self, increasing=True): self.increasing = increasing From 2e8cbf6823e1b20cf3dc517dfac45579c711b5e6 Mon Sep 17 00:00:00 2001 From: Thomas Keefe Date: Sat, 2 Oct 2021 15:21:28 -0400 Subject: [PATCH 5/6] Add LinearEquality constraint and test --- ya_glm/opt/constraint/convex.py | 19 ++++++++++++++++++- ya_glm/opt/constraint/tests.py | 6 ++++++ 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/ya_glm/opt/constraint/convex.py b/ya_glm/opt/constraint/convex.py index 587fb11..ff01b61 100644 --- a/ya_glm/opt/constraint/convex.py +++ b/ya_glm/opt/constraint/convex.py @@ -1,4 +1,5 @@ import numpy as np +from numpy.linalg import norm from sklearn.isotonic import isotonic_regression from ya_glm.opt.base import Func @@ -25,6 +26,22 @@ def _prox(self, x, step=1): def is_proximable(self): return True +class LinearEquality(Constraint): + # credited to PyUNLocBoX + # https://github.com/epfl-lts2/pyunlocbox/ + def __init__(self, A, b): + self.A = A + self.b = b + self.pinvA = np.linalg.pinv(A) + + def _prox(self, x, step=1): + residue = self.A@x - self.b + sol = x - self.pinvA @ residue + return sol + + @property + def is_proximable(self): + return True class Simplex(Constraint): @@ -83,7 +100,7 @@ def __init__(self, mult=1): self.mult = mult def _prox(self, x, step=1): - return x / np.max([np.linalg.norm(x)/self.mult, 1]) + return x / np.max([norm(x)/self.mult, 1]) @property def is_proximable(self): diff --git a/ya_glm/opt/constraint/tests.py b/ya_glm/opt/constraint/tests.py index 4ac39f2..1f8db2f 100644 --- a/ya_glm/opt/constraint/tests.py +++ b/ya_glm/opt/constraint/tests.py @@ -19,6 +19,12 @@ def test_Positive(self): self.assert_arrays_close(cons.prox(v), [0, 0, 2, 3, 0]) self.assertEqual(cons.prox(-2), 0) + def test_LinearEquality(self): + A = np.identity(2) + b = np.array([1,1]) + cons = convex.LinearEquality(A, b) + proj = cons.prox(b) # the proj of b should just be b + self.assert_arrays_close(proj, b) def test_L2Ball(self): cons1 = convex.L2Ball(1) From bb6e60812f9b5912cd5e01e7c66e4493ab935868 Mon Sep 17 00:00:00 2001 From: Thomas Keefe Date: Wed, 10 Nov 2021 19:27:46 -0500 Subject: [PATCH 6/6] Fix simplex projection and linear equality constraints --- yaglm/opt/constraint/convex.py | 13 ++++++++----- yaglm/opt/constraint/tests.py | 2 +- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/yaglm/opt/constraint/convex.py b/yaglm/opt/constraint/convex.py index ff01b61..9b23f87 100644 --- a/yaglm/opt/constraint/convex.py +++ b/yaglm/opt/constraint/convex.py @@ -1,7 +1,7 @@ import numpy as np from numpy.linalg import norm from sklearn.isotonic import isotonic_regression -from ya_glm.opt.base import Func +from yaglm.opt.base import Func class Constraint(Func): @@ -32,13 +32,19 @@ class LinearEquality(Constraint): def __init__(self, A, b): self.A = A self.b = b - self.pinvA = np.linalg.pinv(A) + self._pinvA = None # will be set on first request to .pinvA def _prox(self, x, step=1): residue = self.A@x - self.b sol = x - self.pinvA @ residue return sol + @property + def pinvA(self): + if self._pinvA is None: + self._pinvA = np.linalg.pinv(self.A) + return self._pinvA + @property def is_proximable(self): return True @@ -76,9 +82,6 @@ def is_proximable(self): # for more algorithms. def project_simplex(v, z=1): # z is what the entries need to add up to, e.g. z=1 for probability simplex - if np.sum(v) <= z: # don't we want the simplex to mean sum == z not sum <= z? - return v # also this doesn't work when v has, say, all negative entries - n_features = v.shape[0] u = np.sort(v)[::-1] cssv = np.cumsum(u) - z diff --git a/yaglm/opt/constraint/tests.py b/yaglm/opt/constraint/tests.py index 1f8db2f..a3099bc 100644 --- a/yaglm/opt/constraint/tests.py +++ b/yaglm/opt/constraint/tests.py @@ -1,6 +1,6 @@ import numpy as np from unittest import TestCase -from ya_glm.opt.constraint import convex +from yaglm.opt.constraint import convex class TestProjectionsOnConstraints(TestCase): def setUp(self):