Merge pull request #34 from epfl-lts2/new-projections

feat: add spsd and lineq projection operators

CHANGELOG.rst

@@ -9,6 +9,8 @@ and this project adheres to `Semantic Versioning <https://semver.org>`_.
 Unreleased
 ----------
 
+* New function: proj_linalg.
+* New function: proj_sdsp.
 * New function: proj_positive.
 * New function: structured_sparsity.
 * Continuous integration with Python 3.6, 3.7, 3.8, 3.9. Dropped 2.7, 3.4, 3.5.

pyunlocbox/functions.py

@@ -36,6 +36,8 @@
 
     proj_positive
     proj_b2
+    proj_lineq
+    proj_spsd
 
 **Miscellaneous**
 
@@ -803,24 +805,13 @@ class proj(func):
     See generic attributes descriptions of the
     :class:`pyunlocbox.functions.func` base class.
 
-    Parameters
-    ----------
-    epsilon : float, optional
-        The radius of the ball. Default is 1.
-    method : {'FISTA', 'ISTA'}, optional
-        The method used to solve the problem. It can be 'FISTA' or 'ISTA'.
-        Default is 'FISTA'.
-
     Notes
     -----
     * All indicator functions (projections) evaluate to zero by definition.
 
     """
-
-    def __init__(self, epsilon=1, method='FISTA', **kwargs):
+    def __init__(self, **kwargs):
         super(proj, self).__init__(**kwargs)
-        self.epsilon = epsilon
-        self.method = method
 
     def _eval(self, x):
         # Matlab version returns a small delta to avoid division by 0 when
@@ -835,9 +826,9 @@ class proj_positive(proj):
     r"""
     Projection on the positive octant (eval, prox).
 
-    This function is the indicator function :math:`i_S(z)` of the set S which
-    is zero if `z` is in the set and infinite otherwise. The set S is defined
-    by :math:`\left\{z \in \mathbb{R}^N \mid z \leq 0 \right\}`.
+    This function is the indicator function :math:`i_S(z)` of the set
+    :math:`S = \left\{z \in \mathbb{R}^N \mid z \leq 0 \right\}`
+    that is zero if :math:`z` is in the set and infinite otherwise.
 
     See generic attributes descriptions of the
     :class:`pyunlocbox.functions.proj` base class. Note that the constructor
@@ -867,19 +858,78 @@ def _prox(self, x, T):
         return np.clip(x, 0, np.inf)
 
 
