LW1: Least-Squares Problems

1. Polynomial regression¶

This first exercise will permit to review basic instructions in Python as well as standard functions from NumPy and Matplotlib. The goal is to estimate the parameters of a polynomial model given experimental (i.e. noisy) data using the least-squares method. We’ll first define the true model, generate some noisy synthetic data, and then estimate the model parameters.

Consider the ground truth polynomial model $y_{\text{true}}(t) = 10t^3 + 3t^2 + 2t + 1$ for $t\in \mathbb{R}$.

Questions

1. Represent the function $y_{\text{true}}(t)$ on $[-1, 1]$.

2. Given $N =10$ equispaced time instants $t_1, \ldots, t_N$ over $[-1, 1]$, generate a noisy vector of measurements $\mathbf{y} \in \mathbb{R}^N$ with entries $y_n = y(t_n) = y_{\text{true}}(t) + \epsilon$, where $\epsilon \sim \mathcal{N}(0, 1)$ is Gaussian. Represent noisy measurements and the ground truth model on the same graph.

3. From the data $\mathbf{y}$, we wish to fit a polynomial model of order $d < 10$ of the form

   $$y_{\text{model}}(t) = x_0 + x_1 t + x_2 t^2 + \ldots + x_{d+1}t^d$$

   Write the linear model matrix $\mathbf{A}$ for an arbitrary order $d$. Construct this matrix numerically, given $d$ and the vector of time instants $\mathbf{t} = [t_1, t_2, \ldots, t_N]$. (Hint: check the np.vander function).

4. Given an integer $d <10$, write and solve the associated least squares problem. Display the estimated vector of coefficients. Compare the solution obtained by the explicit expression and that obtained with the np.linalg.pinv function.

5. Represent the estimated polynomial model for several values of $d$ (say $d=1, 3, 7$). Comment. What happens for $d=9$?

6. (Bonus) The choice of the degree of the polynomial can be evaluated using a standard criterion such as Akaike Information Criterion, defined as $\text{AIC} = 2p + 2\log r^2$, where $p$ is the number of parameters to be estimated and $r$ is the norm of the residual (in this context). Plot this criterion for the different choices of $d$. For each value of $d$ is the minimum attained?

7. (Bonus bis) Play around with differents values of $N$, true polynomial $y_{\text{true}}(t)$, or noise level $\sigma$!

2. Denoising¶

This second problem is the denoising problem which corresponds to the reconstruction of a signal $x(t)$ from noisy measurements $y(t) = x(t) + n(t)$, where $n(t)$ is some noise. A typical assumption is that noise is Gaussian, i.e., for a given time instant $t_i$, $n(t_i) \sim \mathcal{N}(0, \sigma^2)$ where $\sigma^2$ is the noise variance.

The denoising problem can be formulated as a least squares problem. Let $t_1, \ldots, t_N$ be $N$ time instants and denote by $\mathbf{x}, \mathbf{y}$ the vectors encoding the signal and the measurements, respectively. The denoising problem seeks to solve

Preliminary questions

1. What is the solution to the least squares problem $\hat{\mathbf{x}}$?

Usually, we have access to extra information about the signal $\mathbf{x}$: e.g. it is non-negative, it belongs to some class of constraints, it exhibits regularity in some representation domain, etc. A classical assumption is that the signal is smooth. At first order, this is encoded by the function

2. Show that $r(\mathbf{x})$ can be expressed as $r(\mathbf{x}) = \Vert \mathbf{D}_1\mathbf{x}\Vert^2_2$

To take into account the smoothness of the solution into the denoising problem, we formulate a new problem, called regularized least squares or penalized least squares. It reads

where

• $\Vert \mathbf{y}-\mathbf{x}\Vert$ is called the data fidelity term;
• $r(\mathbf{x})$ is the regularization (or penalty) and $\lambda\geq 0$ is called the regularization parameter.

3. Compute the solution to the regularized least squares problem above with $r(\mathbf{x}) = \Vert \mathbf{D}_1\mathbf{x}\Vert^2_2$. Comment on the uniqueness of the solution. What happens when $\lambda \rightarrow 0$? When $\lambda \rightarrow \infty$?

Write here your answers to preliminary questions using Markdown

Implementation tasks

1. Generate a smooth signal $x(t) = \cos(2\pi f_0 t)$ where $f_0 = 5$ Hz and $t \in [0, 1]$ s. Choose a sampling frequency that satisfies the Shannon-Nyquist criterion.
2. Generate noisy measurements by corrupting the signal vector $\mathbf{x}$ by additive white Gaussian noise with SNR = 10 dB. (Recall the formula for the SNR!)
3. Solve the regularized least squares for different values of $\lambda$, and represent corresponding solutions. Compare it to the ground truth signal.
4. A common question in penalized optimization problems is to choose the hyperparameter $\lambda$. A classical approach is to use a heuristic called L-curve: for various values of $\lambda$, plot the data fidelity term vs the penalty term in the cost function. From this tool, what would be a good choice of $\lambda$ according to you? Why?
5. Start again the exercise with the second-order penalty function
   $$r(\mathbf{x}) = \sum_{n=2}^{N-1} (x_{n+1}- 2x_{n} + x_{n-1})^2$$
   Compare your results with the first-order difference penalty.