\documentclass[12pt]{report} \setlength{\textwidth}{6.5in} \setlength{\textheight}{8.5in} \setlength{\topmargin}{0.0in} \setlength{\oddsidemargin}{0.0in} \usepackage{graphics} \usepackage{curves} \begin{document} Version: \today \Huge \begin{center} \textbf{Notes for Applied Multivariate Analysis with MATLAB} \end{center} \Large \bigskip These notes were written for use in the Quantitative Psychology courses at the University of Illinois, Champaign. The expectation is that for Psychology 406/7 (Statistical Methods I and II), the material up through Section 0.1.12 be available to a student. For Multivariate Analysis (Psychology 594) and Covariance Structure and Factor Models (Psychology 588), the remainder of the notes are relevant, with particular emphasis on the Singular Value Decomposition (SVD) and the Eigenvector/Eigenvalue Decomposition (Spectral Decomposition). \newpage \tableofcontents \listoffigures \newpage \section{Necessary Matrix Algebra Tools} The strategies of multivariate analysis tend to be confusing unless specified compactly in matrix terms. Therefore, we will spend some significant amount of time on these topics because, in fact, most of multivariate analysis falls out directly once we have these tools under control. Remember the old Saturday Night Live skit with Hans and Franz, ``listen to me now, and believe me later.'' I have a goal in mind of where I would like you all to be --- at the point of understanding and being able to work with what is called the Singular Value Decomposition (SVD) of a matrix, and to understand the matrix topics that lead up to the SVD. Very much like teaching how to use some word-processing program, where we need to learn all the various commands and what they do, an introduction to the matrix tools can seem a little disjointed. But just like word-processing comes together more meaningfully when required to do your own manuscripts from beginning to end, once we proceed into the techniques of multivariate analysis per se, the wisdom of this preliminary matrix excursion will become apparent. \subsection{Preliminaries} A \emph{matrix} is merely an array of numbers; for example, \[ \left( \begin{array}{rrrr} 4 & -1 & 3 & 1 \\ 4 & 6 & 0 & 2 \\ 7 & 2 & 1 & 4 \\ \end{array} \right) \] \noindent is a matrix. In general, we denote a matrix by an uppercase (capital) boldface letter such as $\mathbf{A}$ (or using a proofreader representation on the blackboard, a capital letter with a wavy line underneath to indicate boldface): \[ \mathbf{A} = \left( \begin{array}{ccccc} a_{11} & a_{12} & a_{13} & \cdots & a_{1V} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2V} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{U1} & a_{U2} & a_{U3} & \cdots & a_{UV} \\ \end{array} \right) \] \noindent This matrix has $U$ rows and $V$ columns and is said to have \emph{order} $U \times V$. An arbitrary element $a_{uv}$ refers to the element in the $u^{th}$ row and $v^{th}$ column, with the row index always preceding the column index (and therefore, we might use the notation of $\mathbf{A} = \{a_{uv}\}_{U \times V}$ to indicate the matrix $\mathbf{A}$ as well as its order). \smallskip A $1 \times 1$ matrix such as $(4)_{1 \times 1}$ is just an ordinary number, called a \emph{scalar}. So without loss of any generality, numbers are just matrices. A \emph{vector} is a matrix with a single row or column; we denote a column vector by a lowercase boldface letter, e.g., $\mathbf{x}$, $\mathbf{y}$, $\mathbf{z}$, and so on. The vector \[ \mathbf{x} = \left( \begin{array}{c} x_{1} \\ \vdots \\ x_{U} \\ \end{array} \right)_{U \times 1} \] \noindent is of order $U \times 1$; the column indices are typically omitted since there is only one. A row vector is written as \[ \mathbf{x}' = (x_{1}, \ldots, x_{U})_{1 \times U} \] \noindent with the prime indicating the \emph{transpose} of $\mathbf{x}$, i.e., the interchange of row(s) and column(s). This transpose operation can be applied to any matrix; for example, \[ \mathbf{A} = \left( \begin{array}{rr} 1 & -1 \\ 3 & 7 \\ 4 & 1 \\ \end{array} \right)_{3 \times 2} \] \[ \mathbf{A}' = \left( \begin{array}{rrr} 1 & 3 & 4 \\ -1 & 7 & 1 \\ \end{array} \right)_{2 \times 3} \] If a matrix is \emph{square}, defined by having the same number of rows as columns, say $U$, and if the matrix and its transpose are equal, the matrix is said to be \emph{symmetric}. Thus, in $\mathbf{A} = \{a_{uv}\}_{U \times U}$, $a_{uv} = a_{vu}$ for all $u$ and $v$. As an example, \[ \mathbf{A} = \mathbf{A}' = \left( \begin{array}{rrr} 1 & 4 & 3 \\ 4 & 7 & -1 \\ 3 & -1 & 3 \\ \end{array} \right) \] \noindent For a square matrix $\mathbf{A}_{U \times U}$, the elements $a_{uu}$, $1 \le u \le U$, lie along the \emph{main} or \emph{principal} diagonal. The sum of main diagonal entries of a square matrix is called the \emph{trace}; thus, \[ \mathrm{trace}(\mathbf{A}_{U \times U}) \equiv \mathrm{tr}(\mathbf{A}) = a_{11} + \cdots + a_{UU} \] A number of special matrices appear periodically in the notes to follow. A $U \times V$ matrix of all zeros is called a \emph{null} matrix, and might be denoted by \[ \mathbf{\emptyset} = \left( \begin{array}{ccc} 0 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & 0 \\ \end{array} \right) \] \noindent Similarly, we might at times need an $U \times V$ matrix of all ones, say $\mathbf{E}$: \[ \mathbf{E} = \left( \begin{array}{ccc} 1 & \cdots & 1 \\ \vdots & \ddots & \vdots \\ 1 & \cdots & 1 \\ \end{array} \right) \] \noindent A \emph{diagonal} matrix is square with zeros in all the off main-diagonal positions: \[ \mathbf{D}_{U \times U} = \left( \begin{array}{ccc} a_{1} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & a_{U} \\ \end{array} \right)_{U \times U} \] \noindent Here, we again indicate the main diagonal entries with just one index as $a_{1}, a_{2}, \ldots, a_{U}$. If all of the main diagonal entries in a diagonal matrix are $1$s, we have the \emph{identity} matrix denoted by $\mathbf{I}$: \[ \mathbf{I} = \left( \begin{array}{cccc} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \\ \end{array} \right) \] To introduce some useful operations on matrices, suppose we have two matrices $\mathbf{A}$ and $\mathbf{B}$ of the same $U \times V$ order: \[ \mathbf{A} = \left( \begin{array}{ccc} a_{11} & \cdots & a_{1V} \\ \vdots & \ddots & \vdots \\ a_{U1} & \cdots & a_{UV} \\ \end{array} \right)_{U \times V} \] \[ \mathbf{B} = \left( \begin{array}{ccc} b_{11} & \cdots & b_{1V} \\ \vdots & \ddots & \vdots \\ b_{U1} & \cdots & b_{UV} \\ \end{array} \right)_{U \times V} \] \noindent As a definition for equality of two matrices of the same order (and for which it only makes sense to talk about equality), we have: $\mathbf{A}$ = $\mathbf{B}$ if and only if $a_{uv}$ = $b_{uv}$ for all $u$ and $v$. \noindent Remember, the ``if and only if'' statement (sometimes abbreviated as ``iff'') implies two conditions: if $\mathbf{A}$ = $\mathbf{B}$, then $a_{uv}$ = $b_{uv}$ for all $u$ and $v$; if $a_{uv}$ = $b_{uv}$ for all $u$ and $v$, then $\mathbf{A}$ = $\mathbf{B}$. \noindent Any definition by its very nature implies an ``if and only if'' statement. \smallskip To add two matrices together, they first have to be of the same order (referred to as conformable for addition); we then do the addition component by component: \[ \mathbf{A} + \mathbf{B} = \left( \begin{array}{ccc} a_{11} + b_{11} & \cdots & a_{1V} + b_{1V}\\ \vdots & \ddots & \vdots \\ a_{U1} + b_{U1} & \cdots & a_{UV} + b_{UV} \\ \end{array} \right)_{U \times V} \] \indent To preform scalar multiplication of a matrix $\mathbf{A}$ by, say, a constant $c$, we again do the multiplication component by component: \[ c\mathbf{A} = \left( \begin{array}{ccc} ca_{11} & \cdots & ca_{1V} \\ \vdots & \ddots & \vdots \\ ca_{U1} & \cdots & ca_{UV} \\ \end{array} \right) = c\left( \begin{array}{ccc} a_{11} & \cdots & a_{1V} \\ \vdots & \ddots & \vdots \\ a_{U1} & \cdots & a_{UV} \\ \end{array} \right) \] \noindent Thus, if one wished to define the difference of two matrices, we could proceed rather obviously as follows: \[ \mathbf{A} - \mathbf{B} \equiv \mathbf{A} + (-1)\mathbf{B} = \{a_{uv} - b_{uv}\} \] One of the more important matrix operations is multiplication where two matrices are said to be conformable for multiplication if the number of rows in one matches the number of columns in the second. For example, suppose $\mathbf{A}$ is $U \times V$ and $\mathbf{B}$ is $V \times W$; then, because the number of columns in $\mathbf{A}$ matches the number of rows in $\mathbf{B}$, we can define $\mathbf{AB}$ as $\mathbf{C}_{U \times W}$, where $\{c_{uw}\}$ = $\{\sum_{k=1}^{V} a_{uk}b_{kw} \}$. This process might be referred to as row (of $\mathbf{A}$) by column (of $\mathbf{B}$) multiplication; the following simple example should make this clear: \[ \mathbf{A}_{3 \times 2} = \left( \begin{array}{rr} 1 & 4 \\ 3 & 1 \\ -1 & 0 \\ \end{array} \right) , \; \mathbf{B}_{2 \times 4} = \left( \begin{array}{rrrr} -1 & 2 & 0 & 1 \\ 1 & 0 & 1 & 4 \\ \end{array} \right) ; \] \[ \mathbf{A}\mathbf{B} = \mathbf{C}_{3 \times 4} = \] \[ \left( \begin{array}{rrrr} 1(-1) + 4(1) & 1(2) + 4(0) & 1(0) + 4(1) & 1(1) + 4(4) \\ 3(-1) + 1(1) & 3(2) + 1(0) & 3(0) + 1(1) & 3(1) + 1(4) \\ -1(-1) + 0(1) & -1(2) + 0(0) & -1(0) + 0(1) & -1(1) + 0(4) \\ \end{array} \right) = \] \[ \left( \begin{array}{rrrr} 3 & 2 & 4 & 17 \\ -2 & 6 & 1 & 7 \\ 1 & -2 & 0 & -1 \\ \end{array} \right) \] Some properties of matrix addition and multiplication follow, where the matrices are assumed conformable for the operations given: \smallskip (A) matrix addition is commutative: \[ \mathbf{A} + \mathbf{B} = \mathbf{B} + \mathbf{A} \] (B) matrix addition is associative: \[ \mathbf{A} + (\mathbf{B} + \mathbf{C}) = (\mathbf{A} + \mathbf{B}) + \mathbf{C} \] (C) matrix multiplication is right and left distributive over matrix addition: \[ \mathbf{A}(\mathbf{B} + \mathbf{C}) = \mathbf{A}\mathbf{B} + \mathbf{A}\mathbf{C} \] \[(\mathbf{A} + \mathbf{B})\mathbf{C} = \mathbf{A}\mathbf{C} + \mathbf{B}\mathbf{C} \] (D) matrix multiplication is associative: \[\mathbf{A}(\mathbf{B}\mathbf{C}) = (\mathbf{A}\mathbf{B})\mathbf{C} \] \noindent In general, $\mathbf{A}\mathbf{B} \ne \mathbf{B}\mathbf{A}$ even if both products are defined. Thus, multiplication is not commutative as the following simple example shows: \[ \mathbf{A}_{2 \times 2} = \left( \begin{array}{cc} 0 & 1 \\ 1 & 0 \\ \end{array} \right); \; \mathbf{B}_{2 \times 2} = \left( \begin{array}{cc} 1 & 1 \\ 0 & 1 \\ \end{array} \right); \; \mathbf{A}\mathbf{B} = \left( \begin{array}{cc} 0 & 1 \\ 1 & 1 \\ \end{array} \right); \; \mathbf{B}\mathbf{A} = \left( \begin{array}{cc} 1 & 1 \\ 1 & 0 \\ \end{array} \right) \] In the product $\mathbf{A}\mathbf{B}$, we say that $\mathbf{B}$ is \emph{premultiplied} by $\mathbf{A}$ and $\mathbf{A}$ is \emph{postmultiplied} by $\mathbf{B}$. Thus, if we pre- or postmultiply a matrix by the identity, the same matrix is retrieved: \[ \mathbf{I}_{U \times U}\mathbf{A}_{U \times V} = \mathbf{A}_{U \times V}; \; \mathbf{A}_{U \times V}\mathbf{I}_{V \times V} = \mathbf{A}_{U \times V} \] \noindent If we premultiply $\mathbf{A}$ by a diagonal matrix $\mathbf{D}$, then each row of $\mathbf{A}$ is multiplied by a particular diagonal entry in $\mathbf{D}$: \[ \mathbf{D}_{U \times U}\mathbf{A}_{U \times V} = \left( \begin{array}{ccc} d_{1}a_{11} & \cdots & d_{1}a_{1V} \\ \vdots & \ddots & \vdots \\ d_{U}a_{U1} & \cdots & d_{U}a_{UV} \\ \end{array} \right) \] \noindent If $\mathbf{A}$ is post-multiplied by a diagonal matrix $\mathbf{D}$, then each column of $\mathbf{A}$ is multiplied by a particular diagonal entry in $\mathbf{D}$: \[ \mathbf{A}_{U \times V}\mathbf{D}_{V \times V} = \left( \begin{array}{ccc} d_{1}a_{11} & \cdots & d_{V}a_{1V} \\ \vdots & \ddots & \vdots \\ d_{1}a_{U1} & \cdots & d_{V}a_{UV} \\ \end{array} \right) \] \noindent Finally, we end this section with a few useful results on the transpose operation and matrix multiplication and addition: \[ (\mathbf{A}\mathbf{B})' = \mathbf{B}'\mathbf{A}'; \; (\mathbf{A}\mathbf{B}\mathbf{C})' = \mathbf{C}'\mathbf{B}'\mathbf{A}'; \; \dots \] \[ (\mathbf{A}')' = \mathbf{A}; \; (\mathbf{A} + \mathbf{B})' = \mathbf{A}' + \mathbf{B}' \] \subsection{The Data Matrix} A very common type of matrix encountered in multivariate analysis is what is referred to as a data matrix containing, say, observations for $N$ subjects on $P$ variables. We will typically denote this matrix by $\mathbf{X}_{N \times P}$ = $\{x_{ij}\}$, with a generic element of $x_{ij}$ referring to the observation for subject or row $i$ on variable or column $j$ ($1 \le i \le N$ and $1 \le j \le P$): \[ \mathbf{X}_{N \times P} = \left( \begin{array}{cccc} x_{11} & x_{12} & \cdots & x_{1P} \\ x_{21} & x_{22} & \cdots & x_{2P} \\ \vdots & \vdots & \ddots & \vdots \\ x_{N1} & x_{N2} & \cdots & x_{NP} \\ \end{array} \right) \] \noindent All right-thinking people always list subjects as rows and variables as columns, conforming also to the now-common convention for computer spreadsheets. \smallskip Any matrix in general, including a data matrix, can be viewed either as a collection of its row vectors or of its column vectors, and these interpretations can be generally useful. For a data matrix $\mathbf{X}_{N \times P}$, let $\mathbf{x}_{i}' = (x_{i1}, \ldots, x_{iP})_{1 \times P}$ denote the row vector for subject $i$, $1 \le i \le N$, and let $\mathbf{v}_{j}$ denote the $N \times 1$ column vector for variable $j$: \[ \mathbf{v}_{j} = \left( \begin{array}{c} x_{1j} \\ \vdots \\ x_{Nj} \\ \end{array} \right)_{N \times 1} \] \noindent Thus, each subject could be viewed as providing a vector of coordinates ($1 \times P$) in $P$-dimensional ``variable space,'' where the $P$ axes correspond to the $P$ variables; or each variable could be viewed as providing a vector of coordinates ($N \times 1$) in ``subject space,'' where the $N$ axes correspond to the $N$ subjects: \[ \mathbf{X}_{N \times P} = \left( \begin{array}{c} \mathbf{x}_{1}' \\ \mathbf{x}_{2}' \\ \vdots \\ \mathbf{x}_{N}' \\ \end{array} \right) = \left( \begin{array}{cccc} \mathbf{v}_{1} & \mathbf{v}_{2} & \cdots & \mathbf{v}_{P} \\ \end{array} \right) \] \subsection{Inner Products} The \emph{inner product} (also called the dot or scalar product) of two vectors , $\mathbf{x}_{U \times 1}$ and $\mathbf{y}_{U \times 1}$, is defined as \[ \mathbf{x}'\mathbf{y} = (x_{1}, \ldots, x_{U})\left( \begin{array}{c} y_{1} \\ \vdots \\ y_{U} \\ \end{array} \right) = \sum_{u = 1}^{U} x_{u}y_{u} \] \noindent Thus, the inner product of a vector with itself is merely the sum of squares of the entries in the vector: $\mathbf{x}'\mathbf{x} = \sum_{u = 1}^{U} x_{u}^{2}$. Also, because an inner product is a scalar and must equal it own transpose (i.e., $\mathbf{x}'\mathbf{y} = (\mathbf{x}'\mathbf{y})' = \mathbf{y}'\mathbf{x}$), we have the end result that \[ \mathbf{x}'\mathbf{y} = \mathbf{y}'\mathbf{x} \] \noindent If there is an inner product, there should also be an \emph{outer product} defined as the $U \times U$ matrices given by $\mathbf{x}\mathbf{y}'$ or as $\mathbf{y}\mathbf{x}'$. As indicated by the display equations below, $\mathbf{x}\mathbf{y}'$ is the transpose of $\mathbf{y}\mathbf{x}'$: \[ \mathbf{x}\mathbf{y}' = \left( \begin{array}{c} x_{1} \\ \vdots \\ x_{U} \\ \end{array} \right)(y_{1}, \dots, y_{N}) = \left( \begin{array}{ccc} x_{1}y_{1} & \cdots & x_{1}y_{U} \\ \vdots & \dots & \vdots \\ x_{U}y_{1} & \cdots & x_{U}y_{U} \\ \end{array} \right) \] \[ \mathbf{y}\mathbf{x}' = \left( \begin{array}{c} y_{1} \\ \vdots \\ y_{U} \\ \end{array} \right)(x_{1}, \dots, x_{U}) = \left( \begin{array}{ccc} y_{1}x_{1} & \cdots & y_{1}x_{U} \\ \vdots & \dots & \vdots \\ y_{U}x_{1} & \cdots & y_{U}x_{U} \\ \end{array} \right) \] A vector can be viewed as a geometrical vector in $U$ dimensional space. Thus, the two $2 \times 1$ vectors \[ \mathbf{x} = \left( \begin{array}{c} 3 \\ 4 \\ \end{array} \right) ; \; \mathbf{y} = \left( \begin{array}{c} 4 \\ 1 \\ \end{array} \right) \] \noindent can be represented in the two-dimensional Figure 1 below, with the entries in the vectors defining the coordinates of the endpoints of the arrows. \smallskip \begin{figure} \setlength{\unitlength}{1cm} \begin{picture}(9,9) \put(0,0){\vector(1,0){6}} \put(0,0){\vector(0,1){6}} \put(0,0){\vector(3,4){3}} \put(0,0){\vector(4,1){4}} \put(3,4.1){\makebox(0,0)[br]{(3,4)}} \put(4,1.1){\makebox(0,0)[br]{(4,1)}} \put(-.2,0){\makebox(0,0)[tr]{(0,0)}} \put(3,2){\makebox(0,0)[b]{$\theta$}} \put(0,0){\arc(3,.75){39}} \put(1,-.1){\line(0,1){.2}} \put(2,-.1){\line(0,1){.2}} \put(3,-.1){\line(0,1){.2}} \put(4,-.1){\line(0,1){.2}} \put(-.1,1){\line(1,0){.2}} \put(-.1,2){\line(1,0){.2}} \put(-.1,3){\line(1,0){.2}} \put(-.1,4){\line(1,0){.2}} \put(1,-.2){\makebox(0,0)[t]{1}} \put(2,-.2){\makebox(0,0)[t]{2}} \put(3,-.2){\makebox(0,0)[t]{3}} \put(4,-.2){\makebox(0,0)[t]{4}} \put(-.2,1){\makebox(0,0)[r]{1}} \put(-.2,2){\makebox(0,0)[r]{2}} \put(-.2,3){\makebox(0,0)[r]{3}} \put(-.2,4){\makebox(0,0)[r]{4}} \end{picture} \caption{Two vectors plotted in two-dimensional space} \end{figure} The \emph{Euclidean distance} between two vectors, $\mathbf{x}$ and $\mathbf{y}$, is given as: \[ \sqrt{\sum_{u=1}^{U} (x_{u} - y_{u})^{2}} = \sqrt{(\mathbf{x} - \mathbf{y})'(\mathbf{x} - \mathbf{y})} \] \noindent and the \emph{length} of any vector is the Euclidean distance between the vector and the origin. Thus, in Figure 1, the distance between $\mathbf{x}$ and $\mathbf{y}$ is $\sqrt{10}$ with respective lengths of $5$ and $\sqrt{17}$. \smallskip The cosine of the angle between the two vectors $\mathbf{x}$ and $\mathbf{y}$ is defined by: \[ \cos(\theta) = \frac{\mathbf{x}'\mathbf{y}}{(\mathbf{x}'\mathbf{x})^{1/2}(\mathbf{y}'\mathbf{y})^{1/2}} \] \noindent Thus, in the figure we have \[\cos(\theta) = \frac{\left( \begin{array}{cc} 3 & 4 \\ \end{array} \right) \left( \begin{array}{c} 4 \\ 1 \\ \end{array} \right)}{5\sqrt{17}} = \frac{16}{5\sqrt{17}} = .776 \] \noindent The cosine value of $.776$ corresponds to an angle of $39.1$ degrees or $.68$ radians; these later values can be found with the inverse (or arc) cosine function (on, say, a hand calculator, or using MATLAB as we suggest in the next section). \smallskip When the means of the entries in $\mathbf{x}$ and $\mathbf{y}$ are zero (i.e., deviations from means have been taken), then $\cos(\theta)$ is the correlation between the entries in the two vectors. Vectors at right angles have $\cos(\theta) = 0$, or alternatively, the correlation is zero. \begin{figure} \setlength{\unitlength}{1cm} \begin{picture}(9,9)(-2,-2) \put(0,0){\vector(1,1){4}} \put(0,0){\vector(1,0){6}} \put(0,0){\vector(1,0){4}} \put(4,-1){\vector(0,1){1}} \put(4,0){\line(0,1){4}} \put(2,2.1){\makebox(0,0)[r]{$c^2$}} \put(4.1,2){\makebox(0,0)[l]{$b^2$}} \put(2,-.1){\makebox(0,0)[t]{$a^2$}} \put(0,0){\arc(2,0){45}} \put(4,4.1){\makebox(0,0)[b]{$\mathbf{x}$}} \put(6.1,0){\makebox(0,0)[l]{$\mathbf{y}$}} \put(4,-1){\makebox(0,0)[t]{$d\mathbf{y}$}} \put(2,1){\makebox(0,0){$\theta$}} \end{picture} \caption{Illustration of projecting one vector onto another} \end{figure} \smallskip Figure 2 shows two generic vectors, $\mathbf{x}$ and $\mathbf{y}$, where without loss of any real generality, $\mathbf{y}$ is drawn horizontally in the plane and $\mathbf{x}$ is projected at a right angle onto the vector $\mathbf{y}$ resulting in a point defined as a multiple $d$ of the vector $\mathbf{y}$. The formula for $d$ that we demonstrate below is based on the Pythagorean theorem that $c^2 = b^2 + a^2$: \[ c^2 = b^2 + a^2 \Rightarrow \mathbf{x}'\mathbf{x} = (\mathbf{x} - d\mathbf{y})'(\mathbf{x} - d\mathbf{y}) + d^2\mathbf{y}'\mathbf{y} \Rightarrow \] \[ \mathbf{x}'\mathbf{x} = \mathbf{x}'\mathbf{x} - d\mathbf{x}'\mathbf{y} - d\mathbf{y}'\mathbf{x} + d^2\mathbf{y}'\mathbf{y} + d^2\mathbf{y}'\mathbf{y} \Rightarrow \] \[0 = -2d\mathbf{x}'\mathbf{y} + 2d^2\mathbf{y}'\mathbf{y} \Rightarrow \] \[ d = \frac{\mathbf{x}'\mathbf{y}}{\mathbf{y}'\mathbf{y}} \] \noindent The diagram in Figure 2 is somewhat constricted in the sense that the angle between the vectors shown is less than $90$ degrees; this allows the constant $d$ to be positive. Other angles might lead to negative $d$ when defining the projection of $\mathbf{x}$ onto $\mathbf{y}$, and would merely indicate the need to consider the vector $\mathbf{y}$ oriented in the opposite (negative) direction. Similarly, the vector $\mathbf{y}$ is drawn with a larger length than $\mathbf{x}$ which gives a value for $d$ that is less than $1.0$; otherwise, $d$ would be greater than $1.0$, indicating a need to stretch $\mathbf{y}$ to represent the point of projection onto it. \smallskip There are other formulas possible based on this geometric information: the length of the projection is merely $d$ times the length of $\mathbf{y}$; and $\cos(\theta)$ can be given as the length of $d\mathbf{y}$ divided by the length of $\mathbf{x}$, which is $d\sqrt{\mathbf{y}'\mathbf{y}}/\sqrt{\mathbf{x}'\mathbf{x}} = \mathbf{x}'\mathbf{y}/(\sqrt{\mathbf{x}'\mathbf{x}}\sqrt{\mathbf{y}'\mathbf{y}})$. \subsection{Determinants} To each square matrix, $\mathbf{A}_{U \times U}$, there is an associated scalar called the \emph{determinant} of $\mathbf{A}$ that is denoted by $|\mathbf{A}|$ or $\det(\mathbf{A})$. Determinants up to a $3 \times 3$ can be given by formula: \[ \det(\left( \begin{array}{c} a \\ \end{array} \right)_{1 \times 1}) = a; \; \det(\left( \begin{array}{cc} a & b \\ c & d \\ \end{array} \right)_{2 \times 2}) = ad - bc; \] \[ \det(\left( \begin{array}{ccc} a & b & c \\ d & e & f \\ g & h & i \\ \end{array} \right)_{3 \times 3}) = aei + dhc + gfb -(ceg + fha + idb) \] \noindent Beyond a $3 \times 3$ we can use a recursive process illustrated below. This requires the introduction of a few additional matrix terms that we now give: for a square matrix $\mathbf{A}_{U \times U}$, define $\mathbf{A}_{uv}$ to be the $(n-1) \times (n-1)$ submatrix of $\mathbf{A}$ constructed by deleting the $u^{th}$ row and $v^{th}$ column of $\mathbf{A}$. We call $\det(\mathbf{A}_{uv})$ the \emph{minor} of the entry $a_{uv}$; the signed minor of $(-1)^{u+v}\det(\mathbf{A}_{uv})$ is called the \emph{cofactor} of $a_{uv}$. The recursive algorithm would chose some row or column (rather arbitrarily), and find the cofactors for the entries in it; the cofactors would then be weighted by the relevant entries and summed. \smallskip As an example, consider the $4 \times 4$ matrix \[ \left( \begin{array}{rrrr} 1 & -1 & 3 & 1 \\ -1 & 1 & 0 & -1 \\ 3 & 2 & 1 & 2 \\ 1 & 2 & 4 & 3 \\ \end{array} \right) \] \noindent and choose the second row. The expression below involves the weighted cofactors for $3 \times 3$ submatrices that can be obtained by formulas. Beyond a $4 \times 4$ there will be nesting of the processes: \[ (-1)((-1)^{2+1})\det(\left( \begin{array}{rrr} -1 & 3 & 1 \\ 2 & 1 & 2 \\ 2 & 4 & 3 \\ \end{array} \right)) + (1)((-1)^{2+2})\det(\left( \begin{array}{rrr} 1 & 3 & 1 \\ 3 & 1 & 2 \\ 1 & 4 & 3 \\ \end{array} \right)) + \] \[ (0)((-1)^{2+3})\det(\left( \begin{array}{rrr} 1 & -1 & 1 \\ 3 & 2 & 2 \\ 1 & 2 & 3 \\ \end{array} \right)) + (-1)((-1)^{2+4})\det(\left( \begin{array}{rrr} 1 & -1 & 3 \\ 3 & 2 & 1 \\ 1 & 2 & 4 \\ \end{array} \right)) = \] \[ 5 + (-15) + 0 + (-29) = -39 \] Another strategy to find the determinant of a matrix is to reduce it a form in which we might note the determinant more or less by simple inspection. The reductions could be carried out by operations that have a known effect on the determinant; the form which we might seek is a matrix that is either \emph{upper-triangular} (all entries below the main diagonal are all zero), \emph{lower-triangular} (all entries above the main diagonal are all zero), or diagonal. In these latter cases, the determinant is merely the product of the diagonal elements. Once found, we can note how the determinant might have been changed by the reduction process and carry out the reverse changes to find the desired determinant. \smallskip The properties of determinants that we could rely on in the above iterative process are as follows: \smallskip \noindent (A) if \emph{one} row of $\mathbf{A}$ is multiplied by a constant $c$, the new determinant is $c\det(\mathbf{A})$; the same is true for multiplying a column by $c$; \smallskip \noindent (B) if two rows or two columns of a matrix are interchanged, the sign of the determinant is changed; \smallskip \noindent (C) if two rows or two columns of a matrix are equal, the determinant is zero; \smallskip \noindent (D) the determinant is unchanged by adding a multiple of some row to another row; the same is true for columns; \smallskip \noindent (E) a zero row or column implies a zero determinant; \smallskip \noindent (F) $\det(\mathbf{A}\mathbf{B}) = \det(\mathbf{A})\det(\mathbf{B})$ \subsection{Linear Independence/Dependence of Vectors} Suppose I have a collection of $K$ vectors each of size $U \times 1$, $\mathbf{x}_{1},\ldots,\mathbf{x}_{K}$. If no vector in the set can be written as a linear combination of the remaining ones, the set of vectors is said to be \emph{linearly independent}; otherwise, the vectors are \emph{linearly dependent}. As an example, consider the three vectors: \[ \mathbf{x}_{1} = \left( \begin{array}{r} 1 \\ 4 \\ 0 \\ \end{array} \right); \; \mathbf{x}_{2} = \left( \begin{array}{r} 1 \\ -1 \\ 1 \\ \end{array} \right); \; \mathbf{x}_{3} = \left( \begin{array}{r} 3 \\ 7 \\ 1 \\ \end{array} \right) \] \noindent Because $2\mathbf{x}_{1} + \mathbf{x}_{2} = \mathbf{x}_{3}$, we have a linear dependence among the three vectors; however, $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$, or, $\mathbf{x}_{2}$ and $\mathbf{x}_{3}$, are linearly independent. \smallskip If the $U$ vectors (each of size $U \times 1$), $\mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{U}$, are linearly independent, then the collection defines a \emph{basis}, i.e., any vector can be written as a linear combination of $\mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{U}$. For example, using the \emph{standard basis}, $\mathbf{e}_{1}, \mathbf{e}_{2}, \ldots, \mathbf{e}_{U}$, where $\mathbf{e}_{u}$ is a vector of all zeros except for a single one in the $u^{th}$ position, any vector $\mathbf{x}' = (x_{1},\ldots,x_{U})$ can be written as: \[ \left( \begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{U} \\ \end{array} \right) = x_{1}\left( \begin{array}{c} 1 \\ 0 \\ \vdots \\ 0 \\ \end{array} \right) + x_{2}\left( \begin{array}{c} 0 \\ 1 \\ \vdots \\ 0 \\ \end{array} \right) + \cdots + x_{U}\left( \begin{array}{c} 0 \\ 0 \\ \vdots \\ 1 \\ \end{array} \right) = \] \[ x_{1}\mathbf{e}_{1} + x_{2}\mathbf{e}_{2} + \cdots + x_{U}\mathbf{e}_{U} \] \noindent Bases that consist of \emph{orthogonal} vectors (where all inner products are zero) are important later in what is known as principal components analysis. The standard basis involves orthogonal vectors, and any other basis may always be modified by what is called the Gram-Schmidt orthogonalization process to produce a new basis that does contain all orthogonal vectors. \subsection{Matrix Inverses} Suppose $\mathbf{A}$ and $\mathbf{B}$ are both square and of size $U \times U$. If $\mathbf{A}\mathbf{B} = \mathbf{I}$, then $\mathbf{B}$ is said to be an inverse of $\mathbf{A}$ and is denoted by $\mathbf{A}^{-1} (\equiv \mathbf{B})$. Also, if $\mathbf{A}\mathbf{A}^{-1} = \mathbf{I}$, then $\mathbf{A}^{-1}\mathbf{A} = \mathbf{I}$ holds automatically. If $\mathbf{A}^{-1}$ exists, the matrix $\mathbf{A}$ is said to be \emph{nonsingular}; if $\mathbf{A}^{-1}$ does not exist, $\mathbf{A}$ is \emph{singular}. \smallskip An example: \[ \left( \begin{array}{rr} 1 & 3 \\ 2 & 1 \\ \end{array} \right) \left( \begin{array}{rr} -1/5 & 3/5 \\ 2/5 & -1/5 \\ \end{array} \right) = \left( \begin{array}{rr} 1 & 0 \\ 0 & 1 \\ \end{array} \right) \] \[ \left( \begin{array}{rr} -1/5 & 3/5 \\ 2/5 & -1/5 \\ \end{array} \right) \left( \begin{array}{rr} 1 & 3 \\ 2 & 1 \\ \end{array} \right) = \left( \begin{array}{rr} 1 & 0 \\ 0 & 1 \\ \end{array} \right) \] Given a matrix $\mathbf{A}$, the inverse $\mathbf{A}^{-1}$ can be found using the following four steps: \smallskip (A) form a matrix of the same size as $\mathbf{A}$ containing the minors for all entries of $\mathbf{A}$; \smallskip (B) multiply the matrix of minors by $(-1)^{u+v}$ to produce the matrix of cofactors; \smallskip (C) divide all entries in the cofactors matrix by $\det(\mathbf{A}$); \smallskip (D) the transpose of the matrix found in (C) gives $\mathbf{A}^{-1}$. \smallskip \noindent As a mnemonic device to remember these four steps, we have the phrase ``My Cat Does Tricks'' for Minor, Cofactor, Determinant Division, Transpose'' (I tried to work ``my cat turns tricks'' into the appropriate phrase but failed with the second to the last ``t''). Obviously, an inverse exists for a matrix $\mathbf{A}$ if $\det(\mathbf{A}) \ne 0$, allowing the division in step (C) to take place. \smallskip An example: for \[ \mathbf{A} = \left( \begin{array}{rrr} 1 & 3 & 2 \\ 0 & 1 & 1 \\ 0 & 2 & 1 \\ \end{array} \right); \; \det(\mathbf{A}) = -1 \] Step (A), the matrix of minors: \[ \left( \begin{array}{rrr} -1 & 0 & 0 \\ -1 & 1 & 2 \\ 1 & 1 & 1 \\ \end{array} \right) \] Step (B), the matrix of cofactors: \[ \left( \begin{array}{rrr} -1 & 0 & 0 \\ 1 & 1 & -2 \\ 1 & -1 & 1 \\ \end{array} \right) \] Step (C), determinant division: \[ \left( \begin{array}{rrr} 1 & 0 & 0 \\ -1 & -1 & 2 \\ -1 & 1 & -1 \\ \end{array} \right) \] Step (D), matrix transpose: \[ \mathbf{A}^{-1} = \left( \begin{array}{rrr} 1 & -1 & -1 \\ 0 & -1 & 1 \\ 0 & 2 & -1 \\ \end{array} \right) \] We can easily verify that $\mathbf{A}\mathbf{A}^{-1} = \mathbf{I}$: \[\left( \begin{array}{rrr} 1 & 3 & 2 \\ 0 & 1 & 1 \\ 0 & 2 & 1 \\ \end{array} \right) \left( \begin{array}{rrr} 1 & -1 & -1 \\ 0 & -1 & 1 \\ 0 & 2 & -1 \\ \end{array} \right) = \left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) \] \noindent As a very simple instance of the mnemonic in the case of a $2 \times 2$ matrix with arbitrary entries: \[ \mathbf{A} = \left( \begin{array}{rr} a & b \\ c & d \\ \end{array} \right) \] \noindent the inverse exists if $\det(\mathbf{A}) = ad - bc \ne 0$: \[ \mathbf{A}^{-1} = \frac{1}{ad - bc}\left( \begin{array}{rr} d & -b \\ -c & a \\ \end{array} \right) \] Several properties of inverses are given below that will prove useful in our continuing presentation: \smallskip (A) if $\mathbf{A}$ is symmetric, then so is $\mathbf{A}^{-1}$; \smallskip (B) $(\mathbf{A}')^{-1} = (\mathbf{A}^{-1})'$; or, the inverse of a transpose is the transpose of the inverse; \smallskip (C) $(\mathbf{A}\mathbf{B})^{-1} = \mathbf{B}^{-1}\mathbf{A}^{-1}$; $(\mathbf{A}\mathbf{B}\mathbf{C})^{-1} = \mathbf{C}^{-1}\mathbf{B}^{-1}\mathbf{A}^{-1}$; or, the inverse of a product is the product of inverses in the opposite order; \smallskip (D) $(c\mathbf{A})^{-1} = (\frac{1}{c})\mathbf{A}^{-1}$; or, the inverse of a scalar times a matrix is the scalar inverse times the matrix inverse; \smallskip (E) the inverse of a diagonal matrix, is also diagonal with the entries being the inverses of the entries from the original matrix (assuming none are zero): \[ \left( \begin{array}{ccc} a_{1} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & a_{U} \\ \end{array} \right)^{-1} = \left( \begin{array}{ccc} \frac{1}{a_{1}} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \frac{1}{a_{U}} \\ \end{array} \right) \] \subsection{Matrices as Transformations} Any $U \times V$ matrix $\mathbf{A}$ can be seen as transforming a $V \times 1$ vector $\mathbf{x}_{V \times 1}$ to another $U \times 1$ vector $\mathbf{y}_{U \times 1}$: \[ \mathbf{y}_{U \times 1} = \mathbf{A}_{U \times V}\mathbf{x}_{V \times 1}\] \noindent or, \[ \left( \begin{array}{c} y_{1} \\ \vdots \\ y_{U} \\ \end{array} \right) = \left( \begin{array}{ccc} a_{11} & \cdots & a_{1V} \\ \vdots & \ddots & \vdots \\ a_{U1} & \cdots & a_{UV} \\ \end{array} \right) \left( \begin{array}{c} x_{1} \\ \vdots \\ x_{V} \\ \end{array} \right) \] \noindent where $y_{u} = a_{u1}x_{1} + a_{u2}x_{2} + \cdots + a_{uV}x_{V}$. Alternatively, $\mathbf{y}$ can be written as a linear combination of the columns of $\mathbf{A}$ with weights given by $x_{1},\ldots,x_{V}$: \[ \left( \begin{array}{c} y_{1} \\ \vdots \\ y_{U} \\ \end{array} \right) = x_{1}\left( \begin{array}{c} a_{11} \\ \vdots \\ a_{U1} \\ \end{array} \right) + x_{2}\left( \begin{array}{c} a_{12} \\ \vdots \\ a_{U2} \\ \end{array} \right) + \cdots + x_{V}\left( \begin{array}{c} a_{1V} \\ \vdots \\ a_{UV} \\ \end{array} \right) \] To indicate one common usage for matrix transformation in a data context, suppose we consider our data matrix $\mathbf{X} = \{x_{ij}\}_{N \times P}$, where $x_{ij}$ represents an observation for subject $i$ on variable $j$. We would like to use matrix transformations to produce a standardized matrix $\mathbf{Z} = \{(x_{ij} - \bar{x_{j}})/s_{j}\}_{N \times P}$, where $\bar{x_{j}}$ is the mean of the entries in the $j^{th}$ column and $s_{j}$ is the corresponding standard deviation; thus, the columns of $\mathbf{Z}$ all have mean zero and standard deviation one. A matrix expression for this transformation could be written as follows: \[ \mathbf{Z}_{N \times P} = (\mathbf{I}_{N \times N} - (\frac{1}{N})\mathbf{E}_{N \times N})\mathbf{X}_{N \times P}\mathbf{D}_{P \times P} \] \noindent where $\mathbf{I}$ is the identity matrix, $\mathbf{E}$ contains all ones, and $\mathbf{D}$ is a diagonal matrix containing $\frac{1}{s_{1}}, \frac{1}{s_{2}}, \ldots, \frac{1}{s_{P}}$, along the main diagonal positions. Thus, $(\mathbf{I}_{N \times N} - (\frac{1}{N})\mathbf{E}_{N \times N})\mathbf{X}_{N \times P}$ produces a matrix with columns deviated from the column means; a postmultiplication by $\mathbf{D}$ carries out the within column division by the standard deviations. Finally, if we define the expression $(\frac{1}{N})(\mathbf{Z}'\mathbf{Z})_{P \times P} \equiv \mathbf{R}_{P \times P}$, we have the familiar correlation coefficient matrix among the $P$ variables. \subsection{Matrix and Vector Orthogonality} Two vectors, $\mathbf{x}$ and $\mathbf{y}$, are said to be \emph{orthogonal} if $\mathbf{x}'\mathbf{y} = 0$, and would lie at right angles when graphed. If, in addition, $\mathbf{x}$ and $\mathbf{y}$ are both of unit length (i.e., $\sqrt{\mathbf{x}'\mathbf{x}} = \sqrt{\mathbf{y}'\mathbf{y}} = 1$), then they are said to be \emph{orthonormal}. A square matrix $\mathbf{T}_{U \times U}$ is said to be \emph{orthogonal} if its rows form a set of mutually orthonormal vectors. An example (called a Helmert matrix of order 3) follows: \[ \mathbf{T} = \left( \begin{array}{rrr} 1/\sqrt{3} & 1/\sqrt{3} & 1/\sqrt{3} \\ 1/\sqrt{2} & -1/\sqrt{2} & 0 \\ -1/\sqrt{6} & -1/\sqrt{6} & 2/\sqrt{6} \\ \end{array} \right) \] There are several nice properties of orthogonal matrices that we will see again in our various discussions to follow: \smallskip (A) $\mathbf{T}\mathbf{T}' = \mathbf{T}'\mathbf{T} = \mathbf{I}$; thus, the inverse of $\mathbf{T}$ is $\mathbf{T}'$; \smallskip (B) the columns of $\mathbf{T}$ are orthonormal; \smallskip (C) $\det(\mathbf{T}) = \pm 1$; \smallskip (D) if $\mathbf{T}$ and $\mathbf{R}$ are orthogonal, then so is $\mathbf{T}\mathbf{R}$; \smallskip (E) vectors lengths do not change under an orthogonal transformation: to see this, let $\mathbf{y} = \mathbf{T}\mathbf{x}$; then \[ \mathbf{y}'\mathbf{y} = (\mathbf{T}\mathbf{x})'(\mathbf{T}\mathbf{x}) = \mathbf{x}'\mathbf{T}'\mathbf{T}\mathbf{x} = \mathbf{x}'\mathbf{I}\mathbf{x} = \mathbf{x}'\mathbf{x} \] \subsection{Matrix Rank} An arbitrary matrix, $\mathbf{A}$, of order $U \times V$ can be written either in terms of its $U$ rows, say, $\mathbf{r}'_{1}, \mathbf{r}'_{2}, \ldots, \mathbf{r}'_{U}$ or its $V$ columns, $\mathbf{c}_{1}, \mathbf{c}_{2}, \ldots, \mathbf{c}_{V}$, where \[ \mathbf{r}'_{u} = \left( \begin{array}{ccc} a_{u1} & \cdots & a_{uV} \\ \end{array} \right); \; \mathbf{c}_{v} = \left( \begin{array}{c} a_{1v} \\ \vdots \\ a_{Uv} \\ \end{array} \right) \] \noindent and \[ \mathbf{A}_{U \times V} = \left( \begin{array}{c} \mathbf{r}'_{1} \\ \mathbf{r}'_{2} \\ \vdots \\ \mathbf{r}'_{U} \\ \end{array} \right) = \left( \begin{array}{cccc} \mathbf{c}_{1} & \mathbf{c}_{2} & \cdots & \mathbf{c}_{V} \\ \end{array} \right) \] \noindent The maximum number of linearly independent rows of $\mathbf{A}$ and the maximum number of linearly independent columns is the same; this common number is defined to be the \emph{rank} of $\mathbf{A}$. A matrix is said to be of \emph{full rank} is the rank is equal to the minimum of $U$ and $V$. \smallskip Matrix rank has a number of useful properties: \smallskip (A) $\mathbf{A}$ and $\mathbf{A}'$ have the same rank; \smallskip (B) $\mathbf{A}'\mathbf{A}$, $\mathbf{A}\mathbf{A}'$, and $\mathbf{A}$ have the same rank; \smallskip (C) the rank of a matrix is unchanged by a pre- or postmultiplication by a nonsingular matrix; \smallskip (D) the rank of a matrix is unchanged by what are called elementary row and column operations: (a) interchange of two rows or two columns; (2) multiplication or a row or a column by a scalar; (3) addition of a row (or column) to another row (or column). This is true because any elementary operation can be represented by a premultiplication (if the operation is to be on rows) or a postmultiplication (if the operation is to be on columns) of a nonsingular matrix. \smallskip To give a simple example, suppose we wish to perform some elementary row and column operations on the matrix \[ \left( \begin{array}{ccc} 1 & 1 & 1 \\ 1 & 0 & 2 \\ 3 & 2 & 4 \\ \end{array} \right) \] \noindent To interchange the first two rows of this latter matrix, interchange the first two rows of an identity matrix and premultiply; for the first two columns to be interchanged, carry out the operation on the identity and post-multiply: \[ \left( \begin{array}{ccc} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{ccc} 1 & 1 & 1 \\ 1 & 0 & 2 \\ 3 & 2 & 4 \\ \end{array} \right) = \left( \begin{array}{ccc} 1 & 0 & 2 \\ 1 & 1 & 1 \\ 3 & 2 & 4 \\ \end{array} \right) \] \[ \left( \begin{array}{ccc} 1 & 1 & 1 \\ 1 & 0 & 2 \\ 3 & 2 & 4 \\ \end{array} \right) \left( \begin{array}{ccc} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) = \left( \begin{array}{ccc} 1 & 1 & 1 \\ 0 & 1 & 2 \\ 2 & 3 & 4 \\ \end{array} \right) \] \noindent To multiply a row of our example matrix (e.g., the second row by $5$), multiply the desired row of an identity matrix and premultiply; for multiplying a specific column (e.g., the second column by $5$), carry out the operation of the identity and post-multiply: \[ \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 5 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{ccc} 1 & 1 & 1 \\ 1 & 0 & 2 \\ 3 & 2 & 4 \\ \end{array} \right) = \left( \begin{array}{ccc} 1 & 1 & 1 \\ 5 & 0 & 10 \\ 3 & 2 & 4 \\ \end{array} \right) \] \[ \left( \begin{array}{ccc} 1 & 1 & 1 \\ 1 & 0 & 2 \\ 3 & 2 & 4 \\ \end{array} \right) \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 5 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) = \left( \begin{array}{ccc} 1 & 5 & 1 \\ 1 & 0 & 2 \\ 3 & 10 & 4 \\ \end{array} \right) \] \noindent To add one row to a second (e.g., the first row to the second), carry out the operation on the identity and premultiply; to add one column to a second (e.g., the first column to the second), carry out the operation of the identity and post-multiply: \[ \left( \begin{array}{ccc} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{ccc} 1 & 1 & 1 \\ 1 & 0 & 2 \\ 3 & 2 & 4 \\ \end{array} \right) = \left( \begin{array}{ccc} 1 & 1 & 1 \\ 2 & 1 & 3 \\ 3 & 2 & 4 \\ \end{array} \right) \] \[ \left( \begin{array}{ccc} 1 & 1 & 1 \\ 1 & 0 & 2 \\ 3 & 2 & 4 \\ \end{array} \right) \left( \begin{array}{ccc} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) = \left( \begin{array}{ccc} 1 & 2 & 1 \\ 1 & 1 & 2 \\ 3 & 5 & 4 \\ \end{array} \right) \] In general, by performing elementary row and column operations, any $U \times V$ matrix can be reduced to a \emph{canonical form}: \[\left( \begin{array}{cccccc} 1 & \cdots & 0 & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 1 & 0 & \cdots & 0 \\ 0 & \cdots & 0 & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & 0 & \cdots & 0 \\ \end{array} \right) \] \noindent The rank of a matrix can then be found by counting the number of ones in the above matrix. \smallskip Given a $U \times V$ matrix, $\mathbf{A}$, there exist $s$ nonsingular elementary row operation matrices, $\mathbf{R}_{1},\ldots,\mathbf{R}_{s}$, and $t$ nonsingular elementary column operation matrices, $\mathbf{C}_{1},\ldots,\mathbf{C}_{t}$ such that $\mathbf{R}_{s} \cdots \mathbf{R}_{1}\mathbf{A}\mathbf{C}_{1} \cdots \mathbf{C}_{t}$ is in canonical form. Moreover, if $\mathbf{A}$ is square ($U \times U$) and of full rank (i.e., $\det(\mathbf{A}) \ne 0$), then there are $s$ nonsingular elementary row operation matrices, $\mathbf{R}_{1},\ldots,\mathbf{R}_{s}$, and $t$ nonsingular elementary column operation matrices, $\mathbf{C}_{1},\ldots,\mathbf{C}_{t}$, such that $\mathbf{R}_{s} \cdots \mathbf{R}_{1}\mathbf{A} = \mathbf{I}$ or $\mathbf{A}\mathbf{C}_{1} \cdots \mathbf{C}_{t} = \mathbf{I}$. Thus, $\mathbf{A}^{-1}$ can be found either as $ \mathbf{R}_{s} \cdots \mathbf{R}_{1}$ or as $\mathbf{C}_{1} \cdots \mathbf{C}_{t}$. In fact, a common way in which an inverse is calculated ``by hand'' starts with both $\mathbf{A}$ and $\mathbf{I}$ on the same sheet of paper; when reducing $\mathbf{A}$ step-by-step, the same operations are then applied to $\mathbf{I}$, building up the inverse until the canonical form is reached in the reduction of $\mathbf{A}$. \subsection{Using Matrices to Solve Equations} Suppose we have a set of $U$ equations in $V$ unknowns: \[ \begin{array}{ccccccc} a_{11}x_{1} & + & \cdots & + & a_{1V}x_{1} & = & c_{1} \\ \vdots & \vdots & \cdots & \vdots & \vdots & \vdots & \vdots \\ a_{U1}x_{1} & + & \cdots & + & a_{UV}x_{V} & = & c_{U} \\ \end{array} \] \noindent If we let \[ \mathbf{A} = \left( \begin{array}{ccc} a_{11} & \cdots & a_{1V} \\ \vdots & \ddots & \vdots \\ a_{U1} & \cdots & a_{UV} \\ \end{array} \right); \; \mathbf{x} = \left( \begin{array}{c} x_{1} \\ \vdots \\ x_{V} \\ \end{array} \right); \; \mathbf{c} = \left( \begin{array}{c} c_{1} \\ \vdots \\ c_{U} \\ \end{array} \right) \] \noindent then the equations can be written as follows: $\mathbf{A}_{U \times V}\mathbf{x}_{V \times 1} = \mathbf{c}_{U \times 1}$. In the simplest instance, $\mathbf{A}$ is square and nonsingular, implying that a solution may be given simply as $\mathbf{x} = \mathbf{A}^{-1}\mathbf{c}$. If there are fewer (say, $S < V$ linearly independent) equations than unknowns (so, $S$ is the rank of $\mathbf{A}$), then we can solve for $S$ unknowns in terms of the constants $c_{1},\ldots,c_{U}$ and the remaining $V - S$ unknowns. We will see how this works in our discussion of obtaining eigenvectors that correspond to certain eigenvalues in a section to follow. Generally, the set of equations is said to be \emph{consistent} if a solution exists, i.e., a linear combination of the column vectors of $\mathbf{A}$ can be used to define $\mathbf{c}$: \[ x_{1}\left( \begin{array}{c} a_{11} \\ \vdots \\ a_{U1} \\ \end{array} \right) + \cdots + x_{V}\left( \begin{array}{c} a_{1V} \\ \vdots \\ a_{UV} \\ \end{array} \right) = \left( \begin{array}{c} c_{1} \\ \vdots \\ c_{U} \\ \end{array} \right) \] \noindent or the augmented matrix $(\mathbf{A} \; \mathbf{c})$ has the same rank as $\mathbf{A}$; otherwise no solution exists and the system of equations is said to be \emph{inconsistent}. \subsection{Quadratic Forms} Suppose $\mathbf{A}_{U \times U}$ is symmetric and let $\mathbf{x}' = (x_{1}, \dots, x_{U})$. A \emph{quadratic form} is defined by \[ \mathbf{x}'\mathbf{A}\mathbf{x} = \sum_{u=1}^{U}\sum_{v=1}^{U} a_{uv}x_{u}x_{v} = \] \[ a_{11}x_{1}^{2} + a_{22}x_{2}^{2} + \cdots + a_{UU}x_{U}^{2} + 2a_{12}x_{1}x_{2} + \cdots + 2a_{1U}x_{1}x_{U} + \cdots + 2a_{(U-1)U}x_{U-1}x_{U} \] \noindent For example, $\sum_{u=1}^{U}(x_{u} - \bar{x})^{2}$, where $\bar{x}$ is the mean of the entries in $\mathbf{x}$, is a quadratic form because it can be written as \[ \left( \begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{U} \\ \end{array} \right)' \left( \begin{array}{cccc} (U-1)/U & -1/U & \cdots & -1/U \\ -1/U & (U-1)/U & \cdots & -1/U \\ \vdots & \vdots & \ddots & \vdots \\ -1/U & -1/U & \cdots & (U-1)/U \\ \end{array} \right) \left( \begin{array}{c} x_{1} \\ x_{2}\\ \vdots \\ x_{U} \\ \end{array} \right) \] \noindent Due to the ubiquity of sum-of-squares in statistics, it should be of no surprise that quadratic forms play a central role in multivariate analysis. \smallskip A symmetric matrix $\mathbf{A}$ (and associated quadratic form) are called \emph{positive definite} (p.d.) if $\mathbf{x}'\mathbf{A}\mathbf{x} > 0$ for all $\mathbf{x} \ne \mathbf{0}$ (the zero vector); if $\mathbf{x}'\mathbf{A}\mathbf{x} \ge 0$ for all $\mathbf{x}$, then $\mathbf{A}$ is \emph{positive semi-definite} (p.s.d). We could have negative definite, negative semi-definite, and indefinite forms as well. Note that a correlation or covariance matrix is at least positive semi-definite, and satisfies the stronger condition of being positive definite if the vectors of the variables on which the correlation or covariance matrix is based, are linearly independent. \subsection{Multiple Regression} One of the most common topics in any beginning statistics class is \emph{multiple regression} that we now formulate (in matrix terms) as the relation between a dependent random variable $Y$ and a collection of $K$ independent variables, $X_{1}, X_{2}, \ldots, X_{K}$. Suppose we have $N$ subjects on which we observe $Y$, and arrange these values into an $N \times 1$ vector: \[ \mathbf{Y} = \left( \begin{array}{c} Y_{1} \\ Y_{2} \\ \vdots \\ Y_{N}\\ \end{array} \right) \] \noindent The observations on the $K$ independent variables are also placed in vectors: \[ \mathbf{X}_{1} = \left( \begin{array}{c} X_{11} \\ X_{21} \\ \vdots \\ X_{N1} \\ \end{array} \right) ; \; \mathbf{X}_{2} = \left( \begin{array}{c} X_{12} \\ X_{22} \\ \vdots \\ X_{N2} \\ \end{array} \right) ; \; \ldots ; \; \mathbf{X}_{K} = \left( \begin{array}{c} X_{1K} \\ X_{2K} \\ \vdots \\ X_{NK} \\ \end{array} \right) \] \noindent It would be simple if the vector $\mathbf{Y}$ were linearly dependent on $\mathbf{X}_{1}, \mathbf{X}_{2}, \ldots, \mathbf{X}_{K}$ because then \[ \mathbf{Y} = b_{1}\mathbf{X}_{1} + b_{2}\mathbf{X}_{2} + \cdots + b_{K}\mathbf{X}_{K} \] \noindent for some values $b_{1}, \ldots, b_{K}$. We could always write for \emph{any} values of $b_{1}, \ldots, b_{K}$: \[ \mathbf{Y} = b_{1}\mathbf{X}_{1} + b_{2}\mathbf{X}_{2} + \cdots + b_{K}\mathbf{X}_{K} + \mathbf{e} \] \noindent where \[ \mathbf{e} = \left( \begin{array}{c} e_{1} \\ \vdots \\ e_{N} \\ \end{array} \right) \] \noindent is an error vector. To formulate our task as an optimization problem (least-squares), we wish to find a good set of weights, $b_{1}, \ldots, b_{K}$, so the length of $\mathbf{e}$ is minimized, i.e., $\mathbf{e}'\mathbf{e}$ is made as small as possible. \smallskip As notation, let \[ \mathbf{Y}_{N \times 1} = \mathbf{X}_{N \times K}\mathbf{b}_{K \times 1} + \mathbf{e}_{N \times 1} \] \noindent where \[ \mathbf{X} =\left( \begin{array}{ccc} \mathbf{X}_{1} & \ldots & \mathbf{X}_{K} \\ \end{array} \right); \; \mathbf{b} = \left( \begin{array}{c} b_{1} \\ \vdots \\ b_{K} \\ \end{array} \right) \] \noindent To minimize $\mathbf{e}'\mathbf{e} = (\mathbf{Y} - \mathbf{X}\mathbf{b})'(\mathbf{Y} - \mathbf{X}\mathbf{b})$, we use the vector $\mathbf{b}$ that satisfies what are called the normal equations: \[ \mathbf{X}'\mathbf{X}\mathbf{b} = \mathbf{X}'\mathbf{Y} \] \noindent If $\mathbf{X}'\mathbf{X}$ is nonsingular (i.e., $\det(\mathbf{X}'\mathbf{X}) \ne 0$; or equivalently, $\mathbf{X}_{1} , \ldots , \mathbf{X}_{K}$ are linearly independent), then \[ \mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y} \] \noindent The vector that is ``closest'' to $\mathbf{Y}$ in our least-squares sense, is $\mathbf{X}\mathbf{b}$; this is a linear combination of the columns of $\mathbf{X}$ (or in other jargon, $\mathbf{X}\mathbf{b}$ defines the \emph{projection} of $\mathbf{Y}$ into the space defined by (all linear combinations of) the columns of $\mathbf{X}$. \smallskip In statistical uses of multiple regression, the estimated variance-covariance matrix of the regression coefficients, $b_{1}, \ldots, b_{K}$, is given as $(\frac{1}{N-K})\mathbf{e}'\mathbf{e}(\mathbf{X}'\mathbf{X})^{-1}$, where $(\frac{1}{N-K})\mathbf{e}'\mathbf{e}$ is an (unbiased) estimate of the error variance for the distribution from which the errors are assumed drawn. Also, in multiple regression instances that usually involve an additive constant, the latter is obtained from a weight attached to an independent variable defined to be identically one. \smallskip In multivariate multiple regression where there are, say, $T$ dependent variables (each represented by an $N \times 1$ vector), the dependent vectors are merely concatenated together into an $N \times T$ matrix, $\mathbf{Y}_{N \times T}$; the solution to the normal equations now produces a matrix $\mathbf{B}_{K \times T} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ of regression coefficients. In effect, this general expression just uses each of the dependent variables separately and adjoins all the results. \section{Eigenvectors and Eigenvalues} Suppose we are given a square matrix, $\mathbf{A}_{U \times U}$, and consider the polynomial $\det(\mathbf{A} - \lambda\mathbf{I})$ in the unknown value $\lambda$, referred to as Laplace's expansion: \[ \det(\mathbf{A} - \lambda\mathbf{I}) = (-\lambda)^{U} + S_{1}(-\lambda)^{U-1} + \cdots + S_{U-1}(-\lambda)^{-1} + S_{U}(-\lambda)^{0} \] \noindent where $S_{u}$ is the sum of all $u \times u$ principal minor determinants. A \emph{principal} minor determinant is obtained from a submatrix formed from $\mathbf{A}$ that has $u$ diagonal elements left in it. Thus, $S_{1}$ is the trace of $\mathbf{A}$ and $S_{U}$ is the determinant. \smallskip There are $U$ roots, $\lambda_{1}, \ldots, \lambda_{U}$, of the equation $\det(\mathbf{A} - \lambda\mathbf{I}) = 0$, given that the left-hand-side is a $U^{th}$ degree polynomial. The roots are called the \emph{eigenvalues} of $\mathbf{A}$. There are a number of properties of eigenvalues that prove generally useful: \smallskip (A) $\det{\mathbf{A}} = \prod_{u=1}^{U} \lambda_{u}$; trace($\mathbf{A}) = \sum_{u=1}^{U} \lambda_{u}$; \smallskip (B) if $\mathbf{A}$ is symmetric with real elements, then all $\lambda_{u}$ are real; \smallskip (C) if $\mathbf{A}$ is positive definite, then all $\lambda_{u}$ are positive (strictly greater than zero); if $\mathbf{A}$ is positive semi-definite, then all $\lambda_{u}$ are nonnegative (greater than or equal to zero); \smallskip (D) if $\mathbf{A}$ is symmetric and positive semi-definite with rank $R$, then there are $R$ positive roots and $U-R$ zero roots; \smallskip (E) the nonzero roots of $\mathbf{A}\mathbf{B}$ are equal to those of $\mathbf{B}\mathbf{A}$; thus, the trace of $\mathbf{A}\mathbf{B}$ is equal to the trace of $\mathbf{B}\mathbf{A}$; \smallskip (F) eigenvalues of a diagonal matrix are the diagonal elements themselves; \smallskip (G) for any $U \times V$ matrix $\mathbf{B}$, the ranks of $\mathbf{B}$, $\mathbf{B}'\mathbf{B}$, and $\mathbf{B}\mathbf{B}'$ are all the same. Thus, because $\mathbf{B}'\mathbf{B}$ (and $\mathbf{B}\mathbf{B}'$) are symmetric and positive semi-definite (i.e., $\mathbf{x}'(\mathbf{B}'\mathbf{B})\mathbf{x} \ge 0$ because $(\mathbf{B}\mathbf{x})'(\mathbf{B}\mathbf{x})$ is a sum-of-squares which is always nonnegative), we can use (D) to find the rank of $\mathbf{B}$ by counting the positive roots of $\mathbf{B}'\mathbf{B}$. \smallskip We carry through a small example below: \[ \mathbf{A} = \left( \begin{array}{ccc} 7 & 0 & 1 \\ 0 & 7 & 2 \\ 1 & 2 & 3 \\ \end{array} \right) \] \[ S_{1} = \mathrm{trace}(\mathbf{A}) = 17 \] \[ S_{2} = \det(\left( \begin{array}{cc} 7 & 0 \\ 0 & 7 \\ \end{array} \right)) + \det(\left( \begin{array}{cc} 7 & 1 \\ 1 & 3 \\ \end{array} \right)) + \det(\left( \begin{array}{cc} 7 & 2 \\ 2 & 3 \\ \end{array} \right)) = 49 + 20 + 17 = 86 \] \[ S_{3} = \det(\mathbf{A}) = 147 + 0 + 0 - 7 - 28 - 0 = 112 \] \noindent Thus, \[ \det(\mathbf{A} - \lambda\mathbf{I}) = (-\lambda)^{3} + 17(-\lambda)^{2} + 86(-\lambda)^{1} + 112 = \] \[ -\lambda^{3} + 17\lambda^{2} - 86\lambda + 112 = -(\lambda -2)(\lambda - 8)(\lambda-7) = 0 \] \noindent which gives roots of 2, 8, and 7. \smallskip If $\lambda_{u}$ is an eigenvalue of $\mathbf{A}$, then the equations $[\mathbf{A} - \lambda_{u}\mathbf{I}]\mathbf{x}_{u} = 0$ have a nontrivial solution (i.e., the determinant of $\mathbf{A} - \lambda_{u}\mathbf{I}$ vanishes, and so the inverse of $\mathbf{A} - \lambda_{u}\mathbf{I}$ does not exist). The solution is called an \emph{eigenvector} (associated with the corresponding eigenvalue), and can be characterized by the following condition: \[ \mathbf{A}\mathbf{x}_{u} = \lambda_{u}\mathbf{x}_{u} \] \noindent An eigenvector is determined up to a scale factor only, so typically we normalize to unit length (which then gives a $\pm$ option for the two possible unit length solutions). \smallskip We continue our simple example and find the corresponding eigenvalues: when $\lambda = 2$, we have the equations (for $[\mathbf{A} - \lambda\mathbf{I}]\mathbf{x} = \mathbf{0}$) \[ \left( \begin{array}{ccc} 5 & 0 & 1 \\ 0 & 5 & 2 \\ 1 & 2 & 1 \\ \end{array} \right) \left( \begin{array}{c} x_{1} \\ x_{2} \\ x_{3} \\ \end{array} \right) = \left(% \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array}% \right) \] \noindent with an arbitrary solution of \[ \left(% \begin{array}{r} -\frac{1}{5}a \\ -\frac{2}{5}a \\ a \\ \end{array} \right) \] \noindent Choosing $a$ to be $+\frac{5}{\sqrt{30}}$ to obtain one of the two possible normalized solutions, we have as our final eigenvector for $\lambda = 2$: \[ \left( \begin{array}{r} -\frac{1}{\sqrt{30}} \\ -\frac{2}{\sqrt{30}} \\ \frac{5}{\sqrt{30}} \\ \end{array} \right) \] \noindent For $\lambda = 7$ we will use the normalized eigenvector of \[\left( \begin{array}{r} -\frac{2}{\sqrt{5}} \\ \frac{1}{\sqrt{5}} \\ 0 \\ \end{array} \right) \] \noindent and for $\lambda = 8$, \[ \left( \begin{array}{r} \frac{1}{\sqrt{6}} \\ \frac{2}{\sqrt{6}} \\ \frac{1}{\sqrt{6}} \\ \end{array} \right) \] One of the interesting properties of eigenvalues/eigenvectors for a symmetric matrix $\mathbf{A}$ is that if $\lambda_{u}$ and $\lambda_{v}$ are distinct eigenvalues, then the corresponding eigenvectors, $\mathbf{x}_{u}$ and $\mathbf{x}_{v}$, are orthogonal (i.e., $\mathbf{x}_{u}'\mathbf{x}_{v} = 0$). We can show this in the following way: the defining conditions of \[ \mathbf{A}\mathbf{x}_{u} = \lambda_{u}\mathbf{x}_{u} \] \[ \mathbf{A}\mathbf{x}_{v} = \lambda_{v}\mathbf{x}_{v} \] \noindent lead to \[ \mathbf{x}_{v}'\mathbf{A}\mathbf{x}_{u} = \mathbf{x}_{v}'\lambda_{u}\mathbf{x}_{u} \] \[ \mathbf{x}_{u}'\mathbf{A}\mathbf{x}_{v} = \mathbf{x}_{u}'\lambda_{v}\mathbf{x}_{v} \] \noindent Because $\mathbf{A}$ is symmetric and the left-hand-sides of these two expressions are equal (they are one-by-one matrices and equal to their own transposes), the right-hand-sides must also be equal. Thus, \[ \mathbf{x}_{v}'\lambda_{u}\mathbf{x}_{u} = \mathbf{x}_{u}'\lambda_{v}\mathbf{x}_{v} \Rightarrow \] \[ \mathbf{x}_{v}'\mathbf{x}_{u}\lambda_{u} = \mathbf{x}_{u}'\mathbf{x}_{v}\lambda_{v} \] \noindent Due to the equality of $\mathbf{x}_{v}'\mathbf{x}_{u}$ and $\mathbf{x}_{u}'\mathbf{x}_{v}$, and by assumption, $\lambda_{u} \ne \lambda_{v}$, the inner product $\mathbf{x}_{v}'\mathbf{x}_{u}$ must be zero for the last displayed equality to hold. \smallskip In summary of the above discussion, for every real symmetric matrix $\mathbf{A}_{U \times U}$, there exists an orthogonal matrix $\mathbf{P}$ (i.e., $\mathbf{P}'\mathbf{P} = \mathbf{P}\mathbf{P}' = \mathbf{I}$) such that $\mathbf{P}'\mathbf{A}\mathbf{P} = \mathbf{D}$, where $\mathbf{D}$ is a diagonal matrix containing the eigenvalues of $\mathbf{A}$, and \[ \mathbf{P} = \left( \begin{array}{ccc} \mathbf{p}_{1} & \ldots & \mathbf{p}_{U} \\ \end{array} \right) \] \noindent where $\mathbf{p}_{u}$ is a normalized eigenvector associated with $\lambda_{u}$ for $1 \le u \le U$. If the eigenvalues are not distinct, it is still possible to choose the eigenvectors to be orthogonal. Finally, because $\mathbf{P}$ is an orthogonal matrix (and $\mathbf{P}'\mathbf{A}\mathbf{P} = \mathbf{D} \Rightarrow \mathbf{P}\mathbf{P}'\mathbf{A}\mathbf{}\mathbf{P}\mathbf{P}' = \mathbf{P}\mathbf{D}\mathbf{P}'$), we can finally represent $\mathbf{A}$ as \[ \mathbf{A} = \mathbf{P}\mathbf{D}\mathbf{P}' \] In terms of the small numerical example being used, we have for $\mathbf{P}'\mathbf{A}\mathbf{P} = \mathbf{D}$: \[ \left( \begin{array}{rrr} -\frac{1}{\sqrt{30}} & -\frac{2}{\sqrt{30}} & \frac{5}{\sqrt{30}} \\ -\frac{2}{\sqrt{5}} & \frac{1}{\sqrt{5}} & 0 \\ \frac{1}{\sqrt{6}} & \frac{2}{\sqrt{6}} & \frac{1}{\sqrt{6}} \\ \end{array} \right) \left( \begin{array}{rrr} 7 & 0 & 1 \\ 0 & 7 & 2 \\ 1 & 2 & 3 \\ \end{array} \right) \left( \begin{array}{rrr} -\frac{1}{\sqrt{30}} & -\frac{2}{\sqrt{5}} & \frac{1}{\sqrt{6}} \\ -\frac{2}{\sqrt{30}} & \frac{1}{\sqrt{5}} & \frac{2}{\sqrt{6}} \\ \frac{5}{\sqrt{30}} & 0 & \frac{1}{\sqrt{6}} \\ \end{array} \right) = \] \[ \left( \begin{array}{rrr} 2 & 0 & 0 \\ 0 & 7 & 0 \\ 0 & 0 & 8 \\ \end{array} \right) \] \noindent and for $\mathbf{P}\mathbf{D}\mathbf{P}' = \mathbf{A}$: \[ \left( \begin{array}{rrr} -\frac{1}{\sqrt{30}} & -\frac{2}{\sqrt{5}} & \frac{1}{\sqrt{6}} \\ -\frac{2}{\sqrt{30}} & \frac{1}{\sqrt{5}} & \frac{2}{\sqrt{6}} \\ \frac{5}{\sqrt{30}} & 0 & \frac{1}{\sqrt{6}} \\ \end{array} \right) \left( \begin{array}{rrr} 2 & 0 & 0 \\ 0 & 7 & 0 \\ 0 & 0 & 8 \\ \end{array} \right) \left( \begin{array}{rrr} -\frac{1}{\sqrt{30}} & -\frac{2}{\sqrt{30}} & \frac{5}{\sqrt{30}} \\ -\frac{2}{\sqrt{5}} & \frac{1}{\sqrt{5}} & 0 \\ \frac{1}{\sqrt{6}} & \frac{2}{\sqrt{6}} & \frac{1}{\sqrt{6}} \\ \end{array} \right) = \] \[ \left( \begin{array}{rrr} 7 & 0 & 1 \\ 0 & 7 & 2 \\ 1 & 2 & 3 \\ \end{array} \right) \] The representation of $\mathbf{A}$ as $\mathbf{P}\mathbf{D}\mathbf{P}'$ leads to several rather nice computational ``tricks.'' First, if $\mathbf{A}$ is p.s.d., we can define \[ \mathbf{D}^{1/2} \equiv \left( \begin{array}{ccc} \sqrt{\lambda_{1}} & \ldots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \ldots & \sqrt{\lambda_{U}} \\ \end{array} \right) \] \noindent and represent $\mathbf{A}$ as \[ \mathbf{A} = \mathbf{P}\mathbf{D}^{1/2}\mathbf{D}^{1/2}\mathbf{P}' = \mathbf{P}\mathbf{D}^{1/2}(\mathbf{P}\mathbf{D}^{1/2})' = \mathbf{L}\mathbf{L}', \mathrm{say}. \] \noindent In other words, we have ``factored'' $\mathbf{A}$ into $\mathbf{L}\mathbf{L}'$, for \[ \mathbf{L} = \mathbf{P}\mathbf{D}^{1/2} = \left( \begin{array}{cccc} \sqrt{\lambda_{1}}\mathbf{p}_{1} & \sqrt{\lambda_{2}}\mathbf{p}_{2} & \ldots & \sqrt{\lambda_{U}}\mathbf{p}_{U} \\ \end{array} \right) \] \noindent Secondly, if $\mathbf{A}$ is p.d., we can define \[ \mathbf{D}^{-1} \equiv \left( \begin{array}{ccc} \frac{1}{\lambda_{1}} & \ldots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \ldots & \frac{1}{\lambda_{U}} \\ \end{array} \right) \] \noindent and represent $\mathbf{A}^{-1}$ as \[ \mathbf{A}^{-1} = \mathbf{P}\mathbf{D}^{-1}\mathbf{P}' \] \noindent To verify, \[ \mathbf{A}\mathbf{A}^{-1} = (\mathbf{P}\mathbf{D}\mathbf{P}')(\mathbf{P}\mathbf{D}^{-1}\mathbf{P}') = \mathbf{I} \] \noindent Thirdly, to define a ``square root'' matrix, let $\mathbf{A}^{1/2} \equiv \mathbf{P}\mathbf{D}^{1/2}\mathbf{P}'$. To verify, $\mathbf{A}^{1/2}\mathbf{A}^{1/2} = \mathbf{P}\mathbf{D}\mathbf{P}' = \mathbf{A}$. \smallskip There is a generally interesting way to represent the multiplication of two matrices considered as collections of column and row vectors, respectively, where the final answer is a sum of outer products of vectors. This view will prove particularly useful in our discussion of principal component analysis. Suppose we have two matrices $\mathbf{B}_{U \times V}$, represented as a collection of its $V$ columns: \[ \mathbf{B} = \left( \begin{array}{cccc} \mathbf{b}_{1} & \mathbf{b}_{2} & \ldots & \mathbf{b}_{V} \\ \end{array} \right) \] \noindent and $\mathbf{C}_{V \times W}$, represented as a collection of its $V$ rows: \[ \mathbf{C} = \left( \begin{array}{c} \mathbf{c}_{1}' \\ \mathbf{c}_{2}' \\ \vdots \\ \mathbf{c}_{V}' \\ \end{array} \right) \] \noindent The product $\mathbf{B}\mathbf{C} = \mathbf{D}$ can be written as \[ \mathbf{B}\mathbf{C} = \left( \begin{array}{cccc} \mathbf{b}_{1} & \mathbf{b}_{2} & \ldots & \mathbf{b}_{V} \\ \end{array} \right) \left( \begin{array}{c} \mathbf{c}_{1}' \\ \mathbf{c}_{2}' \\ \vdots \\ \mathbf{c}_{V}' \\ \end{array} \right) = \] \[ \mathbf{b}_{1}\mathbf{c}_{1}' + \mathbf{b}_{2}\mathbf{c}_{2}' + \cdots + \mathbf{b}_{V}\mathbf{c}_{V}' = \mathbf{D} \] As an example, consider the \emph{spectral decomposition} of $\mathbf{A}$ considered above as $\mathbf{P}\mathbf{D}\mathbf{P}'$, and where from now on, without loss of any generality, the diagonal entries in $\mathbf{D}$ are ordered as $\lambda_{1} \ge \lambda_{2} \ge \cdots \ge \lambda_{U}$. We can represent $\mathbf{A}$ as \[ \mathbf{A}_{U \times U} = \left( \begin{array}{ccc} \sqrt{\lambda_{1}}\mathbf{p}_{1} & \ldots & \sqrt{\lambda_{U}}\mathbf{p}_{U} \\ \end{array} \right) \left( \begin{array}{c} \sqrt{\lambda_{1}}\mathbf{p}_{1}' \\ \vdots \\ \sqrt{\lambda_{U}}\mathbf{p}_{U}' \\ \end{array} \right) = \] \[ \lambda_{1}\mathbf{p}_{1}\mathbf{p}_{1}' + \cdots + \lambda_{U}\mathbf{p}_{U}\mathbf{p}_{U}' \] \noindent If $\mathbf{A}$ is p.s.d. and of rank $R$, then the above sum obviously stops at $R$ components. In general, the matrix $\mathbf{B}_{U \times U}$ that is a rank $K$ ($ \le R$) least-squares approximation to $\mathbf{A}$ can be given by \[ \mathbf{B} = \lambda_{1}\mathbf{p}_{1}\mathbf{p}_{1}' + \cdots + \lambda_{k}\mathbf{p}_{K}\mathbf{p}_{K}' \] \noindent and the value of the loss function: \[ \sum_{v=1}^{U} \sum_{u=1}^{U} (a_{uv} - b_{uv})^{2} = \lambda_{K+1}^{2} + \cdots + \lambda_{U}^{2} \] \section{The Singular Value Decomposition of a Matrix} The \emph{singular value decomposition} (SVD) or the \emph{basic structure} of a matrix refers to the representation of \emph{any} rectangular $U \times V$ matrix, say, $\mathbf{A}$, as a triple product: \[ \mathbf{A}_{U \times V} = \mathbf{P}_{U \times R} \mathbf{\Delta}_{R \times R} \mathbf{Q}_{R \times V}' \] \noindent where the $R$ columns of $\mathbf{P}$ are orthonormal; the $R$ rows of $\mathbf{Q}'$ are orthonormal; $\mathbf{\Delta}$ is diagonal with ordered positive entries, $\delta_{1} \ge \delta_{2} \ge \cdots \ge \delta_{R} > 0$; and $R$ is the rank of $\mathbf{A}$. Or, alternatively, we can ``fill up'' this decomposition as \[ \mathbf{A}_{U \times V} = \mathbf{P}_{U \times U}^{*} \mathbf{\Delta}_{U \times V}^{*} \mathbf{Q}_{V \times V}^{*'} \] \noindent where the columns of $\mathbf{P}^{*}$ and rows of $\mathbf{Q}^{*'}$ are still orthonormal, and the diagonal matrix $\mathbf{\Delta}$ forms the upper-left-corner of $\mathbf{\Delta}^{*}$: \[ \mathbf{\Delta}^{*} = \left( \begin{array}{cc} \mathbf{\Delta} & \mathbf{\emptyset} \\ \mathbf{\emptyset} & \mathbf{\emptyset} \\ \end{array} \right) \] \noindent here, $\mathbf{\emptyset}$ represents an appropriately dimensioned matrix of all zeros. In analogy to the least-squares result of the last section, if a rank $K$ ($\le R$) matrix approximation to $\mathbf{A}$ is desired, say $\mathbf{B}_{U \times V}$, the first $K$ ordered entries in $\mathbf{\Delta}$ are taken: \[ \mathbf{B} = \delta_{1}\mathbf{p}_{1}\mathbf{q}_{1}' + \cdots + \delta_{K}\mathbf{p}_{K}\mathbf{q}_{K}' \] \noindent and the value of the loss function: \[ \sum_{v=1}^{V} \sum_{u=1}^{U} (a_{uv} - b_{uv})^{2} = \delta_{K+1}^{2} + \cdots + \delta_{R}^{2} \] \noindent This latter result of approximating one matrix (least-squares) by another of lower rank, is referred to as the Eckart-Young theorem in the psychometric literature. \smallskip Once one has the SVD of a matrix, a lot of representation needs can be expressed in terms of it. For example, suppose $\mathbf{A} = \mathbf{P}\mathbf{\Delta}\mathbf{Q}'$; the spectral decomposition of $\mathbf{A}\mathbf{A}'$ can then be given as \[ (\mathbf{P}\mathbf{\Delta}\mathbf{Q}')( \mathbf{P}\mathbf{\Delta}\mathbf{Q}')' = \mathbf{P}\mathbf{\Delta}\mathbf{Q}'\mathbf{Q}\mathbf{\Delta}\mathbf{P}' = \mathbf{P}\mathbf{\Delta}\mathbf{\Delta}\mathbf{P}' = \mathbf{P}\mathbf{\Delta}^{2}\mathbf{P}' \] \noindent Similarly, the spectral decomposition of $\mathbf{A}'\mathbf{A}$ is expressible as $\mathbf{Q}\mathbf{\Delta}^{2}\mathbf{Q}'$. \section{Common Multivariate Methods in Matrix Terms} In this section we give a very brief overview of some common methods of multivariate analysis in terms of the matrix ideas we have introduced thus far in this chapter. Later chapters (if they ever get written) will come back to these topics and develop them in more detail. \subsection{Principal Components} Suppose we have a data matrix $\mathbf{X}_{N \times P}$ = $\{x_{ij}\}$, with $x_{ij}$ referring as usual to the observation for subject $i$ on variable or column $j$: \[ \mathbf{X}_{N \times P} = \left( \begin{array}{cccc} x_{11} & x_{12} & \cdots & x_{1P} \\ x_{21} & x_{22} & \cdots & x_{2P} \\ \vdots & \vdots & \ddots & \vdots \\ x_{N1} & x_{N2} & \cdots & x_{NP} \\ \end{array} \right) \] \noindent The columns can be viewed as containing $N$ observations on each of $P$ random variables that we denote generically by $X_{1}, X_{2}, \ldots, X_{P}$. We let $\mathbf{A}$ denote the $P \times P$ sample covariance matrix obtained among the variables from $\mathbf{X}$, and let $\lambda_{1} \ge \cdots \ge \lambda_{P} \ge 0$ be its $P$ eigenvalues and $\mathbf{p}_{1}, \ldots, \mathbf{p}_{P}$ the corresponding normalized eigenvectors. Then, the linear combination \[ \mathbf{p}_{k}'\left( \begin{array}{c} X_{1} \\ \vdots \\ X_{P} \\ \end{array} \right) \] \noindent is called the $k^{th}$ (sample) \emph{principal component}. \smallskip There are (at least) two interesting properties of principal components to bring up at this time: \smallskip A) The $k^{th}$ principal component has maximum variance among all linear combinations defined by unit length vectors orthogonal to $\mathbf{p}_{1}, \ldots, \mathbf{p}_{k-1}$; also, it is uncorrelated with the components up to $k-1$; \smallskip B) $\mathbf{A} \approx \lambda_{1}\mathbf{p}_{1}\mathbf{p}_{1}' + \cdots + \lambda_{K}\mathbf{p}_{K}\mathbf{p}_{K}'$ gives a least-squares rank $K$ approximation to $\mathbf{A}$ (a special case of the Eckart-Young theorem for an arbitrary symmetric matrix). \subsection{Discriminant Analysis} Suppose we have a one-way analysis-of-variance (ANOVA) layout with $J$ groups ($n_{j}$ subjects in group $j$, $1 \le j \le J$), and $P$ measurements on each subject. If $x_{ijk}$ denotes person $i$, in group $j$, and the observation of variable $k$ ($1 \le i \le n_{j}$; $1 \le j \le J$; $1 \le k \le P$), then define the Between-Sum-of-Squares matrix \[\mathbf{B}_{P \times P} = \{\sum_{j=1}^{J} n_{j}(\bar{x}_{\cdot jk} - \bar{x}_{\cdot \cdot k})(\bar{x}_{\cdot jk'} - \bar{x}_{\cdot \cdot k'})\}_{P \times P} \] \noindent and the Within-Sum-of-Squares matrix \[\mathbf{W}_{P \times P} = \{ \sum_{j=1}^{J} \sum_{i=1}^{n_{j}} (x_{ijk} - \bar{x}_{\cdot jk})(x_{ijk'} - \bar{x}_{\cdot jk'}) \}_{P \times P} \] For the matrix product $\mathbf{W}^{-1}\mathbf{B}$, let $\lambda_{1}, \ldots, \lambda_{T} \ge 0$ be the eigenvectors ($T = \min(P, J-1)$, and $\mathbf{p}_{1}, \ldots, \mathbf{p}_{T}$ the corresponding normalized eigenvectors. Then, the linear combination \[ \mathbf{p}_{k}'\left( \begin{array}{c} X_{1} \\ \vdots \\ X_{P} \\ \end{array} \right) \] \noindent is called the $k^{th}$ \emph{discriminant function}. It has the valuable property of maximizing the univariate $F$-ratio subject to being uncorrelated with the earlier linear combinations. A variety of applications of discriminant functions exists in classification that we will come back to later. Also, standard multivariate ANOVA significance testing is based on various functions of the eigenvalues $\lambda_{1}, \ldots, \lambda_{T}$ and their derived sampling distributions. \subsection{Canonical Correlation} Suppose the collection of $P$ random variables that we have observed over the $N$ subjects is actually in the form of two ``batteries,'' $X_{1}, \ldots, X_{Q}$ and $X_{Q+1}, \ldots, X_{P}$, and the observed covariance matrix $\mathbf{A}_{P \times P}$ is partitioned into four parts: \[ \mathbf{A}_{P \times P} = \left( \begin{array}{cc} \mathbf{A}_{11} & \mathbf{A}_{12} \\ \mathbf{A}_{12}' & \mathbf{A}_{22} \\ \end{array} \right) \] \noindent where $\mathbf{A}_{11}$ is $Q \times Q$ and represents the observed covariances among the variables in the first battery; $\mathbf{A}_{22}$ is $(P-Q) \times (P-Q)$ and represents the observed covariances among the variables in the second battery; $\mathbf{A}_{12}$ is $Q \times (P-Q)$ and represents the observed covariances between the variables in the first and second batteries. Consider the following two equations in unknown vectors $\mathbf{a}$ and $\mathbf{b}$, and unknown scalar $\lambda$: \[ \mathbf{A}_{11}^{-1} \mathbf{A}_{12} \mathbf{A}_{22}^{-1} \mathbf{A}_{12}' \mathbf{a} = \lambda \mathbf{a} \] \[ \mathbf{A}_{22}^{-1} \mathbf{A}_{12}' \mathbf{A}_{11}^{-1} \mathbf{A}_{12} \mathbf{b} = \lambda \mathbf{b} \] \noindent There are $T$ solutions to these expressions (for $T = \min(Q,(P-Q))$), given by normalized unit-length vectors, $\mathbf{a}_{1}, \ldots, \mathbf{a}_{T}$ and $\mathbf{b}_{1}, \ldots, \mathbf{b}_{T}$; and a set of common $\lambda_{1} \ge \cdots \ge \lambda_{T} \ge 0$. \smallskip The linear combinations of the first and second batteries defined by $\mathbf{a}_{k}$ and $\mathbf{b}_{k}$ are the $k^{th}$ \emph{canonical variates} and have squared correlation of $\lambda_{k}$; they are uncorrelated with all other canonical variates (defined either in the first or second batteries). Thus, $\mathbf{a}_{1}$ and $\mathbf{b}_{1}$ are the first canonical variates with squared correlation of $\lambda_{1}$; among all linear combinations defined by unit-length vectors for the variables in the two batteries, this squared correlation is the highest it can be. (We note that the coefficient matrices $\mathbf{A}_{11}^{-1} \mathbf{A}_{12} \mathbf{A}_{22}^{-1} \mathbf{A}_{12}'$ and $\mathbf{A}_{22}^{-1} \mathbf{A}_{12}' \mathbf{A}_{11}^{-1} \mathbf{A}_{12}$ are not symmetric; thus, special symmetrizing and equivalent equation systems are typically used to obtain the solutions to the original set of expressions.) \subsection{Algebraic Restrictions on Correlations} A matrix $\mathbf{A}_{P \times P}$ that represents a covariance matrix among a collection of random variables, $X_{1}, \ldots, X_{P}$ is p.s.d.; and conversely, any p.s.d. matrix represents the covariance matrix for some collection of random variables. We partition $\mathbf{A}$ to isolate its last row and column as \[ \mathbf{A} = \left( \begin{array}{cc} \mathbf{B}_{(P-1) \times (P-1)} & \mathbf{g}_{(P-1) \times 1} \\ \mathbf{g}' & a_{PP} \\ \end{array} \right) \] \noindent $\mathbf{B}$ is the $(P-1) \times (P-1)$ covariance matrix among the variables $X_{1}, \ldots, X_{P-1}$; $\mathbf{g}$ is $(P-1) \times 1$ and contains the cross-covariance between the the first $P-1$ variables and the $P^{th}$; $a_{PP}$ is the variance for the $P^{th}$ variable. \smallskip Based on the observation that determinants of p.s.d. matrices are nonnegative, and a result on expressing determinants for partitioned matrices (that we do not give here), it must be true that \[ \mathbf{g}'\mathbf{B}^{-1}\mathbf{g} \le a_{PP} \] \noindent or if we think correlations rather than merely covariances (so the main diagonal of $\mathbf{A}$ consists of all ones): \[ \mathbf{g}'\mathbf{B}^{-1}\mathbf{g} \le 1 \] \noindent Given the correlation matrix $\mathbf{B}$, the possible values the correlations in $\mathbf{g}$ could have are in or on the ellipsoid defined in $P-1$ dimensions by $\mathbf{g}'\mathbf{B}^{-1}\mathbf{g} \le 1$. The important point is that we do not have a ``box'' in $P-1$ dimensions containing the correlations with sides extending the whole range of $\pm 1$; instead, some restrictions are placed on the observable correlations that gets defined by the size of the correlation in $\mathbf{B}$. For example, when $P = 3$, a correlation between variables $X_{1}$ and $X_{2}$ of $r_{12} = 0$ gives the ``degenerate'' ellipse of a circle for constraining the correlation values between $X_{1}$ and $X_{2}$ and the third variable $X_{3}$ (in a two-dimensional $r_{13}$ versus $r_{23}$ coordinate system); for $r_{12} = 1$, the ellipse flattens to a line in this same two-dimensional space. \smallskip Another algebraic restriction that can be seen immediately is based on the formula for the partial correlation between two variables, ``holding the third constant'': \[ \frac{r_{12} - r_{13}r_{23}}{\sqrt{(1-r_{13}^{2})(1-r_{23}^{2})}} \] \noindent Bounding the above by $\pm 1$ (because it is a correlation) and ``solving'' for $r_{12}$, gives the algebraic upper and lower bounds of \[ r_{12} \le r_{13}r_{23} + \sqrt{(1-r_{13}^{2})(1-r_{23}^{2})} \] \[ r_{13}r_{23} - \sqrt{(1-r_{13}^{2})(1-r_{23}^{2})} \le r_{12} \] \subsection{The Biplot} Let $\mathbf{A} = \{a_{ij}\}$ be an $n \times m$ matrix of rank $r$. We wish to find a second matrix $\mathbf{B} = \{b_{ij}\}$ of the same size, $n \times m$, but of rank $t$, where $t \le r$, such that the least squares criterion, $\sum_{i,j} (a_{ij} - b_{ij})^{2}$, is as small as possible overall all matrices of rank $t$. \smallskip The solution is to first find the singular value decomposition of $\mathbf{A}$ as $\mathbf{UDV}'$, where $\mathbf{U}$ is $n \times r$ and has orthonormal columns, $\mathbf{V}$ is $m \times r$ and has orthonormal columns, and $\mathbf{D}$ is $r \times r$, diagonal, with positive values $d_{1} \ge d_{2} \ge \cdots \ge d_{r} > 0$ along the main diagonal. Then, $\mathbf{B}$ is defined as $\mathbf{U^*D^*V^*}'$, where we take the first $t$ columns of $\mathbf{U}$ and $\mathbf{V}$ to obtain $\mathbf{U^*}$ and $\mathbf{V^*}$, respectively, and the first $t$ values, $d_{1} \ge \cdots \ge d_{t}$, to form a diagonal matrix $\mathbf{D^*}$. \smallskip The approximation of $\mathbf{A}$ by a rank $t$ matrix $\mathbf{B}$, has been one mechanism for representing the row and column objects defining $\mathbf{A}$ in a low-dimensional space of dimension $t$ through what can be generically labeled as a biplot (the prefix ``bi'' refers to the representation of both the row and column objects together in the same space). Explicitly, the approximation of $\mathbf{A}$ and $\mathbf{B}$ can be written as \[ \mathbf{B} = \mathbf{U^*D^*V^*}' = \mathbf{U^*D^{*\alpha}D^{(1-\alpha)}V^*}' = \mathbf{PQ}' \ ,\] \noindent where $\alpha$ is some chosen number between 0 and 1, $\mathbf{P} = \mathbf{U^*D^{*\alpha}}$ and is $n \times t$, $\mathbf{Q} = (\mathbf{D^{(1-\alpha)}V^*}')'$ and is $m \times t$. \smallskip The entries in $\mathbf{P}$ and $\mathbf{Q}$ define coordinates for the row and column objects in a $t$-dimensional space that, irrespective of the value of $\alpha$ chosen, have the following characteristic: \smallskip If a vector is drawn from the origin through the $i^{th}$ row point and the $m$ column points are projected onto this vector, the collection of such projections is proportional to the $i^{th}$ row of the approximating matrix $\mathbf{B}$. The same is true for projections of row points onto vectors from the origin through each of the column points. \subsection{The Procrustes Problem} Procrustes (the subduer), son of Poseidon, kept an inn benefiting from what he claimed to be a wonderful all-fitting bed. He lopped off excessive limbage from tall guests and either flattened short guests by hammering or stretched them by racking. The victim fitted the bed perfectly but, regrettably, died. To exclude the embarrassment of an initially exact-fitting guest, variants of the legend allow Procrustes two, different-sized beds. Ultimately, in a crackdown on robbers and monsters, the young Theseus fitted Procrustes to his own bed. (Gower and Dijksterhuis, 2004) \smallskip Suppose we have two matrices, $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$, each considered (for convenience) to be of the same size, $n \times p$. If you wish, $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ can be interpreted as two separate $p$-dimensional coordinate sets for the same set of $n$ objects. Our task is to match these two configurations optimally, with the criterion being least-squares: find a transformation matrix, $\mathbf{T}_{p \times p}$, such that $\parallel \mathbf{X}_{1}\mathbf{T} - \mathbf{X}_{2} \parallel$ is minimized, where $\parallel \cdot \parallel$ denotes the sum-of-squares of the incorporated matrix, i.e., if $\mathbf{A} = \{a_{uv}\}$, then $\parallel \mathbf{A} \parallel$ = $ \mathrm{trace}(\mathbf{A}'\mathbf{A})$ = $\sum_{u,v} a_{uv}^2$. For convenience, assume both $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ have been normalized so $\parallel \mathbf{X}_{1} \parallel$ = $\parallel \mathbf{X}_{2} \parallel$ = 1, and the columns of $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ have sums of zero. \smallskip Two results are central: \smallskip (a) When $\mathbf{T}$ is unrestricted, we have the multivariate multiple regression solution \[ \mathbf{T}^{*} = (\mathbf{X}_{1}'\mathbf{X}_{1})^{-1}\mathbf{X}_{1}'\mathbf{X}_{2} \ ; \] (b) When $\mathbf{T}$ is orthogonal, we have the Sch\"{o}nemann solution done for his thesis in the Quantitative Division at Illinois in 1965 (published in \emph{Psychometrika} in 1966): \smallskip for the SVD of $\mathbf{X}_{2}' \mathbf{X}_{1}$ = $\mathbf{U}\mathbf{S}\mathbf{V}'$, we let $\mathbf{T}^{*} = \mathbf{V} \mathbf{U}'$. \subsection{Matrix Rank Reduction} Lagrange's Theorem (as inappropriately named by C. R. Rao, because it should really be attributed to Guttman) can be stated as follows: \smallskip Let $\mathbf{G}$ be a nonnegative-definite (i.e., a symmetric positive semi-definite) matrix of order $n \times n$ and of rank $r > 0$. Let $\mathbf{B}$ be of order $n \times s$ and such that $\mathbf{B'GB}$ is non-singular. Then the residual matrix \begin{eqnarray} \mathbf{G_{1} = G - GB(B'GB)^{-1}B'G} \end{eqnarray} \noindent is of rank $r-s$ and is nonnegative definite. \smallskip \noindent Intuitively, this theorem allows you to ``take out'' ``factors'' from a covariance (or correlation) matrix. \smallskip There are two somewhat more general results (from Guttman) on matrix rank reduction that prove useful: \smallskip \noindent Let $\mathbf{S}$ be any matrix of order $n \times N$ and of rank $r > 0$. Let $\mathbf{X}$ and $\mathbf{Y}$ be of orders $s \times n$ and $s \times N$, respectively (where $s \le r$), and such that $\mathbf{XSY}'$ is nonsingular. Then the residual matrix \[ \mathbf{S}_{1} = \mathbf{S} - \mathbf{S}\mathbf{Y}'(\mathbf{X}\mathbf{S}\mathbf{Y}')^{-1}\mathbf{X}\mathbf{S} \] is exactly of rank $r-s$. \smallskip \noindent If $\mathbf{S}$ is of order $n \times N$ and of rank $r$, $\mathbf{F}$ of order $n \times r$ (and of rank $r$), and $\mathbf{SS}' = \mathbf{FF}'$, then there is a unique matrix $\mathbf{P}$ of order $r \times N$ such that \[ \mathbf{S} = \mathbf{FP} \ .\] The matrix $\mathbf{P} = (\mathbf{F}'\mathbf{F})^{-1}\mathbf{F}'\mathbf{S}$ satisfies $\mathbf{PP}' = \mathbf{I}$ (i.e., $\mathbf{P}$ has orthonormal rows). \subsection{Torgerson Metric Multidimensional Scaling} Let $\mathbf{A}$ be a symmetric matrix of order $n \times n$. Suppose we want to find a matrix $\mathbf{B}$ of rank 1 (of order $n \times n$) in such a way that the sum of the squared discrepancies between the elements of $\mathbf{A}$ and the corresponding elements of $\mathbf{B}$ (i.e., $\sum_{j=1}^{n} \sum_{i=1}^{n} (a_{ij} - b_{ij})^{2}$) is at a minimum. It can be shown that the solution is $\mathbf{B} = \lambda \mathbf{kk'}$ (so all columns in $\mathbf{B}$ are multiples of $\mathbf{k}$), where $\lambda$ is the largest eigenvalue of $\mathbf{A}$\ and $\mathbf{k}$ is the corresponding normalized eigenvector. This theorem can be generalized. Suppose we take the first $r$ largest eigenvalues and the corresponding normalized eigenvectors. The eigenvectors are collected in an $n \times r$ matrix $\mathbf{K} = \{\mathbf{k}_{1},\ldots,\mathbf{k}_{r}\}$ and the eigenvalues in a diagonal matrix $\mathbf{\Lambda}$. Then $\mathbf{K \Lambda K'}$ is an $n \times n$ matrix of rank $r$ and is a least-squares solution for the approximation of $\mathbf{A}$ by a matrix of rank $r$. It is assumed, here, that the eigenvalues are all positive. If $\mathbf{A}$ is of rank $r$ by itself and we take the $r$ eigenvectors for which the eigenvalues are different from zero collected in a matrix $\mathbf{K}$ of order $n \times r$, then $\mathbf{A}$ = $\mathbf{K \Lambda K'}$. Note that $\mathbf{A}$ could also be represented by $\mathbf{A} = \mathbf{LL'}$, where $\mathbf{L} = \mathbf{K \Lambda^{1/2}}$ (we factor the matrix), or as a sum of $r$ $n \times n$ matrices --- $\mathbf{A} = \lambda_{1} \mathbf{k_{1}k_{1}'} + \cdots + \lambda_{r} \mathbf{k_{r}k_{r}'}$. \smallskip \emph{Metric Multidimensional Scaling} -- Torgerson's Model (Gower's Principal Coordinate Analysis) Suppose I have a set of $n$ points that can be perfectly represented spatially in $r$ dimensional space. The $i^{th}$ point has coordinates $(x_{i1},x_{i2},\ldots,x_{ir})$. If $d_{ij} = \sqrt{\sum_{k=1}^{r} (x_{ik} - x_{jk})^{2}}$ represents the Euclidean distance between points $i$ and $j$, then \[ d_{ij}^{*} = \sum_{k=1}^{r} x_{ik}x_{jk}, \mathrm{where} \] \begin{equation} d_{ij}^{*} = - \frac{1}{2} (d_{ij}^{2} - A_{i} - B_{j} + C) ; \end{equation} \[ A_{i} = (1/n)\sum_{j=1}^{n} d_{ij}^{2} ; \] \[ B_{j} = (1/n)\sum_{i=1}^{n} d_{ij}^{2} ; \] \[ C = (1/n^{2})\sum_{i=1}^{n}\sum_{j=1}^{n} d_{ij}^{2} . \] Note that $\{d_{ij}^{*}\}_{n \times n} = \mathbf{XX'}$, where $\mathbf{X}$ is of order $n \times r$ and the entry in the $i^{th}$ row and $k^{th}$ column is $x_{ik}$. So, the Question: If I give you $\mathbf{D} = \{d_{ij}\}_{n \times n}$, find me \emph{a} set of coordinates to do it. The Solution: Find $\mathbf{D}^{*} = \{d_{ij}^{*}\}$, and take its Spectral Decomposition. This is \emph{exact} here. \smallskip To use this result to obtain a spatial representation for a set of $n$ objects given \emph{any} ``distance-like'' measure, $p_{ij}$, between objects $i$ and $j$, we proceed as follows: (a) Assume (i.e., pretend) the Euclidean model holds for $p_{ij}$. (b) Define $p_{ij}^{*}$ from $p_{ij}$ using (2). (c) Obtain a spatial representation for $p_{ij}^{*}$ using a suitable value for $r$, the number of dimensions (at most, $r$ can be no larger than the number of positive eigenvalues for $\{p_{ij}^{*} \}_{n \times n}$): \[ \{p_{ij}^{*} \} \approx \mathbf{XX'} \] (d) Plot the $n$ points in $r$ dimensional space. \subsection{A Guttman Multidimensional Scaling Result} If $\mathbf{B}$ is a symmetric matrix of order $n$, having all its elements non-negative, the following quadratic form defined by the matrix $\mathbf{A}$ must be positive semi-definite: \[ \sum_{i,j} b_{ij}(x_{i} - x_{j})^{2} = \sum_{i,j} x_{i} a_{ij} x_{j} , \] where \[ a_{ij} = \left\{ \begin{array}{cc} \sum_{k=1; k \ne i}^{n} b_{ik} & (i = j) \\ -b_{ij} & (i \ne j) \\ \end{array} \right. \] \noindent If all elements of $\mathbf{B}$ are positive, then $\mathbf{A}$ is of rank $n - 1$, and has one smallest eigenvalue equal to zero with an associated eigenvector having all constant elements. Because all (other) eigenvectors must be orthogonal to the constant eigenvector, the entries in these other eigenvectors must sum to zero. \smallskip This Guttman result can be used for a method of multidimensional scaling (mds), and is one that seems to get reinvented periodically in the literature. Generally, this method has been used to provide rational starting points in iteratively-defined nonmetric mds. More recently, the Guttman strategy (although not attributed to him as such) has been applied to graphs and the corresponding 0/1 adjacency matrix (treated as a similarity measure). In this case, we have what are called Laplacian eigenmaps, where the graphs are imbedded into a space by using the coordinates from the \emph{smallest} nonzero eigenvectors. \subsection{A Few General MATLAB Routines to Know About} For Eigenvector/Eigenvalue Decompositions: \smallskip $ [\mathbf{V},\mathbf{D}] = \mathrm{eig}(\mathbf{A})$, where $\mathbf{A}$ = $\mathbf{V}\mathbf{D}\mathbf{V}'$, for $\mathbf{A}$ square; $\mathbf{V}$ is orthogonal and contains eigenvectors (as columns); $\mathbf{D}$ is diagonal and contains the eigenvalues. \smallskip For Singular Value Decompositions: \smallskip $ [\mathbf{U},\mathbf{S},\mathbf{V}] = \mathrm{svd}(\mathbf{B})$, where $\mathbf{B}$ = $\mathbf{U}\mathbf{S}\mathbf{V}'$; the columns of $\mathbf{U}$ and the rows of $\mathbf{V}'$ are orthonormal; $\mathbf{S}$ is diagonal and contains the non-negative singular values (ordered from \emph{largest to smallest}). \smallskip The help comments for the Procrustes routine in the Statistics Toolbox are given verbatim below. Note the very general transformation provided in the form of a MATLAB Structure that involves optimal rotation, translation, and scaling. \small \begin{verbatim} help procrustes procrustes Procrustes Analysis D = procrustes(X, Y) determines a linear transformation (translation, reflection, orthogonal rotation, and scaling) of the points in the matrix Y to best conform them to the points in the matrix X. The "goodness-of-fit" criterion is the sum of squared errors. procrustes returns the minimized value of this dissimilarity measure in D. D is standardized by a measure of the scale of X, given by sum(sum((X - repmat(mean(X,1), size(X,1), 1)).^2, 1)) i.e., the sum of squared elements of a centered version of X. However, if X comprises repetitions of the same point, the sum of squared errors is not standardized. X and Y are assumed to have the same number of points (rows), and procrustes matches the i'th point in Y to the i'th point in X. Points in Y can have smaller dimension (number of columns) than those in X. In this case, procrustes adds columns of zeros to Y as necessary. [D, Z] = procrustes(X, Y) also returns the transformed Y values. [D, Z, TRANSFORM] = procrustes(X, Y) also returns the transformation that maps Y to Z. TRANSFORM is a structure with fields: c: the translation component T: the orthogonal rotation and reflection component b: the scale component That is, Z = TRANSFORM.b * Y * TRANSFORM.T + TRANSFORM.c. [...] = procrustes(..., 'Scaling',false) computes a procrustes solution that does not include a scale component, that is, TRANSFORM.b == 1. procrustes(..., 'Scaling',true) computes a procrustes solution that does include a scale component, which is the default. [...] = procrustes(..., 'Reflection',false) computes a procrustes solution that does not include a reflection component, that is, DET(TRANSFORM.T) is 1. procrustes(..., 'Reflection','best') computes the best fit procrustes solution, which may or may not include a reflection component, 'best' is the default. procrustes(..., 'Reflection',true) forces the solution to include a reflection component, that is, DET(TRANSFORM.T) is -1. Examples: % Create some random points in two dimensions n = 10; X = normrnd(0, 1, [n 2]); % Those same points, rotated, scaled, translated, plus some noise S = [0.5 -sqrt(3)/2; sqrt(3)/2 0.5]; % rotate 60 degrees Y = normrnd(0.5*X*S + 2, 0.05, n, 2); % Conform Y to X, plot original X and Y, and transformed Y [d, Z, tr] = procrustes(X,Y); plot(X(:,1),X(:,2),'rx', Y(:,1),Y(:,2),'b.', Z(:,1),Z(:,2),'bx'); \end{verbatim} \end{document}