From b06e9b92e9469800cad87b03c6bfac19548e1bb1 Mon Sep 17 00:00:00 2001
From: Megumi1773 <lixiaojie0514@foxmail.com>
Date: Thu, 8 Oct 2026 22:08:56 +0800
Subject: [PATCH] feat: add Lipschitz method to LinearOperator

Add LinearOperator.Lipschitz to compute the Lipschitz constant of an
operator, defined as its largest singular value
(sqrt(max eigenvalue of Op.H * Op)). A square=True flag returns the
squared value directly, which is what gradient-based solvers need for
their step size.

This replaces the recurring Op1 = Op.H * Op; L = abs(Op1.eigs(...))
pattern seen across examples and solvers. For the numpy backend the
largest eigenvalue is estimated with LinearOperator.eigs (symmetric);
for cupy/jax it defers to optimization.eigs.power_iteration.

Closes #812
---
 pylops/linearoperator.py       | 67 ++++++++++++++++++++++++++++++++++
 pytests/test_linearoperator.py | 28 ++++++++++++++
 2 files changed, 95 insertions(+)

diff --git a/pylops/linearoperator.py b/pylops/linearoperator.py
index 7a47ac0d..1a1b4d37 100644
--- a/pylops/linearoperator.py
+++ b/pylops/linearoperator.py
@@ -1211,6 +1211,73 @@ def cond(
 
         return cond
 
+    def Lipschitz(
+        self,
+        square: bool = False,
+        backend: str = "numpy",
+        **kwargs_eig: int | float | str,
+    ) -> float:
+        r"""Lipschitz constant of linear operator.
+
+        Return the Lipschitz constant of the linear operator, defined as the
+        largest singular value of the operator
+
+        .. math::
+            L = \|\mathbf{A}\|_2 = \sqrt{\lambda_{max}(\mathbf{A}^H\mathbf{A})}
+
+        i.e. the square root of the largest eigenvalue of the normal operator
+        :math:`\mathbf{A}^H\mathbf{A}`. This is the constant that appears in
+        the convergence condition of gradient-based solvers, so it is
+        frequently needed to select a step size (for example
+        :math:`\alpha = 1/L^2` in :func:`pylops.optimization.sparsity.ista`).
+
+        Parameters
+        ----------
+        square : :obj:`bool`, optional
+            Return the square of the Lipschitz constant
+            :math:`\lambda_{max}(\mathbf{A}^H\mathbf{A})` (``True``) instead of
+            the Lipschitz constant itself (``False``). The squared value can be
+            used directly where solvers need :math:`\lambda_{max}`.
+        backend : :obj:`str`, optional
+            Backend used to compute the largest eigenvalue (``numpy`` or
+            ``cupy`` or ``jax``). Note that this must be consistent with how the
+            operator has been created.
+        **kwargs_eig
+            Arbitrary keyword arguments for
+            :func:`pylops.LinearOperator.eigs` (``numpy`` backend) or
+            :func:`pylops.optimization.eigs.power_iteration` (``cupy``/``jax``
+            backends)
+
+        Returns
+        -------
+        lipschitz : :obj:`float`
+            Lipschitz constant (or its square if ``square=True``).
+
+        Notes
+        -----
+        For the ``numpy`` backend the largest eigenvalue of the normal operator
+        is estimated with :func:`pylops.LinearOperator.eigs` (which is
+        symmetric, as :math:`\mathbf{A}^H\mathbf{A}` is Hermitian). For the
+        ``cupy`` and ``jax`` backends it is estimated with
+        :func:`pylops.optimization.eigs.power_iteration`.
+
+        """
+        Op1 = self.H * self
+        if backend == "numpy":
+            maxeig = Op1.eigs(neigs=1, symmetric=True, **kwargs_eig)[0]
+        else:
+            from pylops.optimization.eigs import power_iteration
+
+            maxeig = power_iteration(
+                Op1,
+                dtype=Op1.dtype,
+                backend=backend,
+                **kwargs_eig,
+            )[0]
+
+        maxeig = np.abs(maxeig)
+        return float(maxeig) if square else float(np.sqrt(maxeig))
+
     def conj(self) -> LinearOperator:
         """Complex conjugate operator
 
diff --git a/pytests/test_linearoperator.py b/pytests/test_linearoperator.py
index 29eba9c8..5a69f57e 100644
--- a/pytests/test_linearoperator.py
+++ b/pytests/test_linearoperator.py
@@ -178,6 +178,34 @@ def test_eigs(par):
         assert_array_almost_equal(np.real(cond), np.real(cond1), decimal=3)
 
 
+@pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j)])
+def test_lipschitz(par):
+    """Lipschitz constant (and its square) of a linear operator"""
+    diag = np.arange(par["nx"], 0, -1) + par["imag"] * np.arange(par["nx"], 0, -1)
+    Op = Diagonal(diag, dtype=par["dtype"])
+    if par["ny"] > par["nx"]:
+        Op = VStack([Op, Zero(par["ny"] - par["nx"], par["nx"])])
+
+    # largest singular value of a diagonal (padded) matrix is the largest
+    # entry of the diagonal
+    lipschitz = Op.Lipschitz()
+    lipschitz_np = np.abs(diag).max()
+    assert_array_almost_equal(lipschitz, lipschitz_np, decimal=3)
+
+    # square=True returns the square of the Lipschitz constant
+    lipschitz_sq = Op.Lipschitz(square=True)
+    assert_array_almost_equal(lipschitz_sq, lipschitz**2, decimal=3)
+
+    # same result against a non-diagonal operator via the singular values
+    M = np.arange(par["ny"] * par["nx"]).reshape(par["ny"], par["nx"]) + par[
+        "imag"
+    ] * np.ones((par["ny"], par["nx"]))
+    Op = MatrixMult(M, dtype=par["dtype"])
+    lipschitz = Op.Lipschitz()
+    lipschitz_np = np.linalg.svd(M, compute_uv=False)[0]
+    assert_array_almost_equal(lipschitz, lipschitz_np, decimal=3)
+
+
 @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j)])
 def test_conj(par):
     """Complex conjugation operator"""