% Hill4.tex

\section{Linear Transformations, Orthogonal Projections and Least Squares}

\subsection{Matrices as Linear Transformations}

You are expected to know the main facts about functions, and in particular
linear transformations.

Similarly the geometric applications, and interpretations, of matrices and
linear transformations given in Hill should be familiar from your first
course in linear algebra.  Hill emphasizes the example (in Computer Aided
Design) of rotations and reflections in $\DoubleR^3$. If you have graphics
capabilities on your computer terminal, you can see this for yourself using
Maple's {\tt plot3d}\index{plot} function:

\begin{verbatim}
# See Hill Section 4.1 Figure 4.2 p237.
with(linalg):
# Instead of using the race car given by Hill, we will use the following
> surface:= [x*sin(x)*cos(y), x*cos(x)*cos(y), x*sin(y) ]:
> domain:= x=0..Pi, y=0..Pi:
> plot3d(surface,domain):

> Q:= matrix([ [cos(Pi/4),sin(Pi/4),0], [-sin(Pi/4),cos(Pi/4),0], [0,0,1] ]):
> rotated:= convert( evalm(Q&*surface), list):   # To make it possible to plot
> plot3d(rotated,domain):

> Q:= matrix([ [-1,0,0], [0,-1,0], [0,0,1] ]):  # Rotation of Pi in the xy-plane
> rotated:= convert( evalm(Q&*surface), list):
> plot3d(rotated,domain):

> Q1:= matrix([ [1,0,0], [0,cos(Pi/3),sin(Pi/3)], [0,-sin(Pi/3),cos(Pi/3)] ]):
> Q2:= matrix([ [cos(Pi/3),0,sin(Pi/3)], [0,1,0], [-sin(Pi/3),0,cos(Pi/3)] ]):
> rotated:= convert(evalm(Q1&*Q2&*surface),list): #Rotations in xz and yz-planes
> plot3d(rotated,domain):

> Q:= matrix([ [1,0,0], [0,-1,0], [0,0,1] ]):   # Reflection through xz-plane
> reflected:= convert( evalm(Q&*surface), list):
> plot3d(reflected,domain):
\end{verbatim}

You are also expected to know the fact that a function
$T:\DoubleR^n\to\DoubleR^m$ is a linear transformation if and only if it is the
induced function of an $m\times n$ matrix.  This is illustrated below,
to-gether with a proof that differentiation is a linear transformation.
\begin{verbatim}
> with(linalg):
# Hill Section 4.1 Theorem 4.3 p241.
# Show that the induced function of  A  is a linear transformation where
> A:= matrix( [[a11,a12,a13],  
>              [a21,a22,a23]] ):
> u:= vector( [u1,u2,u3] ):                         v:= vector( [v1,v2,v3] ):
> evalm( A&*(u+v) - (A&*u+A&*v) ); 
                                    [ 0, 0 ]

> evalm( A&*(r*u) - r*A&*u ); 
                                    [ 0, 0 ]
\end{verbatim}	\label{Linear}
\begin{verbatim}
# Hill Section 4.1 Theorem 4.4 p242.
# Consider a linear transformation  T  from  R^3  to  R^2
# Ensure that  T  is a linear transformation:
> define( Linear(T) );            # This is the Maple V Release 3 version
> define( linear(T) );            # This is the version prior to Release 3
# Write
> T([1,0,0]):= vector( [a11,a21] ):         T([0,1,0]):= vector( [a12,a22] ):
> T([0,0,1]):= vector( [a13,a23] ):
# Show that the matrix A above is the matrix induced by  T.
> evalm( T(Pi*[1,0,0]+E*[0,1,0]+gamma*[0,0,1]) - A&*vector([Pi,E,gamma]) );
                                    [ 0, 0 ]

# Pi, E and gamma must be used since define/Linear only works with constants
#
# Hill Section 4.1 Example 8 p244.
# Show that differentiation is a linear transformation
> diff(f(x),x) + diff(g(x),x) - diff(f(x)+g(x),x);
                                       0

> diff(r*f(x),x) - r*diff(f(x),x);
                                       0
\end{verbatim}

\medskip
\headq{Note:} It is possible (and probably desirable) to look at Section 4.7 of
Hill at this point.

\subsection{Relationships involving Inner Products}
We again assume that, from an earlier course in linear algebra,
you know the properties of inner products in $\DoubleR^n$.

