Singular value decomposition

Logistics

• General announcements

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

  • 7 lectures left!

• Assignment 5 and Midterm 2:

  • Grading starting, we’ll keep you posted

What’s on the agenda for today?

• Last time:

  • Singular value decomposition

• Today

  • Application to solving least squares

• Reading: lecture notes 13/14

Singular value decomposition

Important properties of the SVD

• The columns of $\mathbf{V}$ ${\{{\mathbf{v}}_i\}}_{i=1}^r$ are eigenvectors of the psd matrix ${\mathbf{A}}^{⊺}\mathbf{A}$. $\{{\sigma }_i:1\le i\le n\text{ and }{\sigma }_i\ne 0\}$ are the square roots of the non-zero eigenvalues of ${\mathbf{A}}^{⊺}\mathbf{A}$.

• The columns of $\mathbf{U}$ ${\{{\mathbf{u}}_i\}}_{i=1}^r$ are eigenvectors of the psd matrix $\mathbf{A}{\mathbf{A}}^{⊺}$. $\{{\sigma }_i:1\le i\le n\text{ and }{\sigma }_i\ne 0\}$ are the square roots of the non-zero eigenvalues of $\mathbf{A}{\mathbf{A}}^{⊺}$.

• The columns of $\mathbf{V}$ form an orthobasis for $\text{row}(\mathbf{A})$

• The columns of $\mathbf{U}$ form an orthobasis for $\text{col}(\mathbf{A})$

• Equivalent form of the SVD: $\mathbf{A}=\tilde{\mathbf{U}}\tilde{\mathbf{\Sigma }}{\tilde{\mathbf{V}}}^T$ where

  • $\tilde{\mathbf{U}}\in {\mathbb{R}}^{m\times m}$ is orthonormal
  • $\tilde{\mathbf{V}}\in {\mathbb{R}}^{n\times n}$ is orthonormal
  • $\tilde{\mathbf{\Sigma }}\in {\mathbb{R}}^{m\times n}$ is

  $$\tilde{\mathbf{\Sigma }}\triangleq \left[\begin{array}{cc}\mathbf{\Sigma }& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right]$$

SVD and least-squares

• When we cannot solve $\mathbf{y}=\mathbf{A}\mathbf{x}$, we solve instead $$\underset{\mathbf{x}\in {\mathbb{R}}^n}{min}{\Vert \mathbf{x}\Vert }_2^2\text{ such that }{\mathbf{A}}^{⊺}\mathbf{A}\mathbf{x}={\mathbf{A}}^{⊺}\mathbf{y}$$

  • This allows us to pick the minimum norm solution among potentially infinitely many solutions of the normal equations.

• Recall: when $\mathbf{A}\in {\mathbb{R}}^{m\times n}$ is of rank $n$, then $\mathbf{x}={\mathbf{A}}^{⊺}(\mathbf{A}{\mathbf{A}}^{⊺}{)}^{-1}\mathbf{y}$

• The solution of $$\underset{\mathbf{x}\in {\mathbb{R}}^n}{min}{\Vert \mathbf{x}\Vert }_2^2\text{ such that }{\mathbf{A}}^{⊺}\mathbf{A}\mathbf{x}={\mathbf{A}}^{⊺}\mathbf{y}$$ is $$\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$$ where $\mathbf{A}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T$ is the SVD of $\mathbf{A}$.

Pseudo inverse

• ${\mathbf{A}}^{+}=\mathbf{V}{\mathbf{\Sigma }}^{-1}{\mathbf{U}}^{⊺}$ is called the pseudo-inverse, Lanczos inverse, or Moore-Penrose inverse of $\mathbf{A}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T$.

• If $\mathbf{A}$ is square invertible then ${\mathbf{A}}^{+}=\mathbf{A}$

• If $m\ge n$ (tall and skinny matrix) of rank $n$ then ${\mathbf{A}}^{+}=({\mathbf{A}}^{⊺}\mathbf{A}{)}^{-1}{\mathbf{A}}^{⊺}$

• If $m\ge m$ (short and fat matrix) of rank $m$ then ${\mathbf{A}}^{+}={\mathbf{A}}^{⊺}(\mathbf{A}{\mathbf{A}}^{⊺}{)}^{-1}$

• Note ${\mathbf{A}}^{+}$ is as “close” to an inverse of $\mathbf{A}$ as possible

Stability of least squares

• 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$$

Stable reconstruction by truncation

• 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\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$$

Stable reconstruction by regularization

• 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}$$

Numerical methods

• 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

Easy systems

• Diagonal system
  • $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ invertible and diagonal
  • $O(n)$ complexity
• Orthogonal system
  • $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ invertible and orthogonal
  • $O(n^2)$ complexity
• Lower triangular system
  • $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ invertible and lower diagonal
  • $O(n^2)$ complexity
• General strategy: factorize $\mathbf{A}$ to recover some of the structures above