Linear algebra primer

Beyond filtering, almost every processing step is linear algebra wearing a costume. Migration applies the adjoint of a modelling operator. Tomography is least squares. Inversion is a regularized solve. FWI is nonlinear optimization whose every step is linear algebra: an adjoint-state gradient and, in Gauss-Newton, an inner regularized solve. If you cannot read $Ax = b$ as a picture, the algorithms will look like black magic.

1. A matrix is a function

Forget rows and columns for a moment. A $2\\times 2$ matrix M = \[\\,v\_1\\;\\, v\_2\\,\] is a rule that takes a vector $x = (x\_1, x\_2)$ and returns $x\_1 v\_1 + x\_2 v\_2$. The matrix is its two columns; everything else is bookkeeping. In Figure 0.7 you set those two columns by hand and watch the whole plane follow, then run the matrix backwards to solve $Mx = b$ for noisy data.

The two ringed arrows in (a) are $M\\hat\\imath$ and $M\\hat\\jmath$, the images of $\\hat\\imath = (1, 0)$ and $\\hat\\jmath = (0, 1)$, and they are the columns of the matrix: $M\\hat\\imath = (m\_{11}, m\_{21})$ and $M\\hat\\jmath = (m\_{12}, m\_{22})$. Everything else in (a) follows from them. The pale grid is carried to the darker one, the dashed unit square to the shaded parallelogram, and the dashed unit circle to the solid ellipse. The figure opens on $M\\hat\\imath = (1.50, 0.30)$ and $M\\hat\\jmath = (1.00, 0.80)$.

2. Determinant

For M = \\begin{pmatrix} m\_{11} & m\_{12} \\\\ m\_{21} & m\_{22} \\end{pmatrix} the determinant $\\det M = m\_{11}m\_{22} - m\_{12}m\_{21}$ is the signed area of the transformed unit square. For the starting matrix it is $1.50 \\cdot 0.80 - 1.00 \\cdot 0.30 = 0.90$, the area of the shaded parallelogram in (a). Consequences:

• $|\\det M|$ is how much the matrix scales areas (in 2D) or volumes (in 3D).
• The sign of $\\det M$ is positive if the matrix preserves orientation, negative if it reflects. The small arc at the origin in (a) turns from $M\\hat\\imath$ to $M\\hat\\jmath$: anticlockwise while $\\det M > 0$, clockwise once the plane is flipped.
• $\\det M = 0$ means the image collapses to a line (or point). The matrix is singular and cannot be inverted. Lower $m\_{22}$ to 0.20, or choose Singular, rank 1 under Start from: the columns line up and the ellipse in (a) flattens to a segment. Keep going to $m\_{22} = -0.40$ and the determinant comes back as $-0.90$, the same area with the plane flipped over.

3. Eigenvalues and eigenvectors

An eigenvector $v$ of $M$ is a nonzero direction that survives the transformation up to scaling: $Mv = \\lambda v$. The scalar $\\lambda$ is the eigenvalue. The eigenvectors are the directions along which the matrix is pure stretching; the eigenvalues tell you by how much (and a negative $\\lambda$ also flips the direction). In (a) the real eigen-directions are dotted: the starting matrix has $\\lambda = 1.80$ and $0.50$, and each dotted line is carried onto itself.

For a $2\\times 2$ matrix, the eigenvalues come from the characteristic polynomial:

so $\\lambda = \\tfrac{1}{2}\\bigl(\\operatorname{tr} M \\pm \\sqrt{(\\operatorname{tr} M)^{2} - 4\\det M}\\bigr)$. When the discriminant is negative the eigenvalues are a complex pair $\\lambda = r e^{\\pm i\\theta}$: no real direction survives, and the matrix acts like a rotation by $\\theta$ combined with a scaling by $r = \\sqrt{\\det M}$, up to a change of basis. The rotation by 36.9 degrees under Start from has $\\lambda = 0.80 \\pm 0.60i$, and no dotted line appears.

4. Singular values and conditioning

The singular values $\\sigma\_1 \\ge \\sigma\_2 \\ge \\dots$ are the square roots of the eigenvalues of $M^{T}M$. They are the semi-axis lengths of the ellipse the matrix makes from the unit circle: 1.94 and 0.46 for the starting matrix, the two solid half-axes in (a). They are not the eigenvalues unless $M$ is symmetric; choose the symmetric matrix under Start from and the dotted lines fall on the ellipse’s axes, with $\\sigma\_i = |\\lambda\_i|$ = 3 and 1. The condition number is

