Motivation

A code pattern that arises in many of our examples/tutorials/notebooks (especially with pyproximal solvers) is along the lines of:

Op1 = Op.H * Op
L = np.abs(Op1.eigs(1, niter=5, tol=5e-2))[0]
tau = ... # some formula using L
...

solver(tau=tau)

This calls for a new method in pylops.LinearOperator called Lipschitz which computes the Lipschitz constant of the operator = sqrt(max eigenvalue of Op.H * Op), and possibly a square parameter that returns the square (to be used directly in solvers). This method should call the eigs method or for cupy/jax arrays defer to pylops.optimization.eigs.power_iteration.

This also enables one to overwrite this method in solvers that have closed-form solution - eg Diagonal,Restriction,Blending, FFT, etc.