Optimization basics

Section 0.7 handed us the normal equations: the best least-squares solution to $Ax = b$ is $\\hat x = (A^{T}A)^{-1}A^{T}b$. That is a closed-form answer to a linear problem. FWI and modern seismic inversion are non-linear: the forward operator is a wave-equation simulator, not a matrix, and there is no closed-form solution. What replaces the matrix inverse is iterative optimization.

1. The cost function

A cost function $f(x)$ measures how bad a proposed solution $x$ is. The optimizer’s job is to find the $x$ that makes $f(x)$ as small as possible. For seismic:

• Least-squares inversion: $f(x) = \\|Ax - b\\|\_2^2$. Convex, so every local minimum is global; the minimum is unique when $A$ has full column rank. Closed-form solution available.
• FWI: $f(m) = \\tfrac{1}{2}\\|\\,\\mathcal{F}(m) - d\\,\\|^{2}$, where $m$ is the velocity model and $\\mathcal{F}(m)$ is a wave-equation forward model. It is not convex in general and has many local minima.
• Regularized inversion: $f(x) = \\|Ax - b\\|\_2^2 + \\epsilon\\,\\|Dx\\|\_2^2$, where $D$ is a roughness operator. With $D = I$ it shrinks the solution norm; with $D$ a derivative it smooths the solution.

2. Gradient descent

Take a step in the direction opposite the gradient. The step size $\\alpha$ (also called the learning rate) controls how far. Too small and you crawl; too large and you overshoot or even diverge. Figure 0.9 runs this update on four costs over two model parameters, $m\_1$ and $m\_2$: (a) draws the path on the cost surface, (b) the cost above the minimum at every iteration, and (c) how much one step multiplies the error along each direction of the surface. Drag the start point anywhere in (a), or set it with the sliders, and change the step. The figure opens on the valley $f = m\_1^2 + 5m\_2^2$ with $\\alpha = 0.19$.

At that opening step every iteration overshoots the steep $m\_2$ direction, so the path zigzags across the valley, yet the error still shrinks by a factor of 0.90 per step and descent converges in 96 iterations. Lower $\\alpha$ to the best fixed step, $\\alpha^{\\ast} = 1/6$, which is 0.166 on the slider, and it converges in 25. Then try:

• Round bowl, $f = m\_1^2 + m\_2^2$, curvature 2 in every direction. At $\\alpha = 0.15$ the path slides smoothly to the origin. At $\\alpha = 0.5$ descent lands on the minimum in one step. Past 0.5 each step overshoots the minimum and the path zigzags across it (at 0.8 it still converges, in 18 iterations); at exactly 1 it bounces between two mirror points forever; beyond 1 every step overshoots by more than it corrects and the iterates diverge.
• Valley and narrow valley, $\\kappa = 5$ and 25: the level sets are long ellipses. Any step with $\\alpha\\lambda\_{\\max} > 1$ (on the valley, $\\alpha > 0.1$) overshoots across the narrow direction, and any step small enough to avoid that crawls along the long one. On the narrow valley even the best fixed step on the slider, 0.0384, needs 141 iterations; momentum, $\\beta = 0.45$ with $\\alpha = 0.0556$, brings that down to 36.
• Two basins: at $\\alpha = 0.15$, from the start (2.20, 2.40) descent reaches the global minimum at (1.27, 0.45) in 18 iterations; from (−1.80, 1.60), on the other side of the dashed ridge through the saddle, it stops in the local minimum at (−1.24, −0.53), 0.41 higher. Gradient descent falls into whichever minimum is downhill from its starting point; it has no way to see the other one. Finding the global minimum of a non-convex function has no efficient general guarantee, and for FWI’s millions of parameters exhaustive search is unaffordable; FWI uses frequency continuation and careful initialization instead.

3. Three pathologies the figure exposes

• Slow convergence. A step too small for the flat direction: a slow, monotone approach. Fix: increase $\\alpha$, or use accelerated methods (momentum, conjugate gradients).
• Zigzag and divergence. A step too large for the steepest direction ($\\alpha\\lambda\_{\\max} > 1$ zigzags, $\\alpha\\lambda\_{\\max} > 2$ diverges), which on an ill-conditioned surface is the only way to make progress along the flat direction. Fix: preconditioning, or quasi-Newton methods that use curvature information, so that one step suits every direction.
• Local minimum. A non-convex surface. Fix: multi-start, annealing, frequency continuation in FWI, or convex relaxations.

