Model fitting by least squares

Analytical solutions

We have found a numerical solution to fitting applications, such as:

$${\bf0} \quad \approx \quad \nabla \tau \ - \ {\bf d}$$ (92)

An analytical solution is much faster. From any regression, we get the least squares solution when we multiply by the transpose of the operator. Thus,

$${\bf0} \eq \nabla\T\;\nabla \tau \ -\ \nabla\T\,{\bf d}$$ (93)

We need to understand what is the transpose of the gradient operator. Recall the finite difference representation of a derivative in Chapter 1. Ignoring end effects, the transpose of a derivative is the negative of a derivative. Because the transpose of a column vector is a row vector, the adjoint of a gradient $\nabla$, namely, $\nabla\T$ is more commonly known as the vector divergence ($\nabla \cdot$). Likewise, $\nabla\T\,\nabla$ is a positive definite matrix, the negative of the Laplacian $\nabla^2$. Thus, in more conventional mathematical notation, the solution $\tau$ is that of Poisson's equation.

$$\nabla^2 \tau \eq - \ \nabla \cdot {\bf d}$$ (94)

In the Fourier domain, we can have an analytic solution. There, $-\nabla^2 = k_x^2 + k_y^2$ where $(k_x , k_y)$ are the Fourier frequencies on the $(x , y )$ axes. Instead of thinking of equation (94) as a convolution in physical space, think of it as a product in Fourier space. Thus, the analytic solution is:

$$\tau(x,y) \eq {\bf FT}^{-1}\ \frac{ {\bf FT} \ \ \nabla \cdot {\bf d} } { k_x^2+k_y^2}$$ (95)

where FT denotes two-dimensional Fourier transform over $x$ and $y$ . Here is a trick from numerical analysis that gives better results: Instead of representing the denominator $k_x^2+k_y^2$ in the most obvious way, let us represent it in a manner consistent with the finite-difference way we expressed the numerator $\nabla \cdot {\bf d}$. Recall that $- i\omega \Delta t \approx i \hat{\omega} \Delta t = 1-Z = 1-{\rm exp}(-i \omega \Delta t)$, which is a Fourier domain way of saying that difference equations tend to differential equations at low frequencies. Likewise, a symmetric second time derivative has a finite-difference representation proportional to $(-2+Z+1/Z)$ and in a two-dimensional space, a finite-difference representation of the Laplacian operator is proportional to $(-4+X+1/X+Y+1/Y)$, where $X = {\rm exp}(i k_x \Delta x )$ and $Y = {\rm exp}(i k_y \Delta y )$. Fourier solutions have peculiarities (periodic boundary conditions) that are not always appropriate in practice, but having these solutions available is often a nice place to start from when solving an application that cannot be solved in Fourier space.

For example, suppose we feel some data values are bad, and we would like to throw out the regression equations involving the bad data points. At Vesuvius, we might consider the strength of the radar return (which we have previously ignored) and use it as a weighting function ${\bf W}$. Now, our regression (92) becomes:

$${\bf0} \quad \approx\quad {\bf W} \ (\nabla \phi -{\bf d}) \eq ({\bf W} \nabla) \phi \ -\ {\bf Wd}$$ (96)

which is a regression with an operator ${\bf W \nabla}$ and data ${\bf Wd}$. The weighted problem is not solvable in the Fourier domain, because the operator $({\bf W} \nabla )\T\; {\bf W} \nabla$ has no simple expression in the Fourier domain. Thus, we would use the analytic solution to the unweighted problem as a starting guess for the iterative solution to the real problem.

With the Vesuvius data, we could construct a weight ${\bf W}$ from the signal strength. We also have available the curl, which should vanish. Vanishing is an indicator of questionable data that should be weighted down relative to other data.

Model fitting by least squares