MATH347DS L07: Least squares

Overview

• Gram-Schmidt algorithm, $𝑸𝑹$ factorization

• Projection onto subspaces

• Orthogonal projectors

• Best approximation in the 2-norm

• Linear regression

• Polynomial approximation

• Polynomial interpolation

Orthonormal vector sets

Definition. The Dirac delta symbol ${\delta }_{ij}$ is defined as

Definition. A set of vectors $\{{𝒖}_1,\dots ,{𝒖}_n\}$ is said to be orthonormal if

• The column vectors of the identity matrix are orthonormal

Gram-Schmidt algorithm

• An arbitrary vector set can be transformed into an orthonormal set by the Gram-Schmidt algorithm

• Idea:

  • Start with an arbitrary direction ${𝒂}_1$

  • Divide by its norm to obtain a unit-norm vector ${𝒒}_1={𝒂}_1/||{𝒂}_1||$

  • Choose another direction ${𝒂}_2$

  • Subtract off its component along previous direction(s) ${𝒂}_2-\left({𝒒}_1^T{𝒂}_2\right){𝒒}_1$

  • Divide by norm ${𝒒}_2=({𝒂}_2-\left({𝒒}_1^T{𝒂}_2\right){𝒒}_1)/||{𝒂}_2-\left({𝒒}_1^T{𝒂}_2\right){𝒒}_1||$

  • Repeat the above

Matrix formulation of Gram-Schmidt ($𝑸𝑹$ factorization)

