Merge 0281a68 into b0eb893

pyunlocbox/functions.py

@@ -36,6 +36,7 @@
 
     proj_positive
     proj_b2
+    proj_spsd
 
 **Miscellaneous**
 
@@ -803,24 +804,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
@@ -867,6 +857,57 @@ def _prox(self, x, T):
         return np.clip(x, 0, np.inf)
 
 
+class proj_spsd(proj):
+    r"""
+    Projection on the Symetric positive semi definite set of matrices
+
+    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\}`.
+
+    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 # Symetrize 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 symetric
+        sol = (x + np.conj(x.T)) / 2
+
+        # 2. make 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).
@@ -880,6 +921,14 @@ class proj_b2(proj):
     :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
@@ -904,10 +953,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):
 

pyunlocbox/tests/test_functions.py

@@ -95,16 +95,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])
@@ -430,6 +430,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 the positive octant.
+
+        """
+        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-15).all())
+
+        # Positive value are unchanged
+        np.testing.assert_allclose(eig2[eig1 > 0], eig1[eig1 > 0])
+
+        # The symetrization 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 +527,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()