Tutorial: Hankel matrix recovery

Author: marcusvaltonen

docs/source/api/index.rst

@@ -20,6 +20,7 @@ Orthogonal projections
     AffineSetProj
     BoxProj
     EuclideanBallProj
+    HankelProj
     HyperPlaneBoxProj
     IntersectionProj
     L0BallProj
@@ -62,6 +63,7 @@ Convex
     Box
     Euclidean
     EuclideanBall
+    Hankel
     Huber
     Intersection
     L0
@@ -88,6 +90,7 @@ Non-Convex
     Geman
     Log
     QuadraticEnvelopeCard
+    QuadraticEnvelopeCardIndicator
     SCAD
 
 
@@ -100,7 +103,7 @@ Matrix-only
     Nuclear
     NuclearBall
     SingularValuePenalty
-
+    QuadraticEnvelopeRankL2
 
 Other
 ^^^^^

docs/source/conf.py

@@ -140,6 +140,7 @@
             "sgn": r"\operatorname{sgn}",
             "argmin": r"\operatorname*{argmin}",
             "diag": r"\operatorname{diag}",
+            "rank": r"\operatorname{rank}",
         }
     }
 }

pyproximal/projection/Hankel.py

@@ -0,0 +1,26 @@
+import numpy as np
+from scipy.linalg import hankel
+
+
+class HankelProj:
+    r"""Hankel matrix projection.
+
+    Solves the least squares problem
+
+        .. math::
+          \min_{X\in\mathcal{H}} \|X-X_0\|_F^2
+
+    where :math:`\mathcal{H}` is the set of Hankel matrices.
+
+    Notes
+    -----
+    The solution to the above-mentioned least squares problem is given by a Hankel matrix,
+    where the (constant) anti-diagonals are the average value along the corresponding
+    anti-diagonals of the original matrix :math:`X_0`.
+    """
+
+    def __call__(self, X):
+        m, n = X.shape
+        ind = hankel(np.arange(m, dtype=np.int32), m - 1 + np.arange(n, dtype=np.int32))
+        mean_values = np.bincount(ind.ravel(), weights=X.ravel()) / np.bincount(ind.ravel())
+        return mean_values[ind]

pyproximal/projection/__init__.py

@@ -13,6 +13,7 @@
     NuclearBallProj	            Projection onto a Nuclear Ball
     IntersectionProj	        Projection onto an Intersection of sets
     AffineSetProj	            Projection onto an Affine set
