Merge pull request #76 from marcusvaltonen/user/marcus/separable-singular-value-penalties

Implement generic singular value penalty

Author: mrava87

docs/source/api/index.rst

@@ -80,6 +80,7 @@ Operators
     ETP
     Geman
     QuadraticEnvelopeCard
+    SingularValuePenalty
 
 
 Other operators

pyproximal/proximal/SingularValuePenalty.py

@@ -0,0 +1,59 @@
+import numpy as np
+
+from pyproximal.ProxOperator import _check_tau
+from pyproximal import ProxOperator
+
+
+class SingularValuePenalty(ProxOperator):
+    r"""Proximal operator of a penalty acting on the singular values.
+
+    Generic regularizer :math:`\mathcal{R}_f` acting on the singular values of a matrix,
+
+    .. math::
+
+        \mathcal{R}_f(\mathbf{X}) = f(\boldsymbol\lambda)
+
+    where :math:`\mathbf{X}` is a matrix of size :math:`M \times N` and
+    :math:`\boldsymbol\lambda` is the corresponding singular value vector.
+
+    Parameters
+    ----------
+    dim : :obj:`tuple`
+        Size of matrix :math:`\mathbf{X}`.
+    penalty : :class:`pyproximal.ProxOperator`
+        Function acting on the singular values.
+
+    Notes
+    -----
+    The pyproximal implementation allows ``penalty`` to be any
+    :class:`pyproximal.ProxOperator` acting on the singular values; however, not all
+    penalties will result in a mathematically accurate proximal operator defined this
+    way. Given a penalty :math:`f`, the proximal operator is assumed to be
+
+    .. math::
+
+        \prox_{\tau \mathcal{R}_f}(\mathbf{X}) =
+        \mathbf{U} \diag\left( \prox_{\tau f}(\boldsymbol\lambda)\right) \mathbf{V}^H
+
+    where :math:`\mathbf{X} = \mathbf{U}\diag(\boldsymbol\lambda)\mathbf{V}^H`, is an
+    SVD of :math:`\mathbf{X}`. It is the user's responsibility to check that this is
+    true for their particular choice of ``penalty``.
+    """
+
+    def __init__(self, dim, penalty):
+        super().__init__(None, False)
+        self.dim = dim
+        self.penalty = penalty
+
+    def __call__(self, x):
+        X = x.reshape(self.dim)
+        eigs = np.linalg.eigvalsh(X.T @ X)
+        eigs[eigs < 0] = 0  # ensure all eigenvalues at positive
+        return np.sum(self.penalty(np.sqrt(eigs)))
+
+    @_check_tau
+    def prox(self, x, tau):
+        X = x.reshape(self.dim)
+        U, S, Vh = np.linalg.svd(X, full_matrices=False)
+        X = np.dot(U * self.penalty.prox(S, tau), Vh)
+        return X.ravel()

pyproximal/proximal/__init__.py

@@ -29,6 +29,7 @@
     ETP                             Exponential-type penalty
     Geman                           Geman penalty
     QuadraticEnvelopeCard           The quadratic envelope of the cardinality function