+class proj_spsd(proj):
+    r"""
+    Projection on symmetric positive semi-definite matrices (eval, prox).
+
+    This function is the indicator function :math:`i_S(M)` of the set
+    :math:`S = \left\{M \in \mathbb{R}^{N \times N}
+    \mid M \succeq 0, M=M^T \right\}`
+    that is zero if :math:`M` is in the set and infinite otherwise.
+
+    See generic attributes descriptions of the
+    :class:`pyunlocbox.functions.proj` base class. Note that the constructor
+    takes keyword-only parameters.
+
+    Notes
+    -----
+    * The evaluation of this function is zero.
+
+    Examples
+    --------
+    >>> from pyunlocbox import functions
+    >>> f = functions.proj_spsd()
+    >>> A = np.array([[0, -1] , [-1, 1]])
+    >>> A = (A + A.T) / 2  # Symmetrize the matrix.
+    >>> np.linalg.eig(A)[0]
+    array([-0.61803399,  1.61803399])
+    >>> f.eval(A)
+    0
+    >>> Aproj = f.prox(A, 0)
+    >>> np.linalg.eig(Aproj)[0]
+    array([0.        , 1.61803399])
+
+    """
+    def __init__(self, **kwargs):
+        # Constructor takes keyword-only parameters to prevent user errors.
+        super(proj_spsd, self).__init__(**kwargs)
+
+    def _prox(self, x, T):
+        isreal = np.isreal(x).all()
+
+        # 1. make it symmetric.
+        sol = (x + np.conj(x.T)) / 2
+
+        # 2. make it semi-positive.
+        D, V = np.linalg.eig(sol)
+        D = np.real(D)
+        if isreal:
+            V = np.real(V)
+        D = np.clip(D, 0, np.inf)
+        sol = V @ np.diag(D) @ np.conj(V.T)
+        return sol
+
+
 class proj_b2(proj):
     r"""
     Projection on the L2-ball (eval, prox).
 
-    This function is the indicator function :math:`i_S(z)` of the set S which
-    is zero if `z` is in the set and infinite otherwise. The set S is defined
-    by :math:`\left\{z \in \mathbb{R}^N \mid \|A(z)-y\|_2 \leq \epsilon
-    \right\}`.
+    This function is the indicator function :math:`i_S(z)` of the set
+    :math:`S= \left\{z \in \mathbb{R}^N \mid \|Az-y\|_2 \leq \epsilon \right\}`
+    that is zero if :math:`z` is in the set and infinite otherwise.
 
     See generic attributes descriptions of the
     :class:`pyunlocbox.functions.proj` base class. Note that the constructor
     takes keyword-only parameters.
 
+    Parameters
+    ----------
+    epsilon : float, optional
+        The radius of the ball. Default is 1.
+    method : {'FISTA', 'ISTA'}, optional
+        The method used to solve the problem. It can be 'FISTA' or 'ISTA'.
+        Default is 'FISTA'.
+
     Notes
     -----
     * The `tol` parameter is defined as the tolerance for the projection on the
@@ -893,6 +943,10 @@ class proj_b2(proj):
       :math:`\|A(z)-y\|_2 \leq \epsilon`. It is thus a projection of the vector
       `x` onto an L2-ball of diameter `epsilon`.
 
+    See Also
+    --------
+    proj_lineq : use instead of ``epsilon=0``
+
     Examples
     --------
     >>> from pyunlocbox import functions
@@ -904,10 +958,11 @@ class proj_b2(proj):
     array([1.70710678, 1.70710678])
 
     """
-
-    def __init__(self, **kwargs):
+    def __init__(self, epsilon=1, method='FISTA', **kwargs):
         # Constructor takes keyword-only parameters to prevent user errors.
         super(proj_b2, self).__init__(**kwargs)
+        self.epsilon = epsilon
+        self.method = method
 
     def _prox(self, x, T):
 
@@ -993,6 +1048,78 @@ def _prox(self, x, T):
         return sol
 
 