+    HankelProj                  Projection onto the set of Hankel matrices
 
 """
 
@@ -24,8 +25,9 @@
 from .Nuclear import *
 from .Intersection import *
 from .AffineSet import *
+from .Hankel import *
 
 
 __all__ = ['BoxProj', 'HyperPlaneBoxProj', 'SimplexProj', 'L0BallProj',
            'L1BallProj', 'EuclideanBallProj', 'NuclearBallProj',
-           'IntersectionProj', 'AffineSetProj']
+           'IntersectionProj', 'AffineSetProj', 'HankelProj']

pyproximal/proximal/Hankel.py

@@ -0,0 +1,36 @@
+import numpy as np
+from pyproximal.ProxOperator import _check_tau
+from pyproximal import ProxOperator
+from pyproximal.projection import HankelProj
+
+
+class Hankel(ProxOperator):
+    r"""Hankel proximal operator.
+
+    Proximal operator of the Hankel matrix indicator function.
+
+    Parameters
+    ----------
+    dim : :obj:`tuple`
+        Dimension of the Hankel matrix.
+
+    Notes
+    -----
+    As the Hankel Operator is an indicator function, the proximal operator corresponds to
+    its orthogonal projection (see :class:`pyproximal.projection.HankelProj` for
+    details).
+
+    """
+    def __init__(self, dim):
+        super().__init__(None, False)
+        self.dim = dim
+        self.hankel_proj = HankelProj()
+
+    def __call__(self, x):
+        X = x.reshape(self.dim)
+        return np.allclose(X, self.hankel_proj(X))
+
+    @_check_tau
+    def prox(self, x, tau):
+        X = x.reshape(self.dim)
+        return self.hankel_proj(X).ravel()

pyproximal/proximal/QuadraticEnvelope.py

@@ -1,5 +1,6 @@
 import numpy as np
 
+from pylops.optimization.sparsity import _hardthreshold
 from pyproximal.ProxOperator import _check_tau
 from pyproximal import ProxOperator
 
@@ -21,6 +22,10 @@ class QuadraticEnvelopeCard(ProxOperator):
     mu : :obj:`float`
         Threshold parameter.
 
+    See Also
+    --------
+    QuadraticEnvelopeCardIndicator: Quadratic envelope of the indicator function of :math:`\ell_0`-penalty
+
     Notes
     -----
     The terminology *quadratic envelope* was coined in [1]_, however, the rationale has
@@ -57,6 +62,8 @@ class QuadraticEnvelopeCard(ProxOperator):
     that this proximal operator is identical to the Minimax Concave Penalty (MCP)
     proposed in [3]_.
 
+    References
+    ----------
     .. [1] Carlsson, M. "On Convex Envelopes and Regularization of Non-convex
         Functionals Without Moving Global Minima", In Journal of Optimization Theory
         and Applications, 183:66–84, 2019.
@@ -86,3 +93,226 @@ def prox(self, x, tau):
         else:
             r[idx] = np.maximum(0, (r[idx] - tau * np.sqrt(2 * self.mu)) / (1 - tau))
         return r * np.sign(x)
+
+
+class QuadraticEnvelopeCardIndicator(ProxOperator):
+    r"""Quadratic envelope of the indicator function of the :math:`\ell_0`-penalty.
+
+    The :math:`\ell_0`-penalty is also known as the *cardinality function*, and the
+    indicator function :math:`\mathcal{I}_{r_0}` is defined as
+
+    .. math::
+
+        \mathcal{I}_{r_0}(\mathbf{x}) =
+        \begin{cases}
+        0, & \mathbf{x}\leq r_0 \\
+        \infty, & \text{otherwise}
+        \end{cases}
+
+    Let :math:`\tilde{\mathbf{x}}` denote the vector :math:`\mathbf{x}` resorted such that the
+    sequence :math:`(\tilde{x}_i)` is non-increasing. The quadratic envelope
+    :math:`\mathcal{Q}(\mathcal{I}_{r_0})` can then be written as
+
+    .. math::
+
+        \mathcal{Q}(\mathcal{I}_{r_0})(x) =
+        \frac{1}{2k^*}\left(\sum_{i>r_0-k^*}|\tilde{x}_i|\right)^2
+        - \frac{1}{2}\left(\sum_{i>r_0-k^*}|\tilde{x}_i|^2
+
+    where :math:`r_0 \geq 0` and :math:`k^* \leq r_0`, see [3]_ for details. There are
+    other, equivalent ways, of expressing this penalty, see e.g. [1]_ and [2]_.
+
+    Parameters
+    ----------
+    r0 : :obj:`int`
+        Threshold parameter.
+
+    See Also
+    --------
+    QuadraticEnvelopeCard: Quadratic envelope of the :math:`\ell_0`-penalty
+
+    Notes
+    -----
+    The terminology *quadratic envelope* was coined in [1]_, however, the rationale has
+    been used earlier, e.g. in [2]_. In a general setting, the quadratic envelope
+    :math:`\mathcal{Q}(f)(x)` is defined such that
+
+    .. math::
+
+        \left(f(x) + \frac{1}{2}\|x-y\|_2^2\right)^{**} = \mathcal{Q}(f)(x) + \frac{1}{2}\|x-y\|_2^2
+
+    where :math:`g^{**}` denotes the bi-conjugate of :math:`g`, which is the l.s.c.
+    convex envelope of :math:`g`.
+
+    There is no closed-form expression for :math:`\mathcal{Q}(f)(x)` given an arbitrary
+    function :math:`f`. However, for certain special cases, such as in the case of the
+    indicator function of the cardinality function, such expressions do exist.
+
+    The proximal operator does not have a closed-form, and we refer to [1]_ for more details.
+    Note that this is a non-separable penalty.
+
+    References
+    ----------
+    .. [1] Carlsson, M. "On Convex Envelopes and Regularization of Non-convex
+        Functionals Without Moving Global Minima", In Journal of Optimization Theory
+        and Applications, 183:66–84, 2019.
+    .. [2] Larsson, V. and Olsson, C. "Convex Low Rank Approximation", In International
+        Journal of Computer Vision (IJCV), 120:194–214, 2016.
+    .. [3] Andersson et al. "Convex envelopes for fixed rank approximation", In
+        Optimization Letters, 11:1783–1795, 2017.
+
+    """
+
+    def __init__(self, r0):
+        super().__init__(None, False)
+        self.r0 = r0
+
+    def __call__(self, x):
+        if x.size <= self.r0 or np.count_nonzero(x) <= self.r0:
+            return 0
+        xs = np.sort(np.abs(x))[::-1]
+        sums = np.cumsum(xs[::-1])
+        sums = sums[-self.r0:] / np.arange(1, self.r0 + 1)
+        tmp = np.diff(sums) > 0
+        k_star = np.argmax(tmp)
+        if k_star == 0 and not tmp[k_star]:
+            k_star = self.r0 - 1
+        return 0.5 * ((k_star + 1) * sums[k_star] ** 2 - np.sum(xs[self.r0-k_star-1:] ** 2))
+
+    @_check_tau
+    def prox(self, y, tau):
+        rho = 1 / tau
+        if rho <= 1:
+            return _hardthreshold(y, tau)
+        if y.size <= self.r0:
+            return y
+
+        r = np.abs(y)
+        theta = np.sign(y)
+        id = np.argsort(-r, kind='quicksort')
+        rsort = r[id]
+        idinv = np.zeros_like(id)
+        idinv[id] = np.arange(r.size)
+        rnew = np.concatenate((rsort[:self.r0], rho * rsort[self.r0:]))
+
+        if rho * rsort[self.r0] < rsort[self.r0 - 1]:
+            x = rnew
+            x = x[idinv]
+            x = x * theta
+        else:
+            j = np.min(np.where(rnew <= rnew[self.r0])[0])
+            l = np.max(np.where(rnew >= rnew[self.r0 - 1])[0])
+            z = np.sort(rnew[j:l + 1])[::-1]
+            z1 = z[0]
+            for z2 in z[1:]:
+                s = (z1 + z2) / 2
+                temp = np.where(rnew <= s)[0]
+                j1 = np.min(temp)
+                temp = np.where(rnew >= s)[0]
+                l1 = np.max(temp)
+                sI = (rho * sum(rsort[j1:l1 + 1])) / ((self.r0 - j1) * rho + (l1 + 1 - self.r0) * 1)
+                if z2 <= sI <= z1:
+                    x = np.concatenate((np.maximum(rnew[:self.r0], sI), np.minimum(rnew[self.r0:], sI)))
+                    x = x[idinv]
+                    x = x * theta
+                    break
+                z1 = z2
+
+        return (rho * y - x) / (rho - 1)
+
+
+class QuadraticEnvelopeRankL2(ProxOperator):
+    r"""Quadratic envelope of the rank function with an L2 misfit term.
+
+    The penalty :math:`p` is given by
+
+     .. math::
+
+        p(X) = \mathcal{R}_{r_0}(X) + \frac{1}{2}\|X - M\|_F^2
+
+    where :math:`\mathcal{R}_{r_0}` is the quadratic envelope of the hard-rank function.
+
+    Parameters
+    ----------
+    dim : :obj:`tuple`
+        Size of input matrix :math:`X`.
+    r0 : :obj:`int`
+        Threshold parameter, encouraging matrices with rank lower than or equal to r0.
+    M : :obj:`numpy.ndarray`
+        L2 misfit term (must be the same size as the input matrix).
+
+    See Also
+    --------
+    SingularValuePenalty: Proximal operator of a penalty acting on the singular values
+    QuadraticEnvelopeCardIndicator: Quadratic envelope of the indicator function of :math:`\ell_0`-penalty
+
+    Notes
+    -----
+    The proximal operator solves the minimization problem
+
+        .. math::
+            \argmin_Z \mathcal{R}_{r_0}(Z) + \frac{1}{2}\|Z - M\|_F^2 + \frac{1}{2\tau}\| Z - X \|_F^2
+
+    which is a convex-concave min-max problem, see [1]_ for details.
+
+    References
+    ----------
+    .. [1] Larsson, V. and Olsson, C. "Convex Low Rank Approximation", In International
+        Journal of Computer Vision (IJCV), 120:194–214, 2016.
+
+    """
+
+    def __init__(self, dim, r0, M):
+        super().__init__(None, False)
+        self.dim = dim
+        self.r0 = r0
+        self.M = M.copy()
+        self.penalty = QuadraticEnvelopeCardIndicator(r0)
+
+    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))) + 0.5 * np.linalg.norm(X - self.M, 'fro')
+
+    @_check_tau
+    def prox(self, x, tau):
+        rho = 1 / tau
+        P = x.reshape(self.dim)
+
+        Y = (self.M + rho * P) / (1 + rho)
+        U, yk, Vh = np.linalg.svd(Y, full_matrices=False)
+        n = yk.size
+        r = np.concatenate((yk[:self.r0], (1 + rho) * yk[self.r0:]))
+        ind = np.argsort(r, kind='quicksort')
+        p = r[ind]
+
+        a = (n - self.r0) / rho
+        b = (rho + 1) / rho * np.sum(yk[self.r0:])
+
+        # Base case
+        zk = yk.copy()
+        zk[self.r0:] = (1 + rho) * yk[self.r0:]
+
+        for k, ii in enumerate(ind):
+
+            if ii < self.r0:
+                a = a + (rho + 1) / rho
+                b = b + (rho + 1) / rho * yk[ii]
+            else:
+                a = a - 1 / rho
+                b = b - (rho + 1) / rho * yk[ii]
+
+            if a == 0:
+                continue
+
+            s = b / a
+
+            if p[k] <= s <= p[k + 1]:
+                zk = np.maximum(s, yk)
+                zk[self.r0:] = np.minimum(s, (1 + rho) * yk[self.r0:])
+                break
+
+        Z = np.dot(U * zk, Vh)
+        X = P + (self.M - Z) / rho
+        return X.ravel()

pyproximal/proximal/SingularValuePenalty.py

@@ -20,7 +20,7 @@ class SingularValuePenalty(ProxOperator):
     ----------
     dim : :obj:`tuple`
         Size of matrix :math:`\mathbf{X}`.
-    penalty : :class:`pyproximal.ProxOperator`
+    penalty : :obj:`pyproximal.ProxOperator`
         Function acting on the singular values.
 
     Notes