+    SingularValuePenalty            Generic singular value penalty
 
 """
 
@@ -53,9 +54,10 @@
 from .ETP import *
 from .Geman import *
 from .QuadraticEnvelope import *
+from .SingularValuePenalty import *
 
 __all__ = ['Box', 'Simplex', 'Intersection', 'AffineSet', 'Quadratic',
            'Euclidean', 'EuclideanBall', 'L0Ball', 'L1', 'L1Ball', 'L2',
            'L2Convolve', 'L21', 'L21_plus_L1', 'Huber', 'Nuclear',
            'NuclearBall', 'Orthogonal', 'VStack', 'Nonlinear', 'SCAD',
-           'Log', 'ETP', 'Geman', 'QuadraticEnvelopeCard']
+           'Log', 'ETP', 'Geman', 'QuadraticEnvelopeCard', 'SingularValuePenalty']

pytests/test_proximal.py

@@ -3,20 +3,21 @@
 import numpy as np
 from numpy.testing import assert_array_equal, assert_array_almost_equal
 from pylops import MatrixMult, Identity
+
+import pyproximal
 from pyproximal.utils import moreau
-from pyproximal.proximal import Quadratic, Nonlinear, \
-    L1, L2, Orthogonal, VStack
+from pyproximal.proximal import L1, L2, Nonlinear, Orthogonal, Quadratic, \
+    SingularValuePenalty, VStack
 
 par1 = {'nx': 10, 'sigma': 1., 'dtype': 'float32'}  # even float32
 par2 = {'nx': 11, 'sigma': 2., 'dtype': 'float64'}  # odd float64
 
-np.random.seed(10)
-
 
 @pytest.mark.parametrize("par", [(par1), (par2)])
 def test_Quadratic(par):
     """Quadratic functional and proximal/dual proximal
     """
+    np.random.seed(10)
     A = np.random.normal(0, 1, (par['nx'], par['nx']))
     A = A.T @ A
     quad = Quadratic(Op=MatrixMult(A), b=np.ones(par['nx']), niter=500)
@@ -31,6 +32,7 @@ def test_Quadratic(par):
 def test_DotProduct(par):
     """Dot product functional and proximal/dual proximal
     """
+    np.random.seed(10)
     quad = Quadratic(b=np.ones(par['nx']))
 
     # prox / dualprox
@@ -43,6 +45,7 @@ def test_DotProduct(par):
 def test_Constant(par):
     """Constant functional and proximal/dual proximal
     """
+    np.random.seed(10)
     quad = Quadratic(c=5.)
 
     # prox / dualprox
@@ -55,6 +58,7 @@ def test_Constant(par):
 def test_SemiOrthogonal(par):
     """L1 functional with Semi-Orthogonal operator and proximal/dual proximal
     """
+    np.random.seed(10)
     l1 = L1()
     orth = Orthogonal(l1, 2*Identity(par['nx']), b=np.arange(par['nx']),
                       partial=True, alpha=4.)
@@ -69,6 +73,7 @@ def test_SemiOrthogonal(par):
 def test_Orthogonal(par):
     """L1 functional with Orthogonal operator and proximal/dual proximal
     """
+    np.random.seed(10)
     l1 = L1()
     orth = Orthogonal(l1, Identity(par['nx']), b=np.arange(par['nx']))
 
@@ -82,6 +87,7 @@ def test_Orthogonal(par):
 def test_VStack(par):
     """L2 functional with VStack operator of multiple L1s
     """
+    np.random.seed(10)
     nxs = [par['nx'] // 4] * 4
     nxs[-1] = par['nx'] - np.sum(nxs[:-1])
     l2 = L2()
@@ -106,6 +112,7 @@ def test_Nonlinear():
     """Nonlinear proximal operator. Since this is a template class simply check
     that errors are raised when not used properly
     """
+    np.random.seed(10)
     Nop = Nonlinear(np.ones(10))
     with pytest.raises(NotImplementedError):
         Nop.fun(np.ones(10))
@@ -115,4 +122,20 @@ def test_Nonlinear():
         Nop.optimize()
 
 
+@pytest.mark.parametrize("par", [(par1), (par2)])
+def test_SingularValuePenalty(par):
+    """Test SingularValuePenalty
+    """
+    np.random.seed(10)
+    f_mu = pyproximal.QuadraticEnvelopeCard(mu=par['sigma'])
+    penalty = SingularValuePenalty((par['nx'], 2 * par['nx']), f_mu)
+
+    # norm, cross-check with svd (use tolerance as two methods don't provide
+    # the exact same eigenvalues)
+    X = np.random.uniform(0., 0.1, (par['nx'], 2 * par['nx'])).astype(par['dtype'])
+    _, S, _ = np.linalg.svd(X)
+    assert (penalty(X.ravel()) - f_mu(S)) < 1e-3
 
+    # prox / dualprox
+    tau = 0.75
+    assert moreau(penalty, X.ravel(), tau)