Authors
Affiliations
Francesco Picetti
Image and Sound Processing Lab, Politecnico di Milano
Ettore Biondi
SeismoLab, California Institute of Technology

@Author: Ettore Biondi - ebiondi@caltech.edu

## Instantiation of vectors and operator¶

For testing the library we will be using a discretized version of the following operator:

\begin{align} y = \frac{d^2f(x)}{dx^2}, \end{align}

in which we simply compute the second-order derivative of a function $f(x)$.

#### Importing necessary libraries¶

import numpy as np
import occamypy as o

# Plotting
from matplotlib import rcParams
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.pyplot as plt
rcParams.update({
'image.cmap'     : 'gray',
'image.aspect'   : 'auto',
'image.interpolation': None,
'axes.grid'      : True,
'figure.figsize' : (6, 4),
'savefig.dpi'    : 300,
'axes.labelsize' : 14,
'axes.titlesize' : 16,
'font.size'      : 14,
'legend.fontsize': 14,
'xtick.labelsize': 14,
'ytick.labelsize': 14,
'text.usetex'    : True,
'font.family'    : 'serif',
'font.serif'     : 'Latin Modern Roman',
})
N = 200  # Number of points of the function f(x)
dx = 1/N  # Sampling of the function

# The stencil used is simply: (f(ix-1)-2f(ix)+f(ix+1))/(dx*dx)
D2 = np.zeros((N, N), dtype=np.float64)
np.fill_diagonal(D2, -2 / (dx * dx))
np.fill_diagonal(D2[1:], 1 / (dx * dx))
np.fill_diagonal(D2[:,1:], 1 / (dx * dx))

f = o.VectorNumpy((N, 1)).zero()
y = f.clone()

D2 = o.Matrix(o.VectorNumpy(D2), f, y)

Before we set any inversion problem, we study some of the properties of the constructed operator Deriv2Op.

# Verifying operator adjointness through dot-product test
D2.dotTest(verbose=True)

Computing maximum and minimum eigenvalues of the operator using the power iteration method and compare them against the ones computed using numpy

egsOp = D2.powerMethod(verbose=False, eval_min=True, tol=1e-300)
egsNp, _ = np.linalg.eig(D2.matrix.getNdArray())
egsNp = egsNp[egsNp.argsort()[::-1]] #Sorting the eigenvalues
print("\nMaximum eigenvalue: %s (Power method), %s (NUMPY)" % (egsOp[0], egsNp[-1]))
print(  "Minimum eigenvalue: %s (Power method), %s (NUMPY)" % (egsOp[1], egsNp[0]))

We can see that the matrix is negative definite. The small mismatch in the estimated eigenvalues is due to the dependence of the power method on the initial random eigenvector.

## Inversion tests¶

We will now focus our attention on inverting a function knowning its second-order derivative. In this case we will assume that $y$ is constant and equal to $1$. Therefore, we expect to obtain a parabola with positive curvature. Given the chosen boundary conditions we know that the matrix is invertible since all eigenvalues have the same sign and are different then zero. We will solve the following objective functions using linear conjugate-gradient methods:

$\begin{equation*} \phi_1(\mathbf{f}) = \frac{1}{2}\|D_2\mathbf{f}-\mathbf{y}\|_2^2 \end{equation*}$

and

$\begin{equation*} \phi_2(\mathbf{f}) = \frac{1}{2}\mathbf{f}^T D_2 \mathbf{f} - \mathbf{f}^{T} \mathbf{y}, \end{equation*}$

where $D_2$ represents the discretized second-order derivative operator, while $\mathbf{f}$ and $\mathbf{y}$ are the discretized representations of $f$ and $y$, respectively.

y.set(1.) # y = 1
# Note that f = 0
Phi1 = o.LeastSquares(f.clone().zero(), y, D2)
Phi2 = o.LeastSquaresSymmetric(f.clone().zero(), y, D2)

### Instantiation of solver objects¶

First, we create two different solver object for solving the two inversion problem stated above.

# Create LCG solver
LCGsolver = o.CG(o.BasicStopper(niter=2000))
LCGsolver.setDefaults(save_obj=True)

# Create LCG solver for symmetric systems
SLCG = o.CGsym(o.BasicStopper(niter=2000))
SLCG.setDefaults(save_obj=True, save_model=True)

Secondly, we run the solvers to minimize the objective functions previously defined.

LCGsolver.run(Phi1, verbose=True)
SLCG.run(Phi2, verbose=True)
fig, ax = plt.subplots()
ax.semilogy((LCGsolver.obj / LCGsolver.obj[0]))
ax.set_xlabel("iteration")
ax.set_xlim(0)
plt.suptitle(r"$\phi_{LCG}$")
plt.show()
fig, ax = plt.subplots()
ax.plot(SLCG.obj)
ax.set_xlabel("iteration")
ax.set_xlim(0), ax.set_ylim(0,10)
plt.suptitle(r"$\phi_{symLCG}$")
plt.show()

Also, let's compare the two inverted functions with the analytical solution. To find the solution for the continuous case we need three conditions:

$$$\frac{d^2f(x)}{dx^2}=1,\\ f(x=0)=0,\\ f(x=x_f)=0.$$$

$x = 0$ and $x = x_f$ are not sampled and lay outside of the interval $\mathbf{x}$.

X = np.linspace(dx, N*dx, N)
alpha = 0.5
beta  = -(X[-1] + dx) * 0.5
gamma = 0.0
f_an  = alpha * X * X + beta * X + gamma
fig, ax = plt.subplots()

ax.plot(X, Phi1.model.plot(), linewidth=3, label="$\hat{f}(x)$ from $\phi_{LCG}$")
ax.plot(X, Phi2.model.plot(), linewidth=2, dashes=[6, 2], label="$\hat{f}(x)$ from $\phi_{SymLCG}$")
ax.plot(X, f_an, 'k', linewidth=2, dashes=[2, 2], label="Analytical solution")

ax.legend()
ax.set_xlabel(r"$x$"), ax.set_ylabel("$f(x)$")
ax.set_ylim(-0.14, 0), ax.set_xlim(0,1)

plt.suptitle("Inverted functions")
plt.show()

Now, let's try to solve both inversions using the inverse of $D_2$ as a preconditioner

# P = [D_2]^-1
P = o.Matrix(o.VectorNumpy(np.linalg.inv(D2.matrix.getNdArray())), f, y)

Phi1Prec = o.LeastSquares(f.clone(), y, D2, prec=P * P)
Phi2Prec = o.LeastSquaresSymmetric(f.clone(), y, D2, prec=P)
LCGsolver.setDefaults() # Re-setting default solver values
SLCG.setDefaults() # Re-setting default solver values
LCGsolver.run(Phi1Prec, verbose=True)
SLCG.run(Phi2Prec, verbose=True)

As expected, we converge to the global minimum in effectively one iteration.