# Polynomial interpolation and data-fitting

Given a set of * **equally spaced** point* $t_1, t_2, \dots, t_m$, and *data* $y_1, y_2, \dots, y_m$, determine a polynomial of degree $m-1$: $p(t) = x_1 + x_2 t + x_3 t^2 + \cdots + x_m t^{m-1}$, such that

$$p(t_i) = y_i, i = 1, 2, \dots , m$$

This problem can be written in terms of solving the following system of linear equations

$$\left\{ 
\begin{array}{ccccccccccc}
x_1 & + & x_2 t_1 & + & x_3 t_1^2 & + & \cdots  & + & x_m t_1 ^{m-1} & = & y_1 \\
x_1 & + & x_2 t_2 & + & x_3 t_2^2 & + & \cdots  & + & x_m t_2 ^{m-1} & = & y_2 \\
x_1 & + & x_2 t_3 & + & x_3 t_3^2 & + & \cdots  & + & x_m t_3 ^{m-1} & = & y_3 \\
&&&&\vdots& &&&&=& \vdots \\
x_1 & + & x_2 t_m & + & x_3 t_m^2 & + & \cdots  & + & x_m t_m ^{m-1} & = & y_m \\
\end{array}
\right .$$

or in matrix form

$$\left[
\begin{array}{ccccc}
1 & t_1 & t_1^2 & \cdots & t_1^{m-1} \\
1 & t_2 & t_2^2 & \cdots & t_2^{m-1} \\
1 & t_3 & t_3^2 & \cdots & t_3^{m-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & t_m & t_m^2 & \cdots & t_m^{m-1} 
\end{array}
\right]\left[
\begin{array}{c}
x_1 \\
x_2 \\
x_3 \\
\vdots \\
x_m 
\end{array}
\right]
=\left[
\begin{array}{c}
y_1 \\
y_2 \\
y_3 \\
\vdots \\
y_m 
\end{array}
\right]
$$

Here $(t_i,y_i)$ are all given for $i = 1,2,\dots,m$. Only the polynomial coefficients $x_1, x_2, \dots, x_m$ are unknown. The matrix above is obviously a Vandermonde matrix, which is invertible as long as $t_i \neq t_j$, for $i \neq j$. It also implies the uniqueness of the solution of the system.

However, when there're many points to fit, the interpolation polynomial $p(t)$ will have a high degree, and become a ill conditioned problems.

And to fix this, we can seek a polynomial of degree much less than $m-1$, like $n-1$ with $n\ll m$

|$$m=20,n=20$$|$$m=20,n=10$$|
|:-:|:-:|
|||

But it is impossible for a polynomial of degree $n-1$ to satisfy $m$ interpolation conditinos, since to satisfy such equations we need to solve:

$$\left[
\begin{array}{ccccc}
1 & t_1 & t_1^2 & \cdots & t_1^{n-1} \\
1 & t_2 & t_2^2 & \cdots & t_2^{n-1} \\
1 & t_3 & t_3^2 & \cdots & t_3^{n-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & t_m & t_m^2 & \cdots & t_m^{n-1} 
\end{array}
\right]\left[
\begin{array}{c}
x_1 \\
x_2 \\
x_3 \\
\vdots \\
x_n 
\end{array}
\right]
=\left[
\begin{array}{c}
y_1 \\
y_2 \\
y_3 \\
\vdots \\
y_m 
\end{array}
\right]
$$

which may be without a solution. What to do now?$\newcommand{\ffrac}{\displaystyle\frac}$

We don't require the $p(t)$ of degree $n-1$ to interpolate exactly given $m$ points, rather, we need a $p(t)$ such that the difference $\displaystyle \sum_{i=1}^{m} \left|p(t_i) - y_i \right|^2$ is minimized. This is the **linear least squares problem**. In other words, we need to find the polynomial with coefficients $x_1, x_2, \dots, x_n$, that minimize the residual of the system of linear equations. After that $x_1, x_2, \dots, x_n$ are called the **linear squares solution**.

# Linear least squares problem

Given matrix $A_{m \times n}$, $m \geq n$, full rank, and a vector $\vec{b} \in \mathbb{R}^{m}$, so generally there's no solution for equation $A\vec{x} = \vec{b}$. So instead we consider the following linear least squares problem:

$$\min_{\vec{x} \in \mathbb{R}^n}\|\vec{b} - A\vec{x}\|_2^2$$

the solution $\vec{x}^*$ would minimize $\|\vec{b} - A\vec{x}\|_2$. Here's a picture to get the initial idea:



So the residual $\vec{r} = \vec{b} - A \vec{x}^*$ must be perpendicular to the $\DeclareMathOperator*{\range}{range}\DeclareMathOperator*{\rank}{rank} \range(A)$

$Theorem$ 1

Given an $m \times n$ matrix $A$, with $m > n$ and $\rank(A) = n$, and another vector $\vec{b}\in \mathbb{R}^{m}$, then a vector $\vec{x}^* \in \mathbb{R}^n$ minimizes $\|\vec{b} - A\vec{x}\|_2^2$, thereby solving the least squares problem, $iff$ one of these three equivalent conditions are satisfied
1. $\vec{b} - A \vec{x}^* \perp \range(A)$, that is $A^{\mathrm{T}}\left( \vec{b} - A \vec{x}^* \right) = \vec{0} $
2. $A^{\mathrm{T}}A\vec{x}^* = A^{\mathrm{T}}\vec{b}$
3. $A \vec{x}^* $ is the orthogonal projection of $\vec{b}$ onto $\range(A)$.

$Proof$

- First equivalence

>First we prove: *condition 1* $\Longrightarrow \vec{x}^*$ minimizes $\|\vec{b} - A\vec{x}\|_2^2$

>We choose one vector in $\range(A)$, which is $A(\vec{x}^* - \vec{x})$, here $\vec{x}$ is any vector in $\mathbb{R}^n$. So that follow the *condition 1* we have $\vec{b} - A \vec{x}^* \perp A(\vec{x}^* - \vec{x})$. So that we have

>$$\begin{align}
&\|\vec{b} - A\vec{x}\|_2^2 = \left\|\vec{b} - A\vec{x}^* + A\left( \vec{x}^* - \vec{x} \right)\right\|_2^2 \\ \stackrel{\;\;\vec{b} - A \vec{x}^* \perp \,A(\vec{x}^* - \vec{x})}{=} \;\;&\|\vec{b} - A\vec{x}^*\|_2^2 + \left\| A\left( \vec{x}^* - \vec{x} \right) \right\|_2^2 \geq \|\vec{b} - A\vec{x}^*\|_2^2 
\end{align}$$

>Therefore, $\vec{x}^*$ minimizes $\|\vec{b} - A\vec{x}\|_2^2$

>Then we prove: *condition 1* $\Longleftarrow \vec{x}^*$ minimizes $\|\vec{b} - A\vec{x}\|_2^2$

>Assume not, then we can find a different vector $\vec{y} \in \mathbb{R}^n$ such that $\vec{b} - A \vec{y}$ is orthogonal to $\range(A)$, by solving the equation:

>$$A\vec{y} = P \vec{b}$$

>where $P\vec{b}$ represents the orthogonal projection of $\vec{b}$ to $\range(A)$. So that 

>$$\begin{array}{rl}
& \vec{b} - A \vec{y} = \vec{b} - P \vec{b} \perp \range(A) \\[0.5em]
\Rightarrow & \vec{b} - A \vec{y} \neq \vec{b} - A \vec{x}^*  \\[0.5em]
\Rightarrow & A\left( \vec{x}^* - \vec{y} \right) \neq \vec{0}  \\[0.6em]
\Rightarrow & \begin{align*}
\left\|\vec{b} - A\vec{x}^*\right\|_2^2 &= \left\|\vec{b} - A\vec{y} - A\left( \vec{x}^* - \vec{y} \right) \right\| _2^2 \\
&= \left\|\vec{b} - A\vec{y}\right\|_2^2 + \left\|A\left( \vec{x}^* - \vec{y} \right)\right\|_2^2 \\
&> \left\|\vec{b} - A\vec{y}\right\|_2^2
\end{align*}
\end{array}$$

>This contradiction shows that $\vec{b} - A \vec{x}^*$ must be orthogonal $\range(A)$.

- First and Second Equivalence

>It's trivial that they are equivalent:

>$$A^{\mathrm{T}}A\vec{x}^* = A^{\mathrm{T}}\vec{b} \Longleftrightarrow A^{\mathrm{T}}\left( \vec{b} - A \vec{x}^* \right) = \vec{0} 
$$

- Second and Third Equivalence

> Still we can prove that they are equivalent:

> $$\begin{array}{rrl}
& A^{\mathrm{T}}A\vec{x}^* \!\!\!\!\!&= A^{\mathrm{T}}\vec{b} \\[0.4em]
\Leftrightarrow \!\!\!& \vec{x}^* \!\!\!\!\!&= \left( A^{\mathrm{T}}A \right) A^{\mathrm{T}}\vec{b} \\[0.4em]
\Leftrightarrow \!\!\!& A\vec{x}^* \!\!\!\!\!&= A\left( A^{\mathrm{T}}A \right) A^{\mathrm{T}}\vec{b}
\end{array}$$

> Since $A$ is full rank, so that the inverse exists. And we've already known that $A\left( A^{\mathrm{T}}A \right) A^{\mathrm{T}}\vec{b}$ is just the orthogonal projection of $\vec{b}$ onto $\range(A)$.

# Conditioning of linear least squares problems
The condition number of solving a linear system $A\vec{x} = \vec{b}$ is $\DeclareMathOperator*{\Cond}{Cond} \Cond(A) = \big\|A\big\|_2\left\|A^{-1}\right\|_2$, when $A$ is a square matrix.

So what's the condition number of solving the least squares problem $\displaystyle \min_{\vec{x}\in\mathbb{R}^n} \left\| \vec{b} - A \vec{x} \right\|_2$ for a given $m \times n$ matrix $A$, with $m > n$ and $\rank(A) = n$. From the **theorem 1**, we can see that the solution would be 

$$\vec{x}^* = \left( A^{\mathrm{T}}A \right)^{-1} A^{\mathrm{T}} \vec{b}$$

where $A$ is column full rank matrix so that $A^{\mathrm{T}}A$ is invertible. To complete the explanation, now we define the **pseudoinverse** of $A$ as $A^+:=\left( A^{\mathrm{T}}A \right)^{-1}A^{\mathrm{T}}$. So that now we can shorten the expression like

$$\vec{x}^* = \left( A^{\mathrm{T}}A \right)^{-1} A^{\mathrm{T}} \vec{b} = A^+ \vec{b}$$

And in this situation the condition number will be

$$\Cond(A) = \big\|A\big\|_2\left\|A^{+}\right\|_2$$

There's another point: the conditioning of projecting $\vec{b}$ to $\range(A)$.

Define $\theta  := \small\langle \vec{b},\range(A) \small\rangle$ as the angle between $\vec{b}$ and its projection on $\range(A)$, ranging from $-\ffrac{\pi} {2}$ to $\ffrac{\pi} {2}$.
- If $\theta \approx 0$, then a small relative change in $\vec{b}$ will just lead to a similarly small relative change in $P\vec{b}$, well-conditioned!
- If $\theta \approx \ffrac{\pi} {2}$, then a small relative change in $\vec{b}$ will just lead to a relatively large change in $P\vec{b}$, ill-conditioned!

Based on the above obeservation, we have the following theorem

$Thorem$ 2

The condition number of solving the linear least squares problem is bounded by

$$\frac{\Cond(A)} {\cos \theta}$$

$Proof$

$$\begin{cases}
\;\vec{x} = A^{+} \vec{b} \\
\left( \vec{x} + \Delta \vec{x} \right) = A^{+} \left( \vec{b} + \Delta \vec{b} \right) \\
\end{cases} \Longrightarrow \Delta \vec{x} = A^{+} \Delta \vec{b}$$

So that the condition number of solution $\vec{x}$ with respect to the perturbations in $\vec{b}$ is

$$\begin{align}
\kappa_{\vec{b} \to \vec{x}} &= \sup_{\Delta \vec{b} \neq \vec{0}} \frac{\left\| \Delta \vec{x} \right\| / \left\| \vec{x} \right\|} {\left\| \Delta \vec{b} \right\|\big/\left\| \vec{b} \right\|} \\
&= \sup_{\Delta \vec{b} \neq \vec{0}} \frac{\left\| A^{+} \cdot \Delta \vec{b} \right\| \big/ \left\| \vec{x} \right\|} {\left\| \Delta \vec{b} \right\|\big/\left\| \vec{b} \right\|} \\
&= \frac{\left\| A^{+} \right\| \cdot \left\| \vec{b} \right\|} {\left\| \vec{x} \right\|} = \left\| A^{+} \right\| \cdot \frac{\left\| \vec{b} \right\|} {\left\| P\vec{b} \right\|} \cdot \frac{ \left\| P\vec{b} \right\|} {\left\| \vec{x} \right\|} \\
&= \left\| A^{+} \right\| \cdot \frac{1} {\cos \theta} \cdot \frac{ \left\| A\vec{x} \right\|} {\left\| \vec{x} \right\|} \\
&\leq \frac{\left\| A^{+} \right\| \cdot \left\| A \right\|} {\cos \theta} = \frac{\Cond(A)} {\cos \theta}
\end{align}$$

# Algorithms for solving linear least squares problems
## Solve normal equation
From $Theorem$ 1, we know that the solution to the least squares problem satisfies the normal equation

$$A^{\mathrm{T}}A \vec{x} = A^{\mathrm{T}}\vec{b}$$

And since $A$ is full rank, so that $A^{\mathrm{T}}A$ is symmetric positive definite, and it has a Cholesky Factorization

$$A^{\mathrm{T}}A = LL^{\mathrm{T}}$$

1. Form $A^{\mathrm{T}}A$ and $A^{\mathrm{T}}\vec{b}$
2. Factorize $A^{\mathrm{T}}A = LL^{\mathrm{T}}$, using the build-in function `chol`
3. Solve $LL^{\mathrm{T}} \vec{x} = A^{\mathrm{T}} \vec{b}$, by backward/forward substitution

## By QR factorization
Denote the reduced $QR$ factorization of $A$ by $A = \hat{Q} \hat{R}$, where $\hat{Q}$ is a $m \times n$ matrix with orthonormal columns, and $\hat{R}$ is a $n \times n$ upper triangular matrix. Then, from $A^{\mathrm{T}}A \vec{x} = A^{\mathrm{T}}\vec{b}$, so that

$$\hat{R}^{\mathrm{T}} \hat{Q}^{\mathrm{T}} \hat{Q} \hat{R} \vec{x} = \hat{R}^{\mathrm{T}} \hat{Q}^{\mathrm{T}} \vec{b}$$

which implies

$$\hat{R} \vec{x} = \hat{Q}^{\mathrm{T}} \vec{b}$$

1. Compute the $QR$ factorization $A = \hat{Q} \hat{R}$, by Housholder/Modified Gram-Schmidt
2. Computer $\hat{Q}^{\mathrm{T}} \vec{b}$
3. Solve $\hat{R} \vec{x} = \hat{Q}^{\mathrm{T}} \vec{b}$, by backward substitution
***

Another way to derive the equation: Since $A\vec{x}$ is the orthogonal projection of $\vec{b}$ to the $\range(A)$:

$$A\vec{x} = \hat{Q}\hat{Q}^{\mathrm{T}} \vec{b} \Longrightarrow \hat{R} \vec{x} = \hat{Q}^{\mathrm{T}} \vec{b}$$


## By singular value decomposition
The reduced singular value decomposition of a full rank $m \times n$ matrix $A$ can be written as

$$A = \hat{U} \hat{\Sigma} V^{\mathrm{T}}$$

where $\hat{U}$ is an $m \times n$ matrix with orthonormal columns, $\hat{\Sigma}$ is an $n \times n$ diagonal matrix with singular values of $A$ on its diagonal, and $V$ is an $n \times n$ orthogonal matrix. So that

$$
A^{\mathrm{T}}A\vec{x} = V\Sigma^{\mathrm{T}}\hat{U}^{\mathrm{T}}\hat{U} \hat{\Sigma} V^{\mathrm{T}}\vec{x} = V \Sigma\hat{U}^{\mathrm{T}} \vec{b} = A^{\mathrm{T}}\vec{b}
\\[0.8em]
\vec{x} = {V^{\mathrm{T}}}^{-1} \hat{\Sigma}^{-1} \hat{U}^{-1} \vec{b} = V \hat{\Sigma}^{-1} \hat{U}^{\mathrm{T}} \vec{b}
$$

1. Compute reduced SVD $A = \hat{U} \hat{\Sigma} V^{\mathrm{T}}$
2. $\vec{x} = V \hat{\Sigma}^{-1} \hat{U}^{\mathrm{T}} \vec{b}$

# Stability of solving least squares problems
Here using $QR$ factorization is more stable than directly solve the linear equation. Let's see from an example.

>**e.g.***
>
>Fix $e^{\sin(4t)}$ from $0$ to $1$; we'll discretize the interval by $100$ points and use a polynomial with degree $14$. The function now looks like
>
>$$\left[
\begin{array}{ccccc}
1 & x_1 & x_1^2 & \cdots & x_1^{14} \\
1 & x_2 & x_2^2 & \cdots & x_2^{14} \\
1 & x_3 & x_3^2 & \cdots & x_3^{14} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_{100} & x_{100}^2 & \cdots & x_{100}^{14} 
\end{array}
\right]\left[
\begin{array}{c}
a_1 \\
a_2 \\
a_3 \\
\vdots \\
a_{15} 
\end{array}
\right]
=c \cdot \left[
\begin{array}{c}
e^{\sin(4x_1)} \\
e^{\sin(4x_2)} \\
e^{\sin(4x_3)} \\
\vdots \\
e^{\sin(4x_{100})} 
\end{array}
\right]
$$

>Here, $c$ is used to scale $a_{15}$ to $1$, as a checking solution, $c = 1/2006.787678808116$. In the last case, where the normal equation is solved, $a_{15} = 0.24664$, completely different from $1$.

$Conclusion$



For [method 1](#Solve-normal-equation), $\Cond(A^{\mathrm{T}}A) = \big\|A^{\mathrm{T}}A\big\|_2\left\|\left(A^{\mathrm{T}}A\right)^{-1}\right\|_2$; then

$$\begin{align}
\big\|A^{\mathrm{T}}A\big\|_2 &= \sup_{\vec{x} \neq \vec{0}}\frac{\big\|A^{\mathrm{T}}A\vec{x}\big\|_2} {\left\|\vec{x}\right\|_2} = \lambda_{\text{max}}^2 \\
\big\|\left(A^{\mathrm{T}}A\right)^{-1}\big\|_2 &= \frac{1} {\lambda_{\text{min}}^2}\\
\Longrightarrow \Cond(A^{\mathrm{T}}A) &= \frac{\lambda_{\text{max}}^2} {\lambda_{\text{min}}^2}
\end{align}$$

So that actually $\Cond(A^{\mathrm{T}}A) = \left( \Cond(A) \right)^2$

For [method 2](#By-QR-factorization), $\Cond(A) = \Cond(\hat{R})$.

$$\begin{align}
\Cond(A) &= \Cond(\hat{Q}\hat{R}) \\
&= \left\| \hat{Q}\hat{R} \right\| \cdot \left\| \left( \hat{Q}\hat{R} \right)^{+} \right\| \\
&= \sup_{\vec{x} \neq \vec{0}} \frac{\left\| \hat{Q}\hat{R}\vec{x} \right\|} {\left\| \vec{x} \right\|} \cdot \sup_{\vec{x} \neq \vec{0}} \frac{\left\| \left( \hat{Q}\hat{R} \right)^{+}\vec{x} \right\|} {\left\| \vec{x} \right\|} \\
&= \left\| \hat{R} \right\| \cdot \left\| \left( \hat{R} \right)^{+} \right\| \\
&= \Cond(\hat{R})
\end{align}$$

$Theorem$ 3

Let the full rank least squares problem $\min_{\vec{x} \in \mathbb{R}^n}\|\vec{b} - A\vec{x}\|_2^2$ be solved by $QR$ factorization with Householder triangularization. This algorithm is stable in the sense that the computed solution $\tilde{\vec{x}}$ is the minimizer of a nearby least squares problem, $i.e.$, $\tilde{\vec{x}}$ minimizes $\left\|\vec{b} - \left(A + \delta A \right)\vec{x}\right\|_2^2$ , where $\left\| \delta A \right\|_2 = \| A \|_2 O(\epsilon_{\text{machine}} )$.

$Theorem$ 4

The solution of the full rank least squares problem $\min_{\vec{x} \in \mathbb{R}^n}\|\vec{b} - A\vec{x}\|_2^2$ via solving the normal equation $A^{\mathrm{T}}A \vec{x} = A^{\mathrm{T}}\vec{b}$ is unstable.

cc8b0e284656b6ef3f22cfe921df3a27ec014c74
ipynb
Jupyter Notebook
Computational/Intro to Numerical Computing/Note_Chap07_Solving Linear Least Squares Problems.ipynb
XavierOwen/Notes
d262a9103b29ee043aa198b475654aabd7a2818d
MIT
2018-11-27T10:31:08.000Z
2019-01-20T03:11:58.000Z
Computational/Intro to Numerical Computing/Note_Chap07_Solving Linear Least Squares Problems.ipynb
XavierOwen/Notes
d262a9103b29ee043aa198b475654aabd7a2818d
MIT
Computational/Intro to Numerical Computing/Note_Chap07_Solving Linear Least Squares Problems.ipynb
XavierOwen/Notes
d262a9103b29ee043aa198b475654aabd7a2818d
MIT
2020-07-14T19:57:23.000Z
2020-07-14T19:57:23.000Z
Qwen/Qwen-72B
1. YES
2. YES