Project 3: Adaptive Splines

Due: 2022-04-22

A polyharmonic spline is a type of function approximator for functions from ${\mathbb{R}}^d\to \mathbb{R}$ of the form

$$s(x)=\sum _{j=1}^m\varphi (\parallel x-u_j\parallel )a_j+{\gamma }_{1:d}^{T}x+{\gamma }_{d+1}$$

where

$$\varphi (\rho )=\{\begin{array}{ll}{\rho }^k\log (\rho ),& k\text{ even},\\ {\rho }^k,& k\text{ odd}.\end{array}$$

Common examples include cubic splines $(k=3)$ and thin plate splines ($k=2)$. We will consider cubic splines here.

We will touch on two optimization questions in this project. First, how do we use splines to help optimization? Second, how do we best choose the coefficients $a$ and $\gamma$ and the centers $\{u_j{\}}_{j=1}^m$?

Code setup

In the interest of making all our lives easier, we have a fair amount of starter code for this project. You are not responsible for all of it in detail, but I certainly recommend at least becoming familiar with the interfaces!

Basis functions and evaluation

As a point of comparison, we will want to look at the ordinary spline procedure (where there are as many centers as there are data points, and data is given at the centers). To start with, we need $\varphi$ and its derivatives.

When evaluating the spline will be convenient to refer to $a$ and $\gamma$ together in one vector, which we will call $c$.

We are frequently going to care about the residual on our test function (over many evaluation points), so we also provide small helpers for that computation and for plotting the spline and the residual.

Fitting the spline

We consider two methods for fitting the spline. For both, we need to compute two matrices:

• The kernel matrix $K_{XU}$ where $(K_{XU}{)}_{ij}=\varphi (\parallel x_i-u_j\parallel )$

• The polynomial part ${\mathrm{\Pi }}_X$ consisting of rows of the form $\left[\begin{array}{cc}{x}_i^T& 1\end{array}\right]$

Later, we will also have need for the Jacobian of $K_{XU}$ with respect to the components of $u$.

Standard spline fitting

In the standard scheme, we have one center at each data point. Fitting the spline involves solving the system

$$\left[\begin{array}{cc}K+{\sigma }^{2}I& \mathrm{\Pi }\\ {\mathrm{\Pi }}^T& 0\end{array}\right]\left[\begin{array}{c}a\\ \gamma \end{array}\right]=\left[\begin{array}{c}y\\ 0\end{array}\right].$$

The first $n$ equations are interpolation conditions. The last $d+1$ equations are sometimes called the "discrete orthogonality" conditions. When the regularization parameter ${\sigma }^2$ is nonzero, we have a smoothing spline.

Least squares fitting

When the number of data points is large, we may want to use a smaller number of centers as the basis for our approximation. That is, we consider an approximator of the form

$$s(x)=\sum _{j=1}^n\varphi (\parallel x-u_j\parallel )a_j+{\gamma }_{1:d}^{T}x+{\gamma }_{d+1}$$

where the locations $u_j$ are decoupled from the $x_i$. We determine the coefficients $a$ and $\gamma$ by a least squares problem

$$\text{minimize }{\parallel \left[\begin{array}{cc}{K}_{XU}& \mathrm{\Pi }\end{array}\right]\left[\begin{array}{c}a\\ \gamma \end{array}\right]-y\parallel }^2+{\sigma }^2\parallel a{\parallel }^2$$

Sampling strategies

Splines are fit to sampled data, so we want a number of ways to draw samples. Our example is 2D, so we will stick to 2D samplers.

We start with samplers on a regular mesh of points $(x_i,y_j)$ (sometimes called a tensor product mesh). These get big very fast in high-dimensional spaces, but are fine in 2D.

In higher-dimensional spaces, it is tempting to choose random samples. But taking independent uniform draws is not an especially effective way of covering a space – random numbers tend to clump up. For this reason, low discrepancy sequences are often a better basis for sampling than (pseudo)random draws. There are many such generators; we use a relatively simple one based on an additive recurrence with a multiplier based on the "generalized golden ratio".

We will provide both random and quasi-random samplers.

Optimizers

Most of the time if you are using standard optimizers like Levenberg-Marquardt or Newton, you should use someone else's code and save your ingenuity for problem formulation and initial guesses. In this spirit, we provide two solvers that we will see in class in a couple weeks: a Levenberg-Marquardt implementation and a Newton implementation.

We will see the Levenberg-Marquardt algorithm for nonlinear least squares as part of our regular lecture material, but we also include it here. This version uses logic for dynamically scaling the damping parameter (courtesy Nocedal and Wright).

The trust-region Newton solver optimizes a quadratic model over some region in which it is deemed trustworthy. The model does not have to be convex!

Incremental QR and forward selection

The forward selection algorithm is a factor selection process for least squares problems in which factors are added to a model in a greedy fashion in order to minimize the residual at each step.

For fast implementation, it will be convenient for us to be able to extend a QR factorization by adding additional columns to an existing factorization. This is easy enough to do in practice, but requires a little poking around in the guts of the Julia QR code.

We also will want to compute residuals for many of our algorithms.

Some sanity checks on a random test case (keep first 5, choose 5 more):

• 1 2 3 4 5 11 12 20 17 6

  selected

• Does A[:,I] = QR? Relerr 3.733265243763775e-16

• Does r = b-A[:,I]*x? Relerr 5.105187043779148e-16

• Does A[:,I]'*r = 0? Relerr 1.2869532583151618e-17

Test function

The Himmelblau function is a standard multi-modal test function used to explore the properties of optimization methods. We will take the log of the Himmelblau function (shifted up a little in order to avoid log-of-zero issues) as our running example.

The Himmelblau function has a local maximum at $(-0.270845,-0.923039)$ and four global minimizers:

$$\begin{array}{rl}g(3.0,2.0)& =0.0,\\ g(-2.805118,3.131312)& =0.0,\\ g(-3.779310,-3.283186)& =0.0,\\ g(3.584428,-1.848126)& =0.0.\end{array}$$

Tasks

1. Express the Himmelblau function as a nonlinear least squares function and use the provided Levenberg-Marquardt solver to find the minimum at $(3,2)$ from a starting guess at $(2,2)$. Use the monitor function to produce a semilog convergence plot.

2. Run ten steps of an (unguarded) Newton solver for the Himmelblau function. It may save you some algebra to remember the relation between the Hessian of a nonlinear least squares function and the Gauss-Newton approximation. What does it converge to from the starting guess $(2,2)$? Plot convergence (gradient norm vs iteration) again.

3. Run the same experiment with the Newton solver, but this time using the trust region version that is provided.

Optimizing on a surrogate

Let's now fit a surrogate to the Himmelblau function sampled on 50 samples drawn from our quasi-random number generator.

Does the spline have minima close to the minima of the true function? The answer depends how "nice" the minimizers are and how good the approximation error is. Concretely, suppose $|s(x)-g(x)|<\delta$ over $x\in \mathrm{\Omega }$ and $x_{\ast }$ is a minimizer of $g$. If $\underset{\parallel u\parallel =1}{min}g(x_{\ast }+\rho u)-g(x_{\ast })>2\delta$ and the ball of radius $\rho$ around $x_{\ast }$ lies entirely inside $\mathrm{\Omega }$, then $\underset{\parallel u\parallel =1}{min}s(x_{\ast }+\rho u)-s(x_{\ast })>0$. Assuming continuity, this is enough to guarantee that $s$ has a minimizer in the unit ball of radius $\rho$ around $x_{\ast }$. If $g$ has a Lipschitz continuous Hessian with Lipschitz constant $M$ and $\rho <{\lambda }_{min}(H_g(x_{\ast }))/M$, then a sufficient condition is that $\frac{1}{2}{\rho }^2({\lambda }_{min}(H_g(x_{\ast }))-M\rho )>2\delta$ (still assuming the ball of radius $\rho$ lies within the domain $\mathrm{\Omega }$).

(We can get better bounds if we can control the error in the derivative approximation, but we will leave this alone in the current assignment.)

Putting aside theoretical considerations for a moment, let's try to find the minimizer (or at least a local minimizer) of the spline approximation $s(x)$. We can make a first pass at this by sampling on a mesh and taking the point where we see the smallest function value.

Tasks

1. Let $x_{\ast }$ be a strong local minimizer of $g$, where $g$ is twice continuously differentiable and the Hessian $H_g$ has a Lipschitz constant $M$ in the operator 2-norm, and that $\rho <{\lambda }_{min}(H_g(x_{\ast }))/M$. Using a Taylor expansion about $x_{\ast }$, show that $\frac{1}{2}{\rho }^2({\lambda }_{min}(H_g(x_{\ast }))-M\rho )>2\delta$ is sufficient to guarantee that $\underset{\parallel u\parallel =1}{min}g(x_{\ast }+\rho u)-g(x_{\ast })>2\delta$.

2. Compute the smallest singular value of the Hessian of the log-transformed Himmelblau function at $(3,2)$. Also compute the approximation error in the spline fit at $(3,2)$ and use that as an approximate $\delta$ and plug into the estimate

$$\tilde{\rho }=2\sqrt{\frac{\delta }{{\lambda }_{min}(H_g(x_{\ast }))}}$$

3. Take ten steps of Newton to find minima the minimum of the spline close to the point found above. Plot the gradient norm on a semilog scale to verify quadratic convergence. What happens with a starting guess of $(3,2)$?

For the first task, you may use without proof that ${\lambda }_{min}(A+E)\ge {\lambda }_{min}(A)-\parallel E{\parallel }_2$ when $A$ and $E$ are symmetric. For the second question, remember that the log transform we used is $g\mapsto \log (g+10)$. You should also not expect to get a tiny value for $\tilde{\rho }$!

Non-adaptive sampling strategies

If we want to use a spline as a surrogate for optimization, we need it to be sufficiently accurate, where what is "sufficient" is determined by the behavior of the Hessian near the desired minima. For the rest of this assignment, we focus on making more accurate splines.

We start by defining a highly refined "ground truth" mesh of ${10}^4$ points laid out in a mesh. Our figure of merit will be the root mean squared approximation error. We also define a coarser "data mesh" that we will use for least squares spline fits.

We would like to choose a sample set that is going to minimize the RMS error. We will talk about adapting the sample later; for now, let's consider different ways that we could lay out a set of about 50 sample points:

• We could choose a 7-by-7 mesh (uniform or Chebyshev)

• We could choose 50 points at random

• We could choose 50 points from a low-discrepancy sequence (quasi-random)

Let's set up an experimental framework to compare the RMS error on the ground truth mesh for each of these approaches.

Tasks

Run experiments with the test_sample functions for both the standard spline fitting and least squares fitting against the 4000-point data set. Comment on anything you observe from your experiments. Feel free to play with the parameters (e.g. number of points).

Forward stepwise regression

Forward stepwise regression is a factor selection method in which we greedily choose the next factor from some set of candidates based on which candidate will most improve the solution. We have provided you with a reference implementation of the algorithm. Your task is to use this method to choose centers from a candidate set of points in order to improve your least squares fit.

Tasks

Complete the spline_forward_regression function below, then run the code with 500 points chosen from our quasirandom sequence. Run the test harness above to see how much the RMS error measure improves.

Levenberg-Marquardt refinement

In the lecture notes from 4/11, we discuss the idea of variable projection algorithms for least squares. We can phrase the problem of optimizing the centers in this way; that is, we let

$$r=(I-A{A}^{†})\overline{y}$$

where

$$A=\left[\begin{array}{cc}{K}_{XU}& {\mathrm{\Pi }}_X\\ \sigma I& 0\end{array}\right],\overline{y}=\left[\begin{array}{c}y\\ 0\end{array}\right]$$

is treated as a function of the center locations. The regularization parameter $\sigma$ plays a critical role in this particular problem, so we are careful to keep it throughout!

To save us all some trouble, I will provide the code to compute $r$ and the Jacobian with respect to the center coordinates.

As usual, we need a finite difference check to prevent programming errors.

• Finite difference relerr in Jacobian computation: 1.6565625055732684e-7

Tasks

Use the provided Levenberg-Marquardt solver to refine the center locations picked by the forward selection algorithm. This is a tricky optimization, and you may need to fiddle with solver parameters to get something you are happy about. As usual, you should also provide a convergence plot.

Newton refinement

Levenberg-Marquardt is popular for nonlinear least squares problems because it only requires first derivatives and often still gives very robust convergence. But the convergence behavior for this problem is disappointing, and we have all the derivatives we need if we only have the fortitude to use them.

As in the Levenberg-Marquardt case, I will not require you to write the function evaluation and derivatives.

Finite difference checks:

• For the gradient: 7.256850871224047e-8

• For the Hessian: 2.4552525270252068e-8