Stability and Numerical Aspects of Least Squares

General announcements

• Assignment 6 to be posted… (grading traffic jam)

• 4 lectures left! (No lecture on Wednesday November 24, 2021)

Midterm 2

• 99% graded, grades released after Thanksiving weekend
• Midterm solution during office hours on Tuesday November 23, 2021

Assignment 5

• Grading finalized

Assignment 6 and 7

• Posted this week, due date on Monday December 6 and Monday December 13

Last time:

• Stability of least squares

Today:

• More on stability and numerical considerations

Reading: lecture notes 14/15/16

What if we observe $\mathbf{y}=\mathbf{A}{\mathbf{x}}_0+\mathbf{e}$ and we apply the pseudo inverse? $\hat{\mathbf{x}}={\mathbf{A}}^{+}\mathbf{y}$

We can separate the error analysis into two components $$\hat{\mathbf{x}}-{\mathbf{x}}_0=\underset{\text{null space error}}{\underbrace{{\mathbf{A}}^{+}\mathbf{A}{\mathbf{x}}_0-{\mathbf{x}}_0}}+\underset{\text{noise error}}{\underbrace{{\mathbf{A}}^{+}\mathbf{e}}}$$

We will express the error in terms of the SVD $\mathbf{A}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{⊺}$ With

• ${\{{\mathbf{v}}_i\}}_{i=1}^r$ orthobasis of $\text{row}(\mathbf{A})$, augmented by ${\{{\mathbf{v}}_i\}}_{i=1}^{r+1}\in \ker \mathbf{A}$ to form an orthobasis of ${\mathbb{R}}^n$
• ${\{{\mathbf{u}}_i\}}_{i=1}^r$ orthobasis of $\text{col}(\mathbf{A})$, augmented by ${\{\mathbf{u}\}}_{i=1}^{r+1}\in \ker {\mathbf{A}}^{⊺}$ to form an orthobasis of ${\mathbb{R}}^m$

The null space error is given by $${\Vert {\mathbf{A}}^{+}\mathbf{A}{\mathbf{x}}_0-{\mathbf{x}}_0\Vert }_2^2=\sum _{i=r+1}^n{\left|{⟨{\mathbf{v}}_i,x_0⟩}_{}\right|}^2$$

The noise error is given by $${\Vert {\mathbf{A}}^{+}\mathbf{e}\Vert }_2^2=\sum _{i=1}^r\frac{1}{{\sigma }_i^2}{\left|{⟨\mathbf{e},{\mathbf{u}}_i⟩}_{}\right|}^2$$

How do we mitigate the effect of small singular values in reconstruction? $$\hat{\mathbf{x}}=\mathbf{V}{\mathbf{\Sigma }}^{-1}{\mathbf{U}}^{⊺}\mathbf{y}=\sum _{i=1}^r\frac{1}{{\sigma }_i}{⟨\mathbf{y},{\mathbf{u}}_i⟩}_{}{\mathbf{v}}_i$$

Truncate the SVD to $r^{\prime }<r$ $${\mathbf{A}}_t\triangleq \sum _{i=1}^{r^{\prime }}{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^{⊺}{\mathbf{A}}_t^{+}=\sum _{i=1}^{r^{\prime }}\frac{1}{{\sigma }_i}{\mathbf{u}}_i{\mathbf{v}}_i^{⊺}$$

Reconstruct $\widehat{{\mathbf{x}}_t}=\sum _{i=1}^{r^{\prime }}\frac{1}{{\sigma }_i}{⟨\mathbf{y},{\mathbf{u}}_i⟩}_{}{\mathbf{v}}_i={\mathbf{A}}_t$

Error analysis: $${\Vert {\hat{\mathbf{x}}}_t-{\mathbf{x}}_0\Vert }_2^2=\sum _{i=r+1}^n{\left|{⟨{\mathbf{x}}_0,{\mathbf{v}}_i⟩}_{}\right|}^2+\sum _{i=r^{\prime }+1}^r{\left|{⟨{\mathbf{x}}_0,{\mathbf{v}}_i⟩}_{}\right|}^2+{\sum _{i=1}^r}^{\prime }\frac{1}{{\sigma }_i^2}{\left|{⟨\mathbf{e},{\mathbf{u}}_i⟩}_{}\right|}^2$$

Regularization means changing the problem to solve $$\underset{\mathbf{x}\in {\mathbb{R}}^n}{min}{\Vert \mathbf{y}-\mathbf{A}\mathbf{x}\Vert }_2^2+\lambda {\Vert \mathbf{x}\Vert }_2^2\lambda >0$$

The solution is $$\hat{\mathbf{x}}=({\mathbf{A}}^{⊺}\mathbf{A}+\lambda \mathbf{I}{)}^{-1}{\mathbf{A}}^{⊺}\mathbf{y}=\mathbf{V}({\mathbf{\Sigma }}^2+\lambda \mathbf{I}{)}^{-1}\mathbf{\Sigma }{\mathbf{U}}^{⊺}\mathbf{y}$$

We have seen several solutions to systems of linear equations $\mathbf{A}\mathbf{x}=\mathbf{y}$ so far

• $\mathbf{A}$ full column rank: $\hat{\mathbf{x}}=({\mathbf{A}}^{⊺}\mathbf{A}{)}^{-1}{\mathbf{A}}^{⊺}\mathbf{y}$
• $\mathbf{A}$ full row rank: $\hat{\mathbf{x}}={\mathbf{A}}^{⊺}(\mathbf{A}{\mathbf{A}}^{⊺}{)}^{-1}\mathbf{y}$
• Ridge regression: $\hat{\mathbf{x}}=({\mathbf{A}}^{⊺}\mathbf{A}+\delta \mathbf{I}{)}^{-1}{\mathbf{A}}^{⊺}\mathbf{y}$
• Kernel regression: $\hat{\mathbf{x}}=(\mathbf{K}+\delta \mathbf{I}{)}^{-1}\mathbf{y}$
• Ridge regression in Hilbert space: $\hat{\mathbf{x}}=({\mathbf{A}}^{⊺}\mathbf{A}+\delta \mathbf{G}{)}^{-1}{\mathbf{A}}^{⊺}\mathbf{y}$

Extension: constrained least-squares $$\underset{\mathbf{x}\in {\mathbb{R}}^n}{min}{\Vert \mathbf{y}-\mathbf{A}\mathbf{x}\Vert }_2^2\text{ s.t. }\mathbf{x}=\mathbf{B}\mathit{\alpha }\text{ for some }\mathit{\alpha }$$

• The solution is $\hat{\mathbf{x}}=\mathbf{B}({\mathbf{B}}^{⊺}{\mathbf{A}}^{⊺}\mathbf{A}\mathbf{B}{)}^{-1}{\mathbf{B}}^{⊺}{\mathbf{A}}^{⊺}\mathbf{y}$

All these problems involve a symmetric positive definite system of equations.

• Many methods to achieve this based on matrix factorization

• $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ invertible and diagonal
• $O(n)$ complexity

• $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ invertible and orthogonal
• $O(n^2)$ complexity

• $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ invertible and lower diagonal
• $O(n^2)$ complexity

LU factorization

Cholesky factorization

QR decomposition

SVD and eigenvalue decompositions