If $\\kappa$ is large, the inverse problem is ill-conditioned: the relative error in the solution $x$ can be up to $\\kappa$ times the relative error in the data $b$. The absolute error is simpler to picture. Running $M$ backwards divides each component of a data error by its singular value, so an error of size $|\\delta b|$ comes back up to $|\\delta b|/\\sigma\_2$ long. Plate (b) shows it: 48 noisy copies of $b$ in a disc of radius 0.20 come back as a needle-shaped ellipse whose longest radius is 2.2 times the disc’s, and which is $\\kappa = 4.2$ times longer than it is wide. Lower $m\_{22}$ to 0.25 and the columns nearly line up: $\\sigma\_2$ falls to 0.041, $\\kappa$ rises to 45 and the noise gain to 25 times. Every inversion we will meet later has a condition number, and regularization is the tool for taming it.

5. Systems of equations: $Ax = b$

Given a known mapping $A$ and an observation $b$, find the unknown $x$. Three cases:

• Square, invertible. $x = A^{-1} b$. Rare in the real world; never invert explicitly, solve instead (LU in general, QR for stability, Cholesky when $A$ is symmetric positive definite).
• Overdetermined (more equations than unknowns). No exact solution. Use least-squares: minimize $\\|Ax - b\\|\_2^2$. The solution is the normal equation $A^{T}A\\,x = A^{T}b$. This is the basic form of tomography, residual statics, and inversion.
• Underdetermined (fewer equations than unknowns). Infinitely many solutions; you must pick one by adding a penalty like $\\|x\\|\_2^2$ (Tikhonov) or $\\|x\\|\_1$ (sparsity). Least-squares migration is often underdetermined; conventional migration applies the adjoint $A^{T}$ instead of an inverse, which acts as an implicit regularization.

6. The least-squares normal equations in one line

That one equation sits underneath residual statics solves, surface-consistent amplitude and deconvolution decomposition, traveltime tomography, model-based inversion, and the inner loop of every iterative imaging algorithm. When $A^{T}A$ is nearly singular, you regularize:

The term $\\epsilon I$ adds $\\epsilon$ to every eigenvalue $\\sigma\_i^2$ of $A^{T}A$, lifting the small ones off the floor so the inverse exists. The regularized solve keeps the share $f\_i = \\sigma\_i^2/(\\sigma\_i^2 + \\epsilon)$ of the exact solution’s component along each direction, one half where $\\epsilon = \\sigma\_i^2$, so it switches off the weak directions first: that is plate (c). Start from the near-singular matrix and raise $\\epsilon$ from 0.000001 to 0.01. The noise gain falls from 25 times to 3.5 times, but the solution moves 5.87 away from the exact one, toward the origin: less noise, more bias. Tuning $\\epsilon$ is the eternal art of inversion.

7. Why this is the workhorse for processing

• Tomographic velocity inversion. $A$ is the ray-path matrix (entry $A\_{ij}$ is the length of ray $i$ in cell $j$); $x$ is the slowness update; $b$ is the traveltime residuals. Least squares with smoothing regularization, solved iteratively (CGLS or LSQR) because $A$ is huge but sparse; $A^{T}A$ is never formed.
• Post-stack inversion for impedance. $A$ is the wavelet-convolution matrix applied after a derivative, $A = \\tfrac12 WD$; $x$ is the log-impedance $\\ln Z$; $b$ is the observed seismic trace. Same equations, different names.
• FWI. Each outer iteration linearizes the misfit: the adjoint-state method gets the gradient $J^{T}r$ from one forward and one adjoint wave simulation, and a Gauss-Newton step solves $(J^{T}J + \\epsilon I)\\,\\Delta m = -J^{T}r$. Linear algebra is in the hot core; calculus on matrices is the outer shell.
• Migration (LS-migration). Matching observed data to predicted data from a reflectivity model, least-squares on a linear forward operator.

Where this goes next

Section 0.8 brings in the other piece of every inverse problem: a model of what the noise does. You cannot set the right $\\epsilon$ without a probabilistic framing, so random variables and noise statistics come next.

References

• Strang, G. (2016). Introduction to Linear Algebra (5th ed.). Wellesley-Cambridge.
• Claerbout, J. F. (1976). Fundamentals of Geophysical Data Processing. McGraw-Hill.
• Yilmaz, Ö. (2001). Seismic Data Analysis (2 vols.). SEG.
• Tarantola, A. (1984). Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49, 1259.