Posted: 09/09/21

Due: 09/15/21, 11:55PM

This assigment addresses one of the fundamental topics within scientific computation: finding economical descriptions of complex objects. Some object is described by $𝒚\in {ℂ}^m$ (with $m$ typically large), and a reduced description is sought by linear combination $𝑨𝒙$, with $𝑨\in {ℂ}^{m\times n}$ ($n<m$, often $n\ll m$). The surprisingly simple Euclidean geometry of Fig. 1 (which should be committed to memory) will be shown to have wide-ranging applicability to many different types of problems. The error (or residual) in approximating $𝒚$ by $𝑨𝒙$ is defined as

and 2-norm minimization defines the least-squares problem

Figure 1. Least squares (2-norm error minimization) problem.

1Track 1

Consider data $𝒟=\left\{(t_i,y_i)|i=1,2,\dots ,m.\right\}$ obtained by sampling a function $f:ℝ\to ℝ$, with $y_i=f\left(t_i\right)$. An approximation is sought by linear combination

Introduce the vector-valued function $A:ℝ\to {ℝ}^n$ (organized as a row vector)

such that

With $𝒕={\left[\begin{array}{llll}{t}_1& t_2& \dots & t_m\end{array}\right]}^T$ a sampling of the function domain, a matrix is defined by

Tasks. In each exercise below, construct the least-squares approximant for the stated range of $n\in 𝒩$, sample points $𝒕$, and choice of $A\left(t\right)$. Plot in a single figure all components of $A\left(t\right)$. Plot the approximants, as well as $f$ in a single figure. Construct a convergence plot of the approximations by representation of point data $ℰ=\left\{(\log n,\log ||𝒚-𝑨𝒙||)|𝑨\in {ℝ}^{m\times n},n\in 𝒩.\right\}$. For the largest value of $n$ within $𝒩$, construct a figure superimposing increasing number of sampling points, $m\in ℳ$. Comment on what you observe in each individual exercise. Also compare results from the different exercises.

1. Start with the classical example due to Runge (1901)

   ${\displaystyle \begin{array}{l}f:[-1,1]\to ℝ,f\left(t\right)=\frac{1}{(1+25t^2)},t_i=\frac{2(i-1)}{m-1}-1,\\ ℳ=\{16,32,64,128,256\},𝒩=\{4,8,16,32\},\\ A\left(t\right)=\left[\begin{array}{lllll}1& t& t^2& \dots & t^{n-1}\end{array}\right].\end{array}}$

   Solution. With $𝒕\in {ℝ}^m$ denoting the sampling point vector, and $𝒚\in {ℝ}^m$, the function values at the sample points, the least squares problem is

   ${\displaystyle {min}_{𝒙\in {ℝ}^n}||𝒚-𝑨𝒙||,}$

   where

   ${\displaystyle 𝑨=\left[\begin{array}{llll}𝟏& 𝒕& \dots & {𝒕}^{n-1}\end{array}\right].}$

   The solution to the least squares problem

   ${\displaystyle 𝒛={\mathrm{argmin}}_{𝒙\in {ℝ}^n}||𝒚-𝑨𝒙||,}$

   furnishes the approximation

   ${\displaystyle f\left(t\right)\cong \stackrel{\sim }{f}\left(t\right)=A\left(t\right)𝒛,}$

   that can be sampled at $M⩾m$ points to assess approximation error.

   When carrying convergence studies such as these, it is convenient to define functions for common tasks:

   $\circ$sample. Returns a sample of $f:[a,b]\to ℝ$ at $m$ equidistant points

   ∴	function sample(a,b,f,m) t = LinRange(a,b,m); y = f.(t) return t,y end;

   ∴

   $\circ$plotLSQ. Constructs a figure with plots of:

   1. The $m$ sample points (i.e., data) represented as black dots;

   2. The function sampled at more points, i.e., $M⩾m$;

   3. The approximation sampled at $M$ points

   ∴	function plotLSQ(a,b,f,Basis,m,n,M) data=sample(a,b,f,m); t=data[1]; y=data[2] Data=sample(a,b,f,M); T=Data[1]; Y=Data[2] A = Basis(t,n); x = A\y; z = Basis(T,n)*x plot(t,y,"ok",T,z,"-r",T,Y,"-b"); grid("on"); xlabel("t"); ylabel("y"); title("Least squares approximation") end;

   ∴

   This problem solution is obtained by:

   1. defining the function $f$

      ∴	Runge(t)=1/(1+25*t^2);

      ∴

   2. defining the basis

      ∴	function MonomialBasis(t,n) m=size(t)[1]; A=ones(m,1); for j=1:n-1 A = [A t.^j] end return A end;

      ∴

   3. Invoking plotLSQ with appropriate parameters

      ∴	clf(); plotLSQ(-1,1,Runge,MonomialBasis,16,4,64);

      ∴	FigPrefix=homedir()*"/courses/MATH661/images/H03";

      ∴	savefig(FigPrefix*"Fig01.eps")

      ∴

      Figure 2. Least squares approximant (red) of Runge function (blue) sampled at (black dots).

   Once the above are defined cycling through the parameter ranges is straightforward (open figure folds to see code).

   $\circ$

   Figure 3. Effect of increasing number of monomial basis functions in least squares approximation of Runge function.

   ∴	for p=2:5 clf(); n=2^p for q=4:8 m=2^q; plotLSQ(-1,1,Runge,MonomialBasis,m,n,1024) end savefig(FigPrefix*"Fig01n="*string(n)*".eps") end

   ∴

