# Source code for pyproximal.proximal.Orthogonal

import numpy as np

from scipy.sparse.linalg import lsqr
from pylops import MatrixMult, Identity
from pyproximal.ProxOperator import _check_tau
from pyproximal import ProxOperator

[docs]class Orthogonal(ProxOperator): r"""Proximal operator of any function of the product between an orthogonal matrix and a vector (plus summation of a vector). Proximal operator of any quadratic function :math:g(\mathbf{x})=f(\mathbf{Qx} + \mathbf{b}) where :math:\mathbf{Q} is an orthogonal operator and :math:\mathbf{b} is a vector in the data space of :math:\mathbf{Q}. Parameters ---------- f : :obj:pyproximal.ProxOperator Proximal operator Q : :obj:pylops.LinearOperator Orthogonal operator partial : :obj:bool, optional Partial (True) of full (False) orthogonality b : :obj:np.ndarray, optional Vector alpha : :obj:float, optional Positive coefficient for partial orthogonality. It will be ignored if partial=False Notes ----- The Proximal operator of any function of the form :math:g(\mathbf{x}) = f(\mathbf{Qx} + \mathbf{b}) with the operator :math:\mathbf{Q} satisfying the following condition :math:\mathbf{Q}\mathbf{Q}^T = \alpha \mathbf{I} (partial orthogonality) is [1]_: .. math:: \prox_{\tau g}(x) = \frac{1}{\alpha} ((\alpha \mathbf{I} - \mathbf{Q}^H \mathbf{Q}) \mathbf{x} + \mathbf{Q}^H (\prox_{\alpha \tau f}(\mathbf{Qx} + \mathbf{b}) - \mathbf{b})) A special case arises when :math:\mathbf{Q}\mathbf{Q}^T = \mathbf{Q}^T\mathbf{Q} = \mathbf{I} (full orthogonality), and the proximal operator reduces to: .. math:: \prox_{\tau g}(x) = \mathbf{Q}^H (\prox_{\tau f}(\mathbf{Qx} + \mathbf{b}) - \mathbf{b})) .. [1] Daniel O'Connor, D., and Vandenberghe, L., "Primal-Dual Decomposition by Operator Splitting and Applications to Image Deblurring", SIAM J. Imaging Sciences, vol. 7, pp. 1724–1754. 2014. """ def __init__(self, f, Q, partial=False, b=None, alpha=1.): super().__init__(None, False) self.f = f self.Q = Q self.partial = partial self.alpha = alpha self.b = b if b is not None else 0 def __call__(self, x): y = self.Q.matvec(x) y += self.b f = self.f(y) return f @_check_tau def prox(self, x, tau): y = self.Q.matvec(x) if self.partial: z = (1. / self.alpha) * \ (self.alpha * x - self.Q.rmatvec(y) + self.Q.rmatvec(self.f.prox(y + self.b, self.alpha * tau) - self.b)) else: y = y + self.b z = self.f.prox(y, tau) z = z - self.b z = self.Q.rmatvec(z) return z