CGLS Solver

This example shows how to use the pylops.optimization.leastsquares.cgls solver to minimize the following cost function:

\[J = || \mathbf{y} - \mathbf{Ax} ||_2^2 + \epsilon || \mathbf{x} ||_2^2\]

import warnings
import numpy as np

import matplotlib.pyplot as plt

import pylops

plt.close('all')
warnings.filterwarnings('ignore')

Let’s define a matrix \(\mathbf{A}\) or size (N and M) and fill the matrix with random numbers

N, M = 20, 10
A = np.random.normal(0, 1, (N, M))
Aop = pylops.MatrixMult(A, dtype='float64')

x = np.ones(M)

We can now use the cgls solver to invert this matrix

y = Aop * x
xest, istop, nit, r1norm, r2norm, cost = \
    pylops.optimization.solver.cgls(Aop, y, x0=np.zeros_like(x),
                                    niter=10, tol=1e-10, show=True)

print('x= %s' % x)
print('cgls solution xest= %s' % xest)

plt.figure(figsize=(12, 3))
plt.plot(cost, 'k', lw=2)
plt.title('Cost function')

Out:

CGLS
-----------------------------------------------------------
The Operator Op has 20 rows and 10 cols
damp = 0.000000e+00     tol = 1.000000e-10      niter = 10
-----------------------------------------------------------------
    Itn           x[0]              r2norm
     1         9.1362e-01         1.7415e+02
     2         1.1328e+00         6.5026e+01
     3         1.1030e+00         1.2698e+01
     4         1.0366e+00         1.4992e+00
     5         1.0086e+00         2.1434e-01
     6         1.0069e+00         9.3148e-02
     7         9.9981e-01         1.8127e-02
     8         9.9936e-01         4.3161e-03
     9         1.0006e+00         5.7191e-05
    10         1.0000e+00         4.8266e-28

Iterations = 10        Total time (s) = 0.00
-----------------------------------------------------------------

x= [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
cgls solution xest= [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]

Text(0.5, 1.0, 'Cost function')

Note that while we used a dense matrix here, any other linear operator can be fed to cgls as is the case for any other PyLops solver.