Model fitting by least squares

Analytical solutions

We have found a numerical solution to fitting applications such as this

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

An analytical solution will be 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}$$ (104)

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. Since 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}$$ (105)

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 (105) 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}$$ (106)

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 their own peculiarities (periodic boundary conditions) which 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  becomes

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

This is a problem we know how to solve, 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 might we construct the weight ${\bf W}$ from the signal strength. We also have available the curl, which should vanish. Its non-vanishing is an indicator of questionable data which could be weighted down relative to other data.

Model fitting by least squares