2. Instead of the equidistant point samples of the Runge example above use the Chebyshev nodes

   ${\displaystyle t_i=\cos \left(\frac{2i-1}{2m}\pi \right),}$

   keeping other parameters as in Problem 1.

3. Instead of the monomial family of the Runge example, use the Fourier basis

   ${\displaystyle A\left(t\right)=\left[\begin{array}{llllll}1& \cos \pi t& \sin \pi t& \dots & \cos \pi nt& \sin \pi nt\end{array}\right]}$

   keeping other parameters as in Problem 1. In this case $𝑨\in {ℝ}^{m\times (2n+1)}$.

4. Instead of the monomial family of the Runge example, use the piecewise linear $B$-spline basis

   ${\displaystyle A\left(t\right)=\left[\begin{array}{llll}{N}_1\left(t\right)& N_2\left(t\right)& \dots & N_n\left(t\right)\end{array}\right],}$ ${\displaystyle N_i\left(t\right)=\{\begin{array}{ll}0,& t<t_{i-1}\\ \frac{t-t_{i-1}}{h}& t_{i-1}⩽t<t_i\\ \frac{t_{i+1}-t}{h}& t_i⩽t<t_{i+1}\\ 0& t_{i+1}<t\end{array}.,h=\frac{2}{m-1},}$

   keeping other parameters as in Problem 1.

2Track 2

1. If $𝑸\in {ℂ}^{m\times n}$ has orthonormal columns, prove that ${𝑷}_{𝑸}=𝑸{𝑸}^{\ast }$ is an orthogonal projector onto $C\left(𝑸\right)$. Determine the expression of ${𝑷}_{𝑨}$, the projector onto $C\left(𝑨\right)$, with $𝑨\in {ℂ}^{m\times n}$. Compare the number of arithmetic operations required to compute $𝒚={𝑷}_{𝑨}𝒙$, by comparison to first determining the $QR$ factorization, $𝑨=𝑸𝑹$, and then computing $𝒚=𝑸{𝑸}^{\ast }𝒙$.

2. Continuing Problem 1, determine ${||{𝑷}_{𝑸}||}_2$, and express ${||{𝑷}_{𝑨}||}_2$ in terms of the singular value decomposition of $𝑨$. Comment the result, considering, say, length of shadows at various times of day.

3. A matrix $𝑨=\left[a_{ij}\right]\in {ℂ}^{m\times n}$ is said to be banded with bandwidth $B$ if $a_{ij}=0$ for $|i-j|>B$. Implement the modified Gram-Schmidt algorithm for $𝑨\in {ℂ}^{m\times n}$ a banded matrix with bandwidth $B$ using as few arithmetic operations as possible.

4. Solve Problem 1, Track 1.

5. Solve Problem 4, Track 1.

6. In Problem 1, Track 1, replace the monomial basis with the Legendre polynomials, whose samples are determined by $QR$ decomposition $𝑸𝑹=𝑨$. The resulting least squares problem is now

   ${\displaystyle {min}_{𝒙\in {ℝ}^n}{||𝒚-𝑸𝒙||}_2.}$