+class proj_lineq(proj):
+    r"""
+    Projection on the plane satisfying the linear equality Az = y (eval, prox).
+
+    This function is the indicator function :math:`i_S(z)` of the set
+    :math:`S = \left\{z \in \mathbb{R}^N \mid Az = y \right\}`
+    that is zero if :math:`z` is in the set and infinite otherwise.
+
+    The proximal operator is
+    :math:`\operatorname{arg\,min}_z \| z - x \|_2 \text{ s.t. } Az = y`.
+
+    See generic attributes descriptions of the
+    :class:`pyunlocbox.functions.proj` base class. Note that the constructor
+    takes keyword-only parameters.
+
+    Notes
+    -----
+    * A parameter `pinvA`, the pseudo-inverse of `A`, must be provided if the
+      parameter `A` is provided as an operator/callable (not a matrix).
+    * The evaluation of this function is zero.
+
+    See Also
+    --------
+    proj_b2 : quadratic case
+
+    Examples
+    --------
+    >>> from pyunlocbox import functions
+    >>> import numpy as np
+    >>> x = np.array([0, 0])
+    >>> A = np.array([[1, 1]])
+    >>> pinvA = np.linalg.pinv(A)
+    >>> y = np.array([1])
+    >>> f = functions.proj_lineq(A=A, pinvA=pinvA, y=y)
+    >>> sol = f.prox(x, 0)
+    >>> sol
+    array([0.5, 0.5])
+    >>> np.abs(A.dot(sol) - y) < 1e-15
+    array([ True])
+
+    """
+    def __init__(self, A=None, pinvA=None, **kwargs):
+        # Constructor takes keyword-only parameters to prevent user errors.
+        super(proj_lineq, self).__init__(A=A, **kwargs)
+
+        if pinvA is None:
+            if A is None:
+                print("Are you sure about the parameters?" +
+                      "The projection will return y.")
+                self.pinvA = lambda x: x
+            else:
+                if callable(A):
+                    raise ValueError(
+                            "Provide A as a numpy array or provide pinvA.")
+                else:
+                    # Transform matrix form to operator form.
+                    self._pinvA = np.linalg.pinv(A)
+                    self.pinvA = lambda x: self._pinvA.dot(x)
+        else:
+            if callable(pinvA):
+                self.pinvA = pinvA
+            else:
+                self.pinvA = lambda x: pinvA.dot(x)
+
+    def _prox(self, x, T):
+        # Applying the projection formula.
+        # (for now, only the non scalable version)
+        residue = self.A(x) - self.y()
+        sol = x - self.pinvA(residue)
+        return sol
+
+
 class structured_sparsity(func):
     r"""
     Structured sparsity (eval, prox).

pyunlocbox/tests/test_functions.py

@@ -51,8 +51,10 @@ def assert_equivalent(param1, param2):
         assert_equivalent({'y': 3.2}, {'y': lambda: 3.2})
         assert_equivalent({'A': None}, {'A': np.identity(3)})
         A = np.array([[-4, 2, 5], [1, 3, -7], [2, -1, 0]])
+        pinvA = np.linalg.pinv(A)  # For proj_linalg.
         assert_equivalent({'A': A}, {'A': A, 'At': A.T})
-        assert_equivalent({'A': lambda x: A.dot(x)}, {'A': A, 'At': A})
+        assert_equivalent({'A': lambda x: A.dot(x), 'pinvA': pinvA},
+                          {'A': A, 'At': A})
 
     def test_dummy(self):
         """
@@ -95,16 +97,16 @@ def test_norm_l2(self):
         self.assertEqual(f.eval([4, 6]), 0)
         self.assertEqual(f.eval([5, -2]), 256 + 4)
         nptest.assert_allclose(f.grad([4, 6]), 0)
-#        nptest.assert_allclose(f.grad([5, -2]), [8, -64])
+        # nptest.assert_allclose(f.grad([5, -2]), [8, -64])
         nptest.assert_allclose(f.prox([4, 6], 1), [4, 6])
 
         f = functions.norm_l2(lambda_=2, y=np.fft.fft([2, 4]) / np.sqrt(2),
                               A=lambda x: np.fft.fft(x) / np.sqrt(x.size),
                               At=lambda x: np.fft.ifft(x) * np.sqrt(x.size))
-#        self.assertEqual(f.eval(np.fft.ifft([2, 4])*np.sqrt(2)), 0)
-#        self.assertEqual(f.eval([3, 5]), 2*np.sqrt(25+81))
+        # self.assertEqual(f.eval(np.fft.ifft([2, 4])*np.sqrt(2)), 0)
+        # self.assertEqual(f.eval([3, 5]), 2*np.sqrt(25+81))
         nptest.assert_allclose(f.grad([2, 4]), 0)
-#        nptest.assert_allclose(f.grad([3, 5]), [4*np.sqrt(5), 4*3])
+        # nptest.assert_allclose(f.grad([3, 5]), [4*np.sqrt(5), 4*3])
         nptest.assert_allclose(f.prox([2, 4], 1), [2, 4])
         nptest.assert_allclose(f.prox([3, 5], 1), [2.2, 4.2])
         nptest.assert_allclose(f.prox([2.2, 4.2], 1), [2.04, 4.04])