\smallskip
\headq{Definition 4.8} An inner product space is a vector space equipped with
  an inner product (which satisfies the properties of the dot product given in
  Hill's Theorem 4.5).

\headq{Definition 4.9} If \b u is a vector in an inner product space,
the {\bf norm} of \b u is $\Vert\b u\Vert = \sqrt{\b u\cdot \b u}$.

Maple's {\tt linalg} package has functions {\tt dotprod}\index{dotprod} and
{\tt norm}\index{dotprod} which calculate dot products and find norms
respectively of vectors in $\DoubleR^n$
(get {\tt help} on these functions). Note that this {\tt norm} is
different from the function {\tt norm} not in the {\tt linalg} package.\quad
To find the norm of the vector \b v, use
\begin{verbatim}
  norm(v,2);    # This is the standard norm  -  there are others !
\end{verbatim}
Here is another example of an inner product.

\begin{verbatim}
# Hill Section 4.2 Theorem 4.7 p247.
# Show that the following is an inner product:
> dot:= proc(fx,gx,lims)        int( expand(fx*gx),x=lims )        end:
> combine( dot(f(x),h(x),a..b) +dot(g(x),h(x),a..b) -dot((f+g)(x),h(x),a..b) );
                                       0

> combine( dot(f(x),g(x),a..b) - dot(g(x),f(x),a..b) );
                                       0

> combine( dot(r*f(x),g(x),a..b) - r*dot(f(x),g(x),a..b) );
                                       0

# combine  was needed above to force the linearity of integration.
#
# An example of this inner product was given in Hill chapter 3:
> f:= x->x:                               # Hill Section 3.3 Example 10 (p157)
> dot(f(x),sin(x),0..2*Pi);
                                     - 2 Pi

# The norm of  sin,  and the projection of  f  onto  sin  are given by
> norm_sin:= sqrt( dot(sin(x),sin(x),0..2*Pi) );
                                             1/2
                               norm_sin := Pi

> f_proj_sin:= (dot(f(x),sin(x),0..2*Pi)/dot(sin(x),sin(x),0..2*Pi)) * sin;
                             f_proj_sin := - 2 sin
\end{verbatim}
This inner product is closely related to Fourier series (see Hill \S 6.1).

\setcounter{theorem}{9}
\begin{theorem}[Cauchy-Schwarz inequality] Suppose $\b u,\b v$ are vectors in an
  inner product space $V$. Then
    $$ \vert\b u\cdot \b v\vert \leq \Vert\b u\Vert \times \Vert\b v\Vert. $$
\end{theorem}

\setcounter{theorem}{13}
\begin{theorem}[Triangle inequality] If $\b u,\b v$ are vectors in an
  inner product space $V$, then
    $$ \Vert\b u + \b v\Vert  \leq \Vert\b u\Vert + \Vert\b v\Vert. $$
\end{theorem}

\medskip\noindent
The following theorems give a convenient representation of dot products in
$\DoubleR^n$.

\begin{theorem} Let $A$ be an $m\times n$ matrix, $B$ an $n\times m$ matrix and
  \b x, \b y (column) vectors in $\DoubleR^n$, $\DoubleR^m$ respectively. Then
    $$ A\b x\cdot \b y = \b x\cdot A^T \b y\qquad{\rm and}\qquad
                                    \b x\cdot B\b y = \b B^T \b x\cdot \b y $$
\end{theorem}
\begin{theorem} Suppose we identify a $1\times 1$ matrix with its (single)
   entry.  If $\b u,\b v\in \DoubleR^n$ are represented as column vectors, then
    $$ \b u\cdot \b v = \b u^T \b v $$
\end{theorem}

\begin{verbatim}
# Hill Section 4.2 Theorems 4.15 and 4.16 pp 250,251.
> with(linalg):
> u:= vector( [u1,u2] ):                             v:= vector( [v1,v2] ):
> dotprod(u,v) - evalm( transpose(u)&*v );
                                       0

# Note that Maple treats vectors as COLUMN vectors (i.e. m by 1 matrices). So
> evalm( u&*transpose(v) );
                                [ u1 v1  u1 v2 ]
                                [              ]
                                [ u2 v1  u2 v2 ]

> A:= matrix(2,2):
> dotprod( evalm(A&*u),v ) - dotprod( u,evalm(transpose(A)&*v) ):
> simplify(");
                                       0
\end{verbatim}

\noindent The following theorem will be found to be useful in the next section.
\setcounter{theorem}{17}
\begin{theorem} Let $A$ be an $m\times n$ matrix. Then
  $$ NS(A^TA) = NS(A)\qquad{\rm and~so}\qquad rank(A^TA) = rank(A).$$
\end{theorem}
\headq{COROLLARY 4.19} {\em If the matrix $A$ has linearly independent columns,
  then $A^T A$ is invertible.}

\smallskip
\headq{EXERCISE} Show that, $A^T A$ is symmetric and, if the matrix $A$ has
   linearly independent columns,\quad $\b x\cdot\b y = \b x^TA^T A\b y$\quad
   is an inner product.

\medskip
\subsection{Least Squares and Orthogonal Projections}
Theorem 3.32 and the definition of a column space show that
\setcounter{theorem}{20}
\begin{theorem} If $A$ is a matrix, then the system\quad $A\b x=\b y$\quad is
  consistent if and only if \b y is in the column space $CS(A)$ of $A$.
\end{theorem}

In this section, Hill considers the situation where\quad $A\b x=\b y$\quad is
inconsistent (i.e. has no solution). In this case, it is natural to
look for an \b x which makes $A\b x$ as close as possible to \b y.

\smallskip
\headq{Definition 4.33} The {\bf least squares problem} for the inconsistent
  system\quad $A\b x=\b y$\quad is to find a vector $\bar\b x$
  (called the {\bf least squares solution} to the system) such that
  $\Vert\b y - A\bar \b x\Vert$ is a minimum.

\smallskip
The simplest, and possibly the most widely used, application is the fitting of
a straight line to data points as shown in Hill's figures 4.1 and 4.10.
Maple has a function {\tt leastsqrs}\index{leastsqrs} in the {\tt linalg}
package (Get {\tt help} on it). If $A$ is a matrix and \b y is a vector, then
{\tt leastsqrs(A,y)} finds the least squares solution of the system
$A\b x=\b y$. For example

\begin{verbatim}
# Hill Section 4.3 Example 6, p265
> with(linalg):
# Find the equation  y = m x + b  for the line of best fit for the points
#   (40,482), (45,467), (50,452), (55,433), (60,421)
> A:= matrix([ [40,1], [45,1], [50,1], [55,1], [60,1] ]):
> y:= vector( [482, 467, 452, 433, 421] ):
> xbar:= leastsqrs(A,y):                               map( evalf,xbar );
                             [ -3.120000000, 607. ]

# So the line of best fit is  y = -3.12 x + 607
#
# The same problem can be done with functions from the stats package.
# This has some advantages in connection with built-in graphics.
> with(stats):
#
# If you are using Maple V Release 3,
# remove the (first)  #  from the next 6 lines:
# dat:= [ [40,45,50,55,60], [482,467,452,433,421] ]:
# map( evalf, fit[leastsquare[[x,y],y=a+b*x, {a,b}]](dat) );
# # Other sorts of least-squares fits are possible.
# # Hill S4.3 Q24 p269 has an example of fitting a quadratic curve
# dat:= [ [1,2,3,4], [-1,0,-1,-1] ]:
# map( evalf, fit[leastsquare[[x,y],y=a0+a1*x+a2*x^2 ]](dat) );
#
# If you are using Maple V Release 1 or 2,
# remove the (first)  #  from the next 10 lines:
> dat:= array([ [y,x], [482,40], [467,45], [452,50], [433,55], [421,60] ]):
> v:= regression( dat,y=m*x+b );
                    v := {b = 607.0000000, m = -3.120000000}

> assign(v);              # This assigns the above values of b and m to them
> statplot( dat,y=m*x+b ):
# Other sorts of least-square fits are possible with the regression function.
# Hill S4.3 Q24 p269 has an example of fitting a quadratic curve
> dat:= array([ [y,x],[-1,1],[0,2],[-1,3],[-1,4] ]):
> v:= regression( dat, y = a0 + a1*x + a2*x^2 );
         v := {a1 = 1.150000000, a0 = -1.750000000, a2 = -.2500000000}

#  Run the above for answers, below for a plot.
> assign(v);                               statplot( dat, y = a0+a1*x+a2*x^2 ):
\end{verbatim}

\noindent Two approaches to the least-squares solution are given in the
following theorems.
\setcounter{theorem}{33}
\begin{theorem} $\Vert\b y - A\b x\Vert$ is a minimum if and only
  if $\b y - A\b x\perp CS(A)$.
\end{theorem}

\begin{theorem} Let $A$ be an $m\times n$ matrix and $\b y\in \DoubleR^m$.
  A least-squares solution to the inconsistent system $A\b x=\b y$
  is a solution to the {\bf normal} equation
    $$ A^T A {\bar\b x} = A^T \b y .$$
  If the columns of $A$ are linearly independent, then\quad
  ${\bar\b x} = (A^T A)^{-1} A^T \b y$\quad {\em (See Corollary 4.19).}
\end{theorem}

However, as mentioned in Section 1.6, the equation $A^T A\b x = A^T \b y$
should not be solved using the inverse.  In fact, even Gaussian elimination
presents problems here as {\em the matrix $A^T A$ is almost always
ill-conditioned}. In the next two sections (and also in Section 6.2),
Hill looks at the problem of calculating the least squares solution of a system.
This is done by a standard vector technique of projection.

\smallskip
\headq{Definition 4.38} Suppose $W$ is a subspace of $\DoubleR^m$ and
  $\b y\in \DoubleR^m$.
  Then the {\bf (perpendicular) projection} of $\b y$ into $W$ is the unique
  vector $\b p\in W$ such that $\b y - \b p\perp W$ and the function which
  maps each $\b y$ onto its projection into $W$ is called the
  {\bf (perpendicular) projection (function)} into $W$.

\smallskip
\headq{EXERCISE} Show that the perpendicular projection function is a linear
  transformation.

\setcounter{theorem}{38}
\begin{theorem} Assume $A$ has linearly independent columns.
  Then the matrix of the perpendicular projection function into $CS(A)$ is\quad
  $P=A (A^T A)^{-1} A^T$.
\end{theorem}

\headq{EXERCISE} Show that $P$ satisfies the following properties:

  \smallskip
  \qquad(a) $P^2 = P$.\qquad\qquad\qquad(b) $P^T = P$.

 \smallskip
 Conversely, any symmetric matrix with $P^2=P$ represents a projection function.

\eject
\begin{verbatim}
# Continuing Example 6 of Hill Section 4.3
> Atr:= transpose(A):                            p:= evalm( A&*xbar):
> equal( convert(evalm(Atr &* p),matrix), convert(evalm(Atr &* y),matrix) );
                                      true

# verifying Hill's  normal equations  4.36.
#
#  p  is the projection of  y  into  CS(A).
> dotprod( y-p, col(A,1) ),                     dotprod( y-p, col(A,2) );
                                      0, 0

# showing that  y - p  is orthogonal to both the columns of A - and so to CS(A)
#
# The matrix of the perpendicular projection function is
> P:= evalm( A &* inverse(Atr&*A) &* Atr );
\end{verbatim}
\begin{verbatim}
                           [  3/5   2/5  1/5    0   -1/5 ]
                           [                             ]
                           [  2/5  3/10  1/5  1/10    0  ]
                           [                             ]
                      P := [  1/5   1/5  1/5   1/5   1/5 ]
                           [                             ]
                           [   0   1/10  1/5  3/10   2/5 ]
                           [                             ]
                           [ -1/5    0   1/5   2/5   3/5 ]
\end{verbatim}

\subsection{Orthogonal Bases and the Gram-Schmidt Process}
This section shows how the idea of the perpendicular projection can be used to
find a special basis for a vector space.

\smallskip
\headq{Definitions 4.45 and  4.48} The vectors $\rowi{\b v}{m}$ are

  \quad {\bf orthogonal}\quad if $\b v_i\cdot \b v_j = 0$ whenever $i\ne j$.

  \quad {\bf orthonormal}\quad if they are orthogonal and each
  $\Vert\b v_i\Vert=1$.

\setcounter{theorem}{45}
\begin{theorem} Any set of nonzero orthogonal vectors is linearly independent.
\end{theorem}

In particular, an orthogonal spanning set is an orthogonal basis.

\setcounter{theorem}{48}
\begin{theorem} Let $B=\{\rowi{\b v}{n}\}$ be an ordered orthonormal basis for
  the vector space $V$ and $\b y\in V$. Then the coordinates of \b y relative
  to $B$ are $[\b y\cdot\b v_1, \b y\cdot\b v_2,\ldots \b y\cdot\b v_n]$.
\end{theorem}

\subsubsection*{The Gram-Schmidt Process}
This process starts with linearly independent vectors $\rowi{\b v}{n}$
and ends with orthonormal vectors $\rowi{\b u}{n}$.  At step $k$
it subtracts from $\b v_k$ its projection onto span$\{\rowi{\b v}{k-1}\}$:
  $$\b w_k = \b v_k  - \sum_{i=1}^{k-1} (\b u_i\cdot \b v_k) \b u_i ,$$
and then $\b u_k$ is defined to be the unit vector $\b w_k/\Vert\b w_k \Vert$.

\index{GramSchmidt}
\begin{verbatim}
# Hill Section 4.4 Example 7, p276
> with(linalg):
# Find an orthonormal basis for the space spanned by
> v1:= vector([1,1,1,1]):    v2:= vector([0,1,1,1]):    v3:= vector([0,0,1,1]):
# This can be done using the function  GramSchmidt  in the  linalg  package:
> GramSchmidt( [v1,v2,v3] );
        [[ 1, 1, 1, 1 ], [ -3/4, 1/4, 1/4, 1/4 ], [ 0, -2/3, 1/3, 1/3 ]]

# Warning: Although these vectors are orthogonal, they are not orthonormal,
# and so each should be divided by its norm to give an orthonormal basis
#
# To help understand the process, we also do it using projections:
> w1:= v1:                                    u1:= evalm( w1/norm(w1,2) ):
> w2:= evalm( v2 - dotprod(u1,v2)*u1 ):       u2:= evalm( w2/norm(w2,2) ):
> w3:= evalm( v3 - dotprod(u1,v3)*u1 - dotprod(u2,v3)*u2 ):
> u3:= evalm( w3/norm(w3,2) ):
> print (u1,u2,u3);
                                      1/2       1/2       1/2       1/2
     [ 1/2, 1/2, 1/2, 1/2 ], [ - 1/2 3   , 1/6 3   , 1/6 3   , 1/6 3    ],

                     1/2  1/2       1/2  1/2       1/2  1/2
         [ 0, - 1/3 2    3   , 1/6 2    3   , 1/6 2    3    ]

# Apart from the norms of the vectors, this agrees with the GramSchmidt result.
\end{verbatim}
As Hill points out (in the remark on p. 277), it is computationally easier
to simply find the vectors $\b w_1,\b w_2,\b w_3$, and then divide by their
norms. This is why Maple's function {\tt GramSchmidt} does so.

\subsection{Orthogonal matrices and QR Decomposition}

\headq{Definition 4.58} An {\bf orthogonal} matrix is
a square matrix with orthonormal columns.

\smallskip
\headq{Warning} A square matrix with orthogonal columns is not called
  orthogonal unless its columns also have norm 1.

\smallskip
For example, any permutation matrix is an orthogonal matrix.\quad
Another example is given below, using the function {\tt orthog}\index{orthog}
(in Maple's {\tt linalg} package) which returns the
value {\tt true} if the given matrix is orthogonal.
\begin{verbatim}
# Hill Section 4.5 Example 2 p282.
> with(linalg):
# Show that the matrix Q is orthogonal, where
> Q:= matrix([ [cos(theta),sin(theta)],
>              [-sin(theta),cos(theta)] ]):
> orthog(Q);
                                      true

# Alternatively,
> equal( map( simplify,evalm(transpose(Q)&*Q) ), diag(1,1) );
                                      true

# verifying that  transpose(Q)  is  inverse(Q)
\end{verbatim}

\pagebreak
\setcounter{theorem}{59}
\begin{theorem} If $A$ is an $m\times n$ matrix with orthonormal columns, then
    $$ A^T A = I$$
\end{theorem}

\begin{theorem}[and 4.62] If $Q$ is a square matrix, then the following are
  equivalent:
  \begin{enumerate}
  \item[]
    \begin{enumerate}
    \item $Q$ is orthogonal:
    \item $Q$ is invertible and $Q^{-1}=Q^T$;
    \item the rows of $Q$ are orthonormal.
    \item $Q$ preserves inner products:\qquad
      $(Q \b x)^T (Q \b y) = \b x^T Q^T Q \b y = \b x^T \b y$\quad for all
      vectors \b x,\b y;
    \item $Q$ preserves norms:\qquad
      $\Vert Q \b x\Vert = \Vert\b x\Vert$\quad for every vector \b x.
    \end{enumerate}
  \end{enumerate}
\end{theorem}

\setcounter{theorem}{62}
\begin{theorem}
  The product (and inverse) of orthogonal matrices is orthogonal.
\end{theorem}

\subsubsection*{The QR Decomposition}
We noted in Section 4.3 that the normal equation $A^T A\b x = A^T \b y$
associated with the least-squares problem is almost always ill-conditioned.
This is dealt with using $QR$-decomposition which is an alternative to
$LU$-decomposition.

\begin{theorem} If $A$ is an $m\times n$ matrix with linearly independent
  columns, then $A$ can be decomposed into the product
    $$ A = QR $$
  of an $m\times n$ matrix $Q$ with orthonormal columns and a non-singular
  upper triangular matrix $R$.\par
  Hence, if $A$ is an invertible (square) matrix, then\quad $A=QR$\quad
  where $Q$ is orthogonal and $R$ is a non-singular upper triangular matrix.
\end{theorem}

\begin{verbatim}
# A procedure to find the QR decomposition of a matrix using  GramSchmidt
with(linalg):
> qrfac:= proc(A,Q)              # See Hill Theorem 4.64, pp 283,284
>     local n,Qcols,R,i;
>     n:= coldim(A);
>     Qcols:= GramSchmidt( [col(A,1..n)] );
>     Q:= transpose( matrix(Qcols) );
#   Normalise Q by dividing each  col(Q,i)  by its norm.
>     for i from 1 to n do
>         Q:= mulcol( Q,i,1/norm(col(Q,i),2) )
>     od;
#   Since  Q &* R = A  and  transpose(Q) &* Q = Id,
>     R:= evalm(transpose(Q)&*A);
> end:
# Hill Section 4.5 Example 4, p284
# Find the QR decomposition of the matrix
> A:= matrix([ [1,0,0],
>              [1,1,0],
>              [1,1,1],
>              [1,1,1] ]):
\end{verbatim}
\pagebreak
\begin{verbatim}
> R:= qrfac( A,'Q'):                                 print ( Q,R);
      [             1/2                  ]
      [ 1/2  - 1/2 3            0        ]
      [                                  ]  [ 2     3/2          1       ]
      [            1/2          1/2  1/2 ]  [                            ]
      [ 1/2   1/6 3      - 1/3 2    3    ]  [         1/2          1/2   ]
      [                                  ], [ 0  1/2 3        1/3 3      ]
      [            1/2         1/2  1/2  ]  [                            ]
      [ 1/2   1/6 3       1/6 2    3     ]  [                   1/2  1/2 ]
      [                                  ]  [ 0      0     1/3 2    3    ]
      [            1/2         1/2  1/2  ]
      [ 1/2   1/6 3       1/6 2    3     ]

> equal( A, evalm(Q&*R) );
                                      true
\end{verbatim}
\noindent $QR$ decomposition is important in an algorithm -- the $QR$-algorithm
-- for numerical approximation of eigenvalues (see Hill Section 5.8).
The following shows its use in least squares.

\setcounter{theorem}{67}
\begin{theorem}
  Suppose $A$ is an $m\times n$ matrix with linearly independent columns.  Then
  the $QR$-decomposition $A=QR$ replaces the normal equation
  $A^T A\b x=A^T \b y$ by
    $$ R \b x = Q^T \b y$$
\end{theorem}

\begin{verbatim}
# Hill Section 4.5 Example 5, p285
> with(linalg):
# Use QR decomposition to find the line  y = mx + b  of best fit for the points
#   (2,3), (3,4), (4,4), (7,7)
#
# We want the least-squares solution for   A c = y   where c = [b,m] and
> A:= matrix([ [1,2], [1,3], [1,4], [1,7] ]):
> y:= vector( [3, 4, 4, 7] ):
> R:= qrfac(A,'Q'):            # See the function  qrfac  given above
> c:= linsolve( R, evalm(transpose(Q)&*y) ):
#
# The line of best fit is
> y = c[2]*x + c[1];
                                    11       19
                               y = ---- x + ----
                                    14       14
\end{verbatim}

\subsection{Encoding the QR Decomposition}
The procedure given above for $QR$-factorisation is inefficient and
{\em unstable} (i.e. if floating point arithmetic is used, then
round-off errors can build up as the calculation proceeds).
One way of avoiding this is the use of
Householder transformations.  Here is Maple code for a procedure giving the
$QR$-factorisation of matrices using these and a test for this procedure.

\begin{verbatim}
# Examples on Householder transformations and QR-decomposition based on them
> with(linalg):
\end{verbatim}
\pagebreak
\begin{verbatim}
> housetrans:= proc(v)       # See Hill Theorem 4.75, p292
>     local Id,a,w;
>     Id:= diag( 1$vectdim(v) ):
>     a:= - sign( evalf(v[1]) ):
#     (earlier versions of Maple had trouble with signs of square roots)
>     w:= copy(v):
>     w[1]:= v[1] - a*norm(v,2):
>     map( simplify, evalm( Id - 2*w&*transpose(w)/dotprod(w,w) ) ):
> end:
# Hill 4.6, Example 2(a), p292/293
> Q:= housetrans( vector([2,1,-2]) );
                                [ -2/3  -1/3   2/3 ]
                                [                  ]
                                [        14        ]
                                [ -1/3  ----  2/15 ]
                           Q := [        15        ]
                                [                  ]
                                [              11  ]
                                [  2/3  2/15  ---- ]
                                [              15  ]

> qrfac2:= proc(A,Q)               # See Hill pp 294,295 - returns the matrix R
>     local m,n,R,j,cj,Qj:
>     m:= rowdim(A):                    n:= coldim(A):
#   Assume that m < n for what follows
>     Q:= housetrans( col(A,1) );
>     R:= evalm( Q&*A );         # R is used in place of Hill's Aj
>     for j from 2 to n do
>         cj:= subvector(R,j..m,j); # The last m-j+1 entries of the j'th column
>         Qj:= diag( 1$(j-1), housetrans(cj) );
>         Q:= map( simplify, evalm(Q&*Qj) );        # update Q
>         R:= map( simplify, evalm(Qj&*R) );        # update R
>     od;
> end:
# Hill Section 4.5 Example 4, p284, (using Householder transformations)
> A:= matrix([ [1,0,0],
>              [1,1,0],
>              [1,1,1],
>              [1,1,1] ]):
> R:= qrfac2( A,'Q'):
# Because the terms get messy, we look at approximations to  Q,R:
> map( evalf,Q );
           [ -.5000000000   .8660254038        0             0      ]
           [                                                        ]
           [ -.5000000000  -.2886751346   .8164965810        0      ]
           [                                                        ]
           [ -.5000000000  -.2886751346  -.4082482905  -.7071067812 ]
           [                                                        ]
           [ -.5000000000  -.2886751346  -.4082482905   .7071067812 ]
\end{verbatim}
\pagebreak
\begin{verbatim}
> map( evalf,R );
                      [ -2.  -1.500000000       -1.     ]
                      [                                 ]
                      [  0   -.8660254038  -.5773502692 ]
                      [                                 ]
                      [  0         0       -.8164965810 ]
                      [                                 ]
                      [  0         0             0      ]

> map( evalf,evalm(Q&*R) );
                        [ 1.       0            0      ]
                        [                              ]
                        [ 1.  1.000000000       0      ]
                        [                              ]
                        [ 1.  1.000000000  .9999999998 ]
                        [                              ]
                        [ 1.  1.000000000  .9999999998 ]

# which is the same as  A  (to within the order of accuracy).
\end{verbatim}

\medskip \goodbreak
\headq{Notes}
\begin{enumerate}
\item This is {\em not} the same $QR$-decomposition as in Section 4.5, where
  the matrix $R$ was square (and invertible). In {\em this} version
  the matrix $Q$ is square (and orthogonal).
\item In the procedure {\tt housetrans}, the calculation was simplified by
  replacing the \b u of Hill's Theorem 4.75 by
  $$\b w = \b v - a\Vert\b v\Vert\b\epsilon_1,\qquad \mbox{and so}\qquad
                     \b u \b u^T = {\b w \b w^T\over \Vert\b w\Vert^2} $$
\item For any vector \b y, it is actually easier to evaluate $Q\b y$ directly
  since
    $$ Q \b y = (I - 2\b u \b u^T)\b y = \b y - 2\b u \b u^T\b y
              = \b y - 2(\b u\cdot \b y)\b u $$
  This is also useful for storage reasons. For more details, see page 296 of
  Hill.
\end{enumerate}
\headq{EXERCISE} (Optional, harder Maple) Create a new version of 
  {\tt housetrans} which returns the vector \b w instead of the matrix $Q$.
  Use this to write a new version of {\tt qrfac} which determines the matrix
  $R$ by evaluating each of its columns according to the formula given in (3)
  above.

\subsection{General Matrices of Linear Transformations;\ \ Similarity}
The concept (in Section 3.8) of the coordinates of a vector in a vector space
$V$ says in effect that {\em any} vector space of dimension $n$ may be regarded
as $\DoubleR^n$ (after a basis has been chosen). The same argument shows that,
relative to a chosen basis for the vector spaces $V$ and $W$, every linear
transformation $T:V\to W$ may be thought of as a matrix.

\setcounter{theorem}{81}
\begin{theorem}\quad Let $V$ and $W$ be two vector spaces with bases\quad
  $B=\{\rowi{\b v}{n}\}$\quad and\quad $C=\{\rowi{\b w}{m}\}$ respectively.
  Let $T:V\to W$ be a linear transformation. Then there is an $m\times n$
  matrix $[T]_{BC}$ -- whose columns are the coordinate vectors
  $[T(\b v_i)]_C$ -- such that
    $$ [T]_{BC} [\b x]_B = [T(\b x)]_C\qquad \mbox{for all}~~\b x\in V. $$
  $[T]_{BC}$ {\em is called} {\bf the matrix induced by $T$ relative to the
  bases $B$ and $C$}.
\end{theorem}

\begin{verbatim}
# Hill Section 4.7 Example 1, p298 and Example 5, p303
> with(linalg):
# Consider the linear transformation  T  from R^2 to R^2 where
> T:= proc(v) vector([ 2*v[1]+3*v[2], 4*v[1]+3*v[2] ]) end: # v[1]:=x: v[2]:=y: 
#
# Example 1:  Find the matrix of  T  relative to the (standard) basis  B  where
> e1:= vector([1,0]):      e2:= vector([0,1]):     B:= augment( e1,e2 ):
#
#   This can be done using the coordinates of  T(e1), T(e2)  as in Section 3.8:
> T_BB:= augment( linsolve(B,T(e1)), linsolve(B,T(e2)) );
                                        [ 2  3 ]
                                T_BB := [      ]
                                        [ 4  3 ]

# Now find the matrix of  T  relative to the basis  C  where
> v1:= vector([1,-1]):      v2:= vector([3,4]):     C:= augment( v1,v2 ):
> T_CC:= augment( linsolve(C,T(v1)), linsolve(C,T(v2)) );
                                       [ -1  0 ]
                               T_CC := [       ]
                                       [  0  6 ]

# Example 5:    The transition matrix from  C  to  B  is given by
> gaussjord( augment(B,C) ):
> P:= submatrix( ", 1..2, 3..4 );
                                      [  1  3 ]
                                 P := [       ]
                                      [ -1  4 ]

> equal( evalm(P^(-1) &* T_BB &* P), T_CC );
                                      true
\end{verbatim}
This illustrates
\setcounter{theorem}{83}
\begin{theorem} Let $V$ be a vector space with bases
  $B=\{\rowi{\b v}{n}\}$ and $C=\{\rowi{\b w}{n}\}$. Let
  $T:V\to V$ be a linear transformation. Then
    $$ [T]_{CC} = P^{-1} [T]_{BB} P $$
  where $P$ is the transition matrix from $C$ to $B$.
\end{theorem}

\headq{Warning} $P$ is the transition matrix from $C$ to $B$ --  not the
  transition matrix from $B$ to $C$. This is because function composition is
  from right to left, and not from left to right.

\begin{verbatim}
# Hill Section 4.7 Example 6, p304
> with(linalg):
# Consider the linear transformation  T  from P_2 to P_2 where
> T:= proc(f)      local a,b,c;
>     a:= coeff(f,x,0);    b:= coeff(f,x,1);    c:= coeff(f,x,2);
>     (5*a+b+c) + (b-a-c)*x + (a+b)*x^2;
> end:
\end{verbatim}
\pagebreak
\begin{verbatim}
# The matrix of  T  relative to the standard basis is
> T_S:= matrix([ [ 5, 1, 1],
>                [-1, 1,-1],
>                [ 1, 1, 0] ]):
# If
> f1:= 1-x:        f2:= 1+x:       f3:= 1+x+x^2:      C:= [ f1,f2,f3 ]:
# then the transition matrix from  C  to the standard basis is given by
> P:= matrix([ [ 1, 1, 1],
>              [-1, 1, 1],
>              [ 0, 0, 1] ]):
#
# and the matrix of  T  relative to the basis  C  is
> T_C:= evalm( P^(-1) &* T_S &* P );
                                      [ 3  3  4 ]
                                      [         ]
                               T_C := [ 1  1  1 ]
                                      [         ]
                                      [ 0  2  2 ]

#   This can be checked using the coordinates of
> T(f1),`    `, T(f2),`    `, T(f3);
                                        2                   2
                  4 - 2 x,     , 6 + 2 x ,     , 7 - x + 2 x

> simplify( T_C[1,1]*f1+T_C[2,1]*f2+T_C[3,1]*f3 - T(f1) ),
> simplify( T_C[1,2]*f1+T_C[2,2]*f2+T_C[3,2]*f3 - T(f2) ),
> simplify( T_C[1,3]*f1+T_C[2,3]*f2+T_C[3,3]*f3 - T(f3) );
                                    0, 0, 0
\end{verbatim}
\headq{Definition 4.90} Two $n\times n$ matrices $A$ and $B$ are {\bf similar}
  if there is an invertible $n\times n$ matrix $P$ such that $B = P^{-1} A P$.

\smallskip
\headq{EXERCISE} Let $A$, $B$, $C$ be three $n\times n$ matrices. Show that
\begin{enumerate}\item[]
  \begin{enumerate}
  \item $A$ is similar to $A$;
  \item If $A$ is similar to $B$, then $B$ is similar to $A$;
  \item If $A$ is similar to $B$, and $B$ is similar to $C$,
    then $A$ is similar to $C$.
  \end{enumerate}
\end{enumerate}

\setcounter{theorem}{90}
\begin{theorem}  Two $n\times n$ matrices are similar if and only if they
  represent the same linear transformation on an $n$-dimensional vector space
  (possibly relative to different bases).
\end{theorem}
Further examples of similar matrices will be given in the next chapter.

\vfill
\eject