4. The Lipschitz / condition-number connection

For a quadratic $f(x) = \\tfrac{1}{2}x^{T}Ax$ the Hessian is $A$, and along its eigenvector with eigenvalue $\\lambda\_i$ each step multiplies the error by $1 - \\alpha\\lambda\_i$. Descent therefore diverges once $\\alpha > 2/\\lambda\_{\\max}$, the optimal fixed step is $\\alpha^{\\ast} = 2/(\\lambda\_{\\max} + \\lambda\_{\\min})$, and the convergence rate depends on the condition number

If $\\kappa = 1$ and $\\alpha = 1/\\lambda$, gradient descent converges in one step. At the optimal step the error shrinks by a factor $(\\kappa - 1)/(\\kappa + 1)$ per iteration, so $\\kappa = 100$ needs about 50 iterations for each factor of $e$. The valley in Figure 0.9 is $f = m\_1^2 + 5m\_2^2$, with Hessian eigenvalues 2 and 10, so $\\kappa = 5$ and the best factor is $2/3$. Plate (c) draws both multipliers, $|1 - \\alpha\\lambda\_{\\min}|$ and $|1 - \\alpha\\lambda\_{\\max}|$, against $\\alpha$ and shades the larger of the two, the one that sets the pace: the best step sits at the lowest point of that shaded V, where the two curves cross, and the steps that diverge are hatched.

Every practical inversion is ill-conditioned: the Jacobian’s singular values span many orders of magnitude. The solutions people actually use are conjugate gradients, preconditioning, and quasi-Newton updates like L-BFGS. These are all in the same family as gradient descent; they spend more memory or compute per step to bend the trajectory toward the minimum. Momentum in Figure 0.9 is the smallest example: it cuts the multiplier at the best step from $(\\kappa - 1)/(\\kappa + 1)$ to about $(\\sqrt{\\kappa} - 1)/(\\sqrt{\\kappa} + 1)$.

5. Regularization: adding a penalty

If the cost has many equally good solutions, or the gradient is unstable near the optimum, add a penalty that breaks ties toward desirable solutions:

• Tikhonov (L2): $+\\epsilon\\|x\\|\_2^{2}$ shrinks the norm toward zero. It adds $\\epsilon I$ to $A^{T}A$, which bounds the smallest eigenvalue away from zero and stabilizes the inversion: every eigenvalue rises by $\\epsilon$, so $\\kappa$ falls to $(\\lambda\_{\\max} + \\epsilon)/(\\lambda\_{\\min} + \\epsilon)$.
• Smoothing: $+\\epsilon\\|Dx\\|\_2^{2}$, where $D$ is a discrete derivative. Prefers smooth solutions. Tomographic velocity inversion uses this heavily.
• Sparsity (L1): $+\\epsilon\\|x\\|\_{1}$ prefers sparse solutions with few non-zero entries. Compressed-sensing reconstruction lives here.

6. From one iteration to FWI

FWI is the outer loop above the figure. Every FWI iteration:

1. Run a wave-equation simulator forward to get predicted data.
2. Subtract from observed data to get the data residual.
3. Back-project the residual via the adjoint-state method to get a gradient in model space.
4. Take a gradient step: $m\_{k+1} = m\_k - \\alpha\\, g\_k$.
5. Check the stopping criterion; loop.

One FWI iteration on a modern 3D dataset can take thousands of CPU hours. The efficiency of the outer optimizer (gradient vs. L-BFGS vs. truncated-Newton) directly determines how many iterations you can afford. The basins of FWI’s misfit are cycle skips: a starting model whose traveltimes are off by more than half a period descends to the wrong wiggle, exactly as a start on the wrong side of the ridge does in Figure 0.9.

Where this goes next

Part 0 closes with Section 0.10: the wave equation. The physics behind the simulator inside FWI, the propagation model that migration inverts, and the starting point for modern imaging. Thirty minutes of PDEs, then we are out of prerequisites and into processing proper.

References

• Strang, G. (2016). Introduction to Linear Algebra (5th ed.). Wellesley-Cambridge.
• Tarantola, A. (1984). Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49, 1259.
• Virieux, J., Operto, S. (2009). An overview of full-waveform inversion in exploration geophysics. Geophysics, 74, WCC1.
• Pratt, R. G. (1999). Seismic waveform inversion in the frequency domain, Part 1. Geophysics, 64, 888.