@@ -417,6 +419,46 @@ def test_proj_b2(self):
         f.method = 'NOT_A_VALID_METHOD'
         self.assertRaises(ValueError, f.prox, x, 0)
 
+    def test_proj_lineq(self):
+        """
+        Test the projection on Ax = y
+
+        """
+        x = np.zeros([10])
+        A = np.ones([1, 10])
+        y = np.array([10])
+        f = functions.proj_lineq(A=A, y=y)
+        sol = f.prox(x, 0)
+        np.testing.assert_allclose(sol, np.ones([10]))
+        np.testing.assert_allclose(A.dot(sol), y)
+
+        f = functions.proj_lineq(A=A)
+        sol = f.prox(x, 0)
+        np.testing.assert_allclose(sol, np.zeros([10]))
+
+        for i in range(1, 15):
+            x = np.random.randn(10)
+            y = np.random.randn(i)
+            A = np.random.randn(i, 10)
+            pinvA = np.linalg.pinv(A)
+            f1 = functions.proj_lineq(A=A, y=y)
+            f2 = functions.proj_lineq(A=lambda x: A.dot(x), pinvA=pinvA, y=y)
+            f3 = functions.proj_lineq(A=A, pinvA=lambda x: pinvA.dot(x), y=y)
+            f4 = functions.proj_lineq(A=A, pinvA=pinvA, y=y)
+            sol1 = f1.prox(x, 0)
+            sol2 = f2.prox(x, 0)
+            sol3 = f3.prox(x, 0)
+            sol4 = f4.prox(x, 0)
+            np.testing.assert_allclose(sol1, sol2)
+            np.testing.assert_allclose(sol1, sol3)
+            np.testing.assert_allclose(sol1, sol4)
+            if i <= x.size:
+                np.testing.assert_allclose(A.dot(sol1), y)
+            if i >= x.size:
+                np.testing.assert_allclose(sol1, pinvA.dot(y))
+
+        self.assertRaises(ValueError, functions.proj_lineq, A=lambda x: x)
+
     def test_proj_positive(self):
         """
         Test the projection on the positive octant.
@@ -430,6 +472,30 @@ def test_proj_positive(self):
         nptest.assert_equal(res[x > 0], x[x > 0])  # Positives are unchanged.
         self.assertEqual(fpos.eval(x), 0)
 
+    def test_proj_spsd(self):
+        """
+        Test the projection on symmetric positive semi-definite matrices.
+
+        """
+        f_spds = functions.proj_spsd()
+        A = np.random.randn(10, 10)
+        A = A + A.T
+        eig1 = np.sort(np.real(np.linalg.eig(A)[0]))
+        res = f_spds.prox(A, T=1)
+        eig2 = np.sort(np.real(np.linalg.eig(res)[0]))
+        # All eigenvalues are positive.
+        assert ((eig2 > -1e-13).all())
+
+        # Positive value are unchanged.
+        np.testing.assert_allclose(eig2[eig1 > 0], eig1[eig1 > 0])
+
+        # The symmetrization works.
+        A = np.random.rand(10, 10) + 10 * np.eye(10)
+        res = f_spds.prox(A, T=1)
+        np.testing.assert_allclose(res, (A + A.T) / 2)
+
+        self.assertEqual(f_spds.eval(A), 0)
+
     def test_structured_sparsity(self):
         """
         Test the structured sparsity function.
@@ -503,8 +569,9 @@ def test_independent_problems(self):
             if name == 'norm_tv':
                 # Each column is one-dimensional.
                 f = func(dim=1, maxit=20, tol=0)
-            elif name == 'norm_nuclear':
-                # TODO: make this test two dimensional for the norm nuclear?
+            elif name in ['norm_nuclear', 'proj_spsd']:
+                # TODO: make this test two dimensional for the norm nuclear
+                # and the spsd projection?
                 continue
             else:
                 f = func()