• Consider $𝑨\in {ℝ}^{m\times n}$ with linearly independent columns. By linear combinations of the columns of $𝑨$ a set of orthonormal vectors ${𝒒}_1,\dots ,{𝒒}_n$ will be obtained. This can be expressed as a matrix product

  ${\displaystyle 𝑨=\left(\begin{array}{cccc}{𝒂}_1& {𝒂}_2& \dots & {𝒂}_n\end{array}\right)=\left(\begin{array}{cccc}{𝒒}_1& {𝒒}_2& \dots & {𝒒}_n\end{array}\right)\left(\begin{array}{ccccc}{r}_{11}& r_{12}& r_{13}& \dots & r_{1n}\\ 0& r_{22}& r_{23}& \dots & r_{2n}\\ 0& 0& r_{33}& \dots & r_{3n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \dots & r_{mn}\end{array}\right)=𝑸𝑹}$

  with $𝑸\in {ℝ}^{m\times n}$, $𝑹\in {ℝ}^{n\times n}$. The matrix $𝑹$ is upper-triangular (also referred to as right-triangular) since to find vector ${𝒒}_1$ only vector ${𝒂}_1$ is used, to find vector ${𝒒}_2$ only vectors ${𝒂}_1,{𝒂}_2$ are used

• The above is equivalent to the system

  ${\displaystyle \{\begin{array}{l}{𝒂}_1=r_{11}{𝒒}_1\\ {𝒂}_2=r_{12}{𝒒}_1+r_{22}{𝒒}_2\\ ⋮\\ {𝒂}_n=r_{1n}{𝒒}_1+r_{2n}{𝒒}_2+\dots +r_{nn}{𝒒}_n\end{array}}$

Matrix formulation of Gram-Schmidt ($𝑸𝑹$ factorization)

• The system can be solved to find $𝑹,𝑸$ by:

  1. Imposing $||{𝒒}_1||=1\Rightarrow r_{11}=||{𝒂}_1||$, ${𝒒}_1=a_1/r_{11}$

  2. Computing projections of ${𝒂}_2,\dots ,{𝒂}_n$ along ${𝒒}_1$

     ${\displaystyle r_{12}={𝒒}_1^T{𝒂}_2,\dots ,r_{1n}={𝒒}_1^T{𝒂}_n}$

  3. Subtracting components along ${𝒒}_1$ from ${𝒂}_2,\dots ,{𝒂}_n$

     ${\displaystyle \{\begin{array}{l}{𝒂}_2-r_{12}{𝒒}_1=r_{22}{𝒒}_2\\ ⋮\\ {𝒂}_n-r_{1n}{𝒒}_1=r_{2n}{𝒒}_2+\dots +r_{nn}{𝒒}_n\end{array}}$

  4. The above steps reduced the size of the system by 1. Repeating the steps completes the solution. The overall process is known as the Gram-Schmidt algorithm

$𝑸𝑹$ factorization

• For $𝑨\in {ℝ}^{m\times n}$ with linearly independent columns, the Gram-Schmidt algorithm furnishes a factorization

  ${\displaystyle 𝑸𝑹=𝑨}$

  with $𝑸\in {ℝ}^{m\times n}$ with orthonormal columns and $𝑹\in {ℝ}^{n\times n}$ an upper triangular matrix.

• Since the column vectors within $𝑸$ were obtained through linear combinations of the column vectors of $𝑨$ we have

  ${\displaystyle C\left(𝑨\right)=C\left(𝑸\right)}$

Orthogonal projection of a vector along another vector

• Consider a vector $𝒖\in {ℝ}^m$, and a unit-norm vector ${𝒒}_1\in {ℝ}^m$

Definition. The orthogonal projection of $𝒖\in {ℝ}^m$ along direction ${𝒒}_1\in {ℝ}^m$, $||{𝒒}_1||=1$ is the vector $\left({𝒒}_1^T𝒖\right){𝒒}_1$.

• Scalar-vector multiplication commutativity: $\left({𝒒}_1^T𝒖\right){𝒒}_1={𝒒}_1\left({𝒒}_1^T𝒖\right)$

• Matrix multiplication associativity: ${𝒒}_1\left({𝒒}_1^T𝒖\right)=\left({𝒒}_1{𝒒}_1^T\right)𝒖={𝑷}_1𝒖$, with ${𝑷}_1\in {ℝ}^{m\times m}$

Definition. The matrix ${𝑷}_1={𝒒}_1{𝒒}_1^T\in {ℝ}^{m\times m}$ is the orthogonal projector along direction ${𝒒}_1\in {ℝ}^m$, $||{𝒒}_1||=1$.

Orthogonal projection onto a subspace

• Consider $n$ orthonormal vectors grouped into a matrix $𝑸=\left(\begin{array}{lll}{𝒒}_1& \dots & {𝒒}_n\end{array}\right)\in {ℝ}^{m\times n}$

• The orthogonal projection of $𝒖$ onto the subspace spanned by ${𝒒}_1,\dots ,{𝒒}_n$ is

  ${\displaystyle 𝑷𝒖={𝑷}_1𝒖+\dots +{𝑷}_n𝒖=\left({𝒒}_1{𝒒}_1^T\right)𝒖+\dots {+}_{}\left({𝒒}_n{𝒒}_n^T\right)𝒖\Rightarrow }$ ${\displaystyle 𝑷={𝒒}_1{𝒒}_1^T+\dots +{𝒒}_n{𝒒}_n^T=\left(\begin{array}{ccc}{𝒒}_1& \dots & {𝒒}_n\end{array}\right)\left(\begin{array}{c}{𝒒}_1^T\\ ⋮\\ {𝒒}_n^T\end{array}\right)=𝑸{𝑸}^T}$

  Definition. The orthogonal projector onto $C\left(𝑸\right)$, $𝑸\in {ℝ}^{m\times n}$ with orthonormal column vectors is $𝑷=𝑸{𝑸}^T$

Complementary orthogonal projector

• Given $𝒖\in {ℝ}^m$ and $𝑸=\left(\begin{array}{lll}{𝒒}_1& \dots & {𝒒}_n\end{array}\right)\in {ℝ}^{m\times n}$ with orthonormal columns

Definition. The complementary orthogonal projector to $𝑷=𝑸{𝑸}^T$ is $𝑰-𝑷$, where $𝑸\in {ℝ}^{m\times n}$ is a matrix with orthonormal columns.

• The complementary orthogonal projector projects a vector onto the left null space, $N\left({𝑸}^T\right)$

Orthogonal projectors and linear systems

• Consider the linear system $𝑨𝒙=𝒃$ with $𝑨\in {ℝ}^{m\times n}$, $𝒙\in {ℝ}^n$, $𝒃\in {ℝ}^m$. Orthogonal projectors and knowledge of the four fundamental matrix subspaces allows us to succintly express whether there exist no solutions, a single solution of an infinite number of solutions:

  • Consider the factorization $𝑸𝑹=𝑨$, the orthogonal projector $𝑷=𝑸{𝑸}^T$, and the complementary orthogonal projector $𝑰-𝑷$

  • If $||(𝑰-𝑷)𝒃||\ne 0$, then $𝒃$ has a component outside the column space of $𝑨$, and $𝑨𝒙=𝒃$ has no solution

  • If $||(𝑰-𝑷)𝒃||=0$, then $𝒃\in C\left(𝑸\right)=C\left(𝑨\right)$ and the system has at least one solution

  • If $N\left(𝑨\right)=\left\{𝟎\right\}$ (null space only contains the zero vector, i.e., null space of dimension 0) the system has a unique solution

  • If $\dim N\left(𝑨\right)=n-r>0$, then a vector $𝒚\in N\left(𝑨\right)$ in the null space is written as

    ${\displaystyle 𝒚=c_1{𝒛}_1+\dots +c_{n-r}{𝒛}_{n-r}}$

    and if $𝒙$ is a solution of $𝑨𝒙=𝒃$, so is $𝒙+𝒚$, since

    ${\displaystyle 𝑨(𝒙+𝒚)=𝑨𝒙+c_1𝑨{𝒛}_1+\dots +c_{n-r}𝑨{𝒛}_{n-r}=𝒃+𝟎+\dots +𝟎=𝒃}$

    The linear system has an $(n-r)$-parameter family of solutions

Best approximation in the 2-norm (least squares)

Figure 1. Least squares problem: find $𝒗\in C\left(𝑨\right)$, $𝑨\in {ℝ}^{m\times n}$ closest to some given $𝒚$ in the 2-norm

• Mathematical statement: solve the minimization problem ${min}_{𝒄\in {ℝ}^n}||𝒚-𝑨𝒄||$

• Approach: project $𝒚$ onto the column space of $𝑨$:

  1. Find an orthonormal basis for column space of $𝑨$ by $𝑸𝑹$ factorization, $𝑸𝑹=𝑨$

  2. State that $𝒗$ is the projection of $𝒚$, $𝒗={𝑷}_{C\left(𝑨\right)}𝒚={𝑷}_{𝑸}𝒚=𝑸{𝑸}^T𝒚$

  3. State that $𝒗$ is within the column space of $𝑨$, $𝒗=𝑨𝒄=𝑸𝑹𝒄$

  4. Set equal the two expressions of $𝒗$, $$$𝑸{𝑸}^T𝒚=𝑸𝑹𝒄\Rightarrow 𝑹𝒄={𝑸}^T𝒚$

  5. Solve the triangular system to find $𝒄$ (in Julia, Matlab, Octave: c=R\(Q'y))

Least squares: linear regression calculus approach

• In many scientific fields the problem of determining the straight line $y\left(x\right)=c_0+c_{1}x$, that best approximate data $𝒟=\{(x_i,y_i),i=1,\dots ,m\}$ arises. The problem is to find the coefficients $c_0,c_1$, and this is referred to as the linear regression problem.

• The calculus approach: Form sum of squared differences between $y\left(x_i\right)$ and $y_i$

  ${\displaystyle S(c_0,c_1)=\sum _{i=1}^m{(y\left(x_i\right)-y_i)}^2=\sum _{i=1}^m{(c_0+c_1{x}_i-y_i)}^2}$

  and seek $(c_0,c_1)$ that minimize $S(c_0,c_1)$ by solving the equations

  ${\displaystyle \frac{\partial S}{\partial c_0}=0\Rightarrow 2\sum _{i=1}^m(c_0+c_1{x}_i-y_i)=0\iff m{c}_0+\left(\sum _{i=1}^m{x}_i\right)c_1=\sum _{i=1}^m{y}_i}$ ${\displaystyle \frac{\partial S}{\partial c_1}=0\Rightarrow 2\sum _{i=1}^m(c_0+c_1{x}_i-y_i)x_i=0\iff \left(\sum _{i=1}^m{x}_i\right)c_0+\left(\sum _{i=1}^m{x}_i^2\right)c_1=\sum _{i=1}^m{x}_i{y}_i}$

Geometry of linear regression

• Form a vector of errors with components $e_i=y\left(x_i\right)-x_i$. Recognize that $y\left(x_i\right)$ is a linear combination of $1$ and $x_i$ with coefficients $a_0,a_1$, or in vector form

  ${\displaystyle 𝒆=\left(\begin{array}{cc}1& x_1\\ ⋮& ⋮\\ 1& x_m\end{array}\right)\left(\begin{array}{c}{c}_0\\ c_1\end{array}\right)-𝒚=\left(\begin{array}{cc}𝟏& 𝒙\end{array}\right)𝒄-𝒚=𝑨𝒄-𝒚}$

• The norm of the error vector $||𝒆||$ is smallest when $𝑨𝒄$ is as close as possible to $𝒚$. Since $𝑨𝒄$ is within the column space of $C\left(𝑨\right)$, $𝑨𝒄\in C\left(𝑨\right)$, the required condition is for $𝒆$ to be orthogonal to the column space

  ${\displaystyle 𝒆\perp C\left(𝑨\right)\Rightarrow {𝑨}^T𝒆=\left(\begin{array}{c}{𝟏}^T\\ {𝒙}^T\end{array}\right)𝒆=\left(\begin{array}{c}{𝟏}^T𝒆\\ {𝒙}^T𝒆\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)=𝟎}$ ${\displaystyle {𝑨}^T𝒆=𝟎\iff {𝑨}^T(𝑨𝒄-𝒚)=0\iff \left({𝑨}^T𝑨\right)𝒄={𝑨}^T𝒚}$