Linalg

An Introduction to Linear Algebra Barry M. Wise and Neal B. Gallagher Eigenvector Research, Inc. 830 Wapato Lake Road Manson, WA 98831 USA [email protected] Linear algebra is the language of chemometrics. One cannot expect to truly understand most chemometric techniques without a basic understanding of linear algebra. This article reviews the basics of linear algebra and provides the reader with the foundation required for understanding most chemometrics literature. It is presented in a rather dense fashion: no proofs are given and there is little discussion of the theoretical implications of the theorems and results presented. The goal has been to condense into as few pages as possible the aspects of linear algebra used in most chemometric methods. Readers who are somewhat familiar with linear algebra may find this article to be a good quick review. Those totally unfamiliar with linear algebra should consider spending some time with a linear algebra text. In particular, those by Gilbert Strang are particularly easy to read and understand. Several of the numerical examples in this section are adapted from Strang’s Linear Algebra and Its Applications, Second Edition (Academic Press, 1980). MATLAB (The MathWorks, Inc., Natick MA) commands for performing the operations listed are also included; the reader is encouraged to run the examples presented in the text. Those unfamiliar with MATLAB may wish to read the first few sections of the tutorial chapter of the MATLAB User’s Guide. Scalars, Vectors and Matrices A scalar is a mathematical quantity that is completely described by a magnitude, i.e. a single number. Scalar variables are generally denoted by lowercase letters, e.g. a. Examples of scalar variables include temperature, density, pressure and flow. In MATLAB, a value can be assigned to a scalar at the command line, e.g.

»a = 5; Here we have used the semicolon operator to suppress the echo of the result. Without this semicolon MATLAB would display the result of the assignment:

»a = 5 a= 5 A vector is a mathematical quantity that is completely described by its magnitude and direction. An example of a three dimensional column vector might be  4  b= 3   5 

(1)

Vectors are generally denoted by bold lowercase letters. (In MATLAB no distinction is made between the notation for scalars, vectors and matrices: they can be upper or lower case and bold letters are not used.) In MATLAB, this vector could be entered at the command in one of several ways. One way would be to enter the vector with an element on each line, like this:

»b = [4 3 5] b= 4 3 5 Another way to enter the vector would be to use the semicolon to tell MATLAB that each line was completed. For instance

»b = [4; 3; 5]; produces the same result as above. This vector is represented geometrically in Figure 1, where the three components 4, 3 and 5 are the coordinates of a point in three-dimensional space. Any vector b can be represented by a point in space; there is a perfect match between points and vectors. One can choose to think of a vector as the arrow, the point in space, or as the three numbers which describe the point. In a problem with 400 dimensions (such as in spectroscopy), it is probably easiest to consider the 400 numbers. The transpose of a column vector is a row vector and vice-versa. The transpose is generally indicated by a superscript T, i.e. T, though in some instances, including MATLAB, an apostrophe (') will be used. For example bT = [ 4

3

5 ]

(2)

In MATLAB, we could easily assign bT to another variable c, as follows:

»c = b' c= 4

3

5

A matrix is a rectangular array of scalars, or in some instances, algebraic expressions which evaluate to scalars. Matrices are said to be m by n, where m is the number of rows in the matrix and n is the number of columns. A 3 by 4 matrix is shown here  2 A= 7  5

5 3 2

3 2 0

6  1  3 

(3)

This matrix could be entered in MATLAB as follows:

»A = [2 5 3 6; 7 3 2 1; 5 2 0 3];

0 0 5 4 3 5 4 0 0

0 3 0

Figure 1. Geometric Representation of the Vector [4 3 5]T Matrices are usually denoted by bold uppercase letters. The elements of a matrix can be indicated by their row and column indices, for instance, A2,4 = 1. We can index individual matrix elements in MATLAB in a similar way, for instance:

»A(2,4) ans = 1 The transpose operator “flips” a matrix along its diagonal elements, creating a new matrix with the ith row being equal to the jth column of the original matrix, e.g.  T A =  

2 5 3 6

7 3 2 1

5 2 0 3

   

(4)

This would be done in MATLAB:

»A' ans = 2 5 3 6

7 3 2 1

5 2 0 3

Vectors and scalars are special cases of matrices. A vector is a matrix with either one row or column. A scalar is a matrix with a single row and column. Vector and Matrix Addition and Multiplication Matrices can be added together provided that the dimensions are consistent, i.e. both matrices must have the same number of rows and columns. To add two matrices together simply add the corresponding elements. For instance:

[ 15

4 4

3 0

] + [ 22

4 6

1 3

] = [ 73

8 4 10 3

]

(5)

In MATLAB we might do this as follows:

»x = [1 4 3; 5 4 0]; y = [2 4 1; 2 6 3]; »x + y ans = 3 7

8 10

4 3

Note that if the matrices are not of consistent sizes, the addition will not be defined:

[

1 5

4 4

3 0

]

 2 + 1  6

4  2  = ?? 3 

(6)

If you try this in MATLAB, you will get the following:

»x = [1 4 3; 5 4 0]; y = [2 4; 1 2; 6 3]; »x + y ??? Error using ==> + Matrix dimensions must agree. This is one of the most common error messages you will see when using MATLAB. It tells you that the operation you are attempting is not defined.

Vector and matrix addition is commutative and associative, so provided that A, B and C are matrices with the same number of rows and columns, then: A+B=B+A

(7)

(A + B) + C = A + (B + C)

(8)

Vectors and matrices can be multiplied by a scalar. For instance, if c = 2, then (using our previously defined A):  4 10 6 12  cA =  1 4 6 4 2   10 4 0 6  Multiplication of a matrix by a scalar in MATLAB is done using the * operator:

»c = 2; »c*A ans = 4 14 10

10 6 4

6 4 0

12 2 6

Note that scalar multiplication is commutative and associative, so if c and d are scalars: cA = Ac

(9)

(c + d)A = cA + dA

(10)

Vectors can be multiplied together in two different ways. The most common vector multiplication is the inner product. In order for the inner product to be defined, the two vectors must have the same number of elements or dimension. The inner product is the sum over the products of corresponding elements. For instance, for two 3 element column vectors a and b with  4   2  a =  5  and b =  3   5   1 

(11)

the inner product is aTb = [ 2

5

 4  1 ]  3  = [2*4 + 5*3 + 1*5] = 28  5 

(12)

Note that the inner product of two vectors produces a scalar result. The inner product occurs often in the physical sciences and is often referred to as the dot product. We will see it again in regression problems where we model a system property or output (such as a chemical concentration or quality variable) as the weighted sum of a number of different measurements. In MATLAB, this example looks like:

»a = [2; 5; 1]; b = [4; 3; 5]; »a'*b ans = 28 The inner product is also used when calculating the length of a vector. The length of a vector is square root of the sum of squares of the vector coefficients. The length of a vector a, denoted ||a|| (also known as the 2-norm of the vector), is the square root of the inner product of the vector with itself: ||a|| =  aTa √

(13)

We could calculate the norm of the vector a above explicitly in MATLAB with the sqrt (square root) function. Alternately, we could use the norm function. Both methods are shown here:

»sqrt(a'*a) ans = 5.4772 »norm(a) ans = 5.4772 The vector outer product is encountered less often than the inner product, but will be important in many instances in chemometrics. Any two vectors of arbitrary dimension (number of elements) can be multiplied together in an outer product. The result will be an m by n matrix were m is the number of elements in the first vector and n is the number of elements in the second vector. As an example, take 4  2 3   a =  5  and b =  5  1   7 9

 

(14)

The outer product of a and b is then  2  abT =  5  [ 4  1 

3

5

7

 2*4 2*3 2*5 2*7 2*9  9 ] =  5*4 5*3 5*5 5*7 5*9   1*4 1*3 1*5 1*7 1*9 

6 10 14 18   8 2 0 1 5 25 35 45  =  4 3 5 7 9 

(15)

In MATLAB this would be:

»a = [2 5 1]'; b = [4 3 5 7 9]'; »a*b' ans = 8 20 4

6 15 3

10 25 5

14 35 7

18 45 9

Here we used a slightly different method of entering a column vector: it was entered as the transpose of a row vector. Orthogonal and Orthonormal Vectors Vectors are said to be orthogonal if their inner product is zero. The geometric interpretation of this is that they are at right angles to each other or perpendicular. Vectors are orthonormal if they are orthogonal and of unit length, i.e. if their inner product with themselves is unity. For an orthonormal set of column vectors vi, with i = 1, 2, ... n, then v iT v j =

for i ≠ j { 10 for i = j

(16)

Note that it is impossible to have more than n orthonormal vectors in an n-dimensional space. For instance, in three dimensions, one can only have 3 vectors which are orthogonal, and thus, only three vectors can be orthonormal (although there are an infinite number of sets of 3 orthogonal vectors). Sets of orthonormal vectors can be thought of as new basis vectors for the space in which they are contained. Our conventional coordinate system in 3 dimensions consisting of the vectors [1 0 0]T, [0 1 0]T and [0 0 1]T is the most common orthonormal basis set, however, as we shall see, it is not the only basis set, nor is it always the most convenient for describing real data. Two matrices A and B can be multiplied together provided that the number of columns in A is equal to the number of rows in B. The result will be a matrix C which has as many rows as A and columns as B. Thus, one can multiply Amxn by Bnxp and the result will be Cmxp. The ijth element of C consists of the inner product of the ith row of A with the jth column of B. As an example, if A=

[

2 4

5 5

1 3

]

 4 and B =  9  5

3 5 3

5 3 6

7  4  7 

(17)

then AB =

[ 2*4+5*9+1*5 4*4+5*9+3*5

2*3+5*5+1*3 2*5+5*3+1*6 2*7+5*4+1*7 4*3+5*5+3*3 4*5+5*3+3*6 4*7+5*4+3*7

]

[ 57 86

=

34 31 41 46 53 69

]

(18)

In MATLAB:

»A = [2 5 1; 4 5 3]; B = [4 3 5 7; 9 5 3 4; 5 3 6 7]; »A*B ans = 58 76

34 46

31 53

41 69

Matrix multiplication is distributive and associative A(B + C) = AB + AC

(19)

(AB)C = A(BC)

(20)

but is not, in general, commutative AB

≠ BA

(21)

Other useful matrix identities include the following (AT)T = A

(22)

(A + B)T = AT + BT

(23)

(AB)T = BTAT

(24)

Special Matrices There are several special matrices which deserve some attention. The first is the diagonal matrix, which is a matrix whose only non-zero elements lie along the matrix diagonal, i.e. the elements for which the row and column indices are equal. An example diagonal matrix D is shown here  4 D= 0  0

0 3 0

0 0 7

0  0  0 

(25)

Diagonal matrices need not have the same number of rows as columns, the only requirement is that only the diagonal elements be non-zero. We shall see that diagonal matrices have additional properties which we can exploit. Another special matrix of interest is the identity matrix, typically denoted as I. The identity matrix is always square (number of rows equals number of columns) and contains ones on the diagonal and zeros everywhere else. A 4 by 4 identity matrix is shown here

 I4x4 =   

1 0 0 0

0 1 0 0

0 0 1 0

0 0 0 1

   

(26)

Identity matrices can be created easily in MATLAB using the eye command. Here we create the 3 by 3 identity matrix:

»id = eye(3) id = 1 0 0

0 1 0

0 0 1

Any matrix multiplied by the identify matrix returns the original matrix (provided, of course, that the multiplication is defined). Thus, for a matrix Amxn and identity matrix I of appropriate size AmxnInxn = ImxmAmxn = Amxn

(27)

The reader is encouraged to try this with some matrices in MATLAB. A third type of special matrix is the symmetric matrix. Symmetric matrices must be square and have the property that they are equal to their transpose: A = AT

(28)

Another way to say this is that Aij = Aji. We will see symmetric matrices turn up in many places in chemometrics. Gaussian Elimination: Solving Systems of Equations Systems of equations arise in many scientific and engineering applications. Linear algebra was largely developed to solve these important problems. Gaussian elimination is the method which is consistently used to solve systems of equations. Consider the following system of 3 equations in 3 unknowns 2b1 4b1 -4b1

+ + +

b2 3b2 2b2

+ + +

b3 4b3 2b3

= = =

1 2 -6

(29)

This system can also be written using matrix multiplication as follows:  2 1  4 3 -4 2

1   b1   1  4   b2  =  2  2   b3  -6 

(30)

or simply as Xb = y

(31)

where  2 1 X= 4 3 -4 2

1 4 2

 b1    1  b , b = , and y =  2    2   b3   -6 

(32)

The problem is to find the values of b1, b2 and b3 that make the system of equations hold. To do this, we can apply Gaussian elimination. This process starts by subtracting a multiple of the first equation from the remaining equations so that b1 is eliminated from them. We see here that subtracting 2 times the first equation from the second equation should eliminate b1 . Similarly, subtracting -2 times the first equation from the third equation will again eliminate b1. The result after these first two operations is:  2  0  0

1 1 4

1   b1   1  2   b2  =  0  4   b3  -4 

(33)

The first coefficient in the first equation, the 2, is known as the pivot in the first elimination step. We continue the process by eliminating b2 in all equations after the second, which in this case is just the third equation. This is done by subtracting 4 times the second equation from the third equation to produce:  2  0  0

1 1   b1   1  1 2   b2  =  0  0 -4   b3  -4 

(34)

In this step the 1, the second coefficient in the second equation, was the pivot. We have now created an upper triangular matrix, one where all elements below the main diagonal are zero. This system can now be easily solved by back-substitution. It is apparent from the last equation that b3 = 1. Substituting this into the second equation yields b2 = -2. Substituting b3 = 1 and b2 = -2 into the first equation we find that b1 = 1. Gaussian elimination is done in MATLAB with the left division operator \ (for more information about \ type help mldivide at the command line). Our example problem can be easily solved as follows:

»X = [2 1 1; 4 3 4; -4 2 2]; y = [1; 2; -6]; »b = X\y b= 1 -2 1 It is easy to see how this algorithm can be extended to a system of n equation in n unknowns. Multiples of the first equation are subtracted from the remaining equations to produce zero coefficients below the first pivot. Multiples of the second equation are then used to clear the coefficients below the second pivot, and so on. This is repeated until the last equation contains only the last unknown and can be easily solved. This result is substituted into the second to last equation to get the second to last unknown. This is repeated until the entire system is solved.

Singular Matrices and the Concept of Rank Gaussian elimination appears straightforward enough, however, are there some conditions under which the algorithm may break down? It turns out that, yes, there are. Some of these conditions are easily fixed while others are more serious. The most common problem occurs when one of the pivots is found to be zero. If this is the case, the coefficients below the pivot cannot be eliminated by subtracting a multiple of the pivot equation. Sometimes this is easily fixed: an equation from below the pivot equation with a non-zero coefficient in the pivot column can be swapped with the pivot equation and elimination can proceed. However, if all of the coefficients below the zero pivot are also zero, elimination cannot proceed. In this case, the matrix is called singular and there is no hope for a unique solution to the problem. As an example, consider the following system of equations:  1  2  3

2 4 6

1   b1   1  5   b2  =  3  1   b3   7 

(35)

Elementary row operations can be used to produce the system:  1  0  0

1 1   b1   1  0 3   b2  =  1  0 -2   b3   4 

(36)

This system has no solution as the second equation requires that b3 = 1/3 while the third equation requires that b3 = -2. If we tried this problem in MATLAB we would see the following:

»X = [1 2 1; 2 4 5; 3 6 1]; y = [1; 1; 4]; »b = X\y Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.745057e-18. ans = 1.0e+15 * 2.0016 -1.0008 -0.0000 We will discuss the meaning of the warning message shortly. On the other hand, if the right hand side of the system were changed so that we started with:  1  2  3

2 4 6

1   b1   1  5   b2  = -4  1   b3   7 

(37)

This would be reduced by elementary row operations to:  1  0  0

2 1   b1   1  0 3   b2  = -6  0 -2   b3   4 

(38)

Note that the final two equations have the same solution, b3 = -2. Substitution of this result into the first equation yields b1 + 2b2 = 3, which has infinitely many solutions. We can see that singular matrices can cause us to have no solution or infinitely many solutions. If we do this problem in MATLAB we obtain:

»X = [1 2 1; 2 4 5; 3 6 1]; y = [1; -4; 7]; »b = X\y Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.745057e-18. ans = -14.3333 8.6667 -2.0000 The solution obtained agrees with the one we calculated above, i.e. b3 = -2. From the infinite number of possible solutions for b1 and b2 MATLAB chose the values -14.3333 and 8.6667, respectively. Consider once again the matrix from our previous example and its reduced form:  1  2  3

2 4 6

1   1 5 ➔ 0 1   0

2 1  0 3  0 -2 

(39)

We could have taken the elimination one step further by subtracting -2/3 of the second equation from the third to produce:  1  0  0

3 0 0

2  3  0 

(40)

This is known as the echelon form of a matrix. It is necessarily upper triangular, that is, all non-zero entries are confined to be on, or above, the main diagonal (the pivots are the 1 in the first row and the 3 in the second row). The non-zero rows come first, and below each pivot is a column of zeros. Also, each pivot lies to the right of the pivot in the row above. The number of non-zero rows obtained by reducing a matrix to echelon form is the rank of the matrix. Note that this procedure can be done on any matrix; it need not be square. It can be shown that the rank of a m by n matrix has to be less than, or equal to, the lessor of m or n, i.e. rank(X) ≤ min(m,n)

(41)

A matrix whose rank is equal to the lessor of m or n is said to be of full rank. If the rank is less than min(m,n), the matrix is rank deficient or collinear. MATLAB can be used to determine the rank of a matrix. In this case we get:

»rank(X) ans = 2 which is just what we expect. Matrix Inverses The inverse of matrix A, A-1, is the matrix such that AA-1 = I and A-1A = I. It is assumed here that A, A-1 and I are all square and of the same dimension. Note, however, that A-1 might not exist. If it does, A is said to be invertible, and its inverse A-1 will be unique. It can also be shown that the product of invertible matrices has an inverse, which is the product of the inverses in reverse order: (AB)-1 = B-1A-1

(42)

This can be extended to hold with three or more matrices, for example: (ABC)-1 = C-1B-1A-1

(43)

How do we find matrix inverses? We have already seen that we can perform elementary row operations to reduce a matrix A to an upper triangular echelon form. However, provided that the matrix A is nonsingular, this process can be continued to produce a diagonal matrix, and finally, the identity matrix. This same set of row operations, if performed on the identity matrix produces the inverse of A, A-1. Said another way, the same set of row operations that transform A to I transform I to A -1. This process is known as the Gauss-Jordan method. Let us consider an example. Suppose that we want to find the inverse of  2 1 A= 4 3 -4 2

1  4  2 

(44)

We will start by augmenting A with the identity matrix, then we will perform row operations to reduce A to echelon form:  2 1 1 | 1 0 0   2 1 1 | 1 0 0   4 3 4 | 0 1 0  ➔  0 1 2 |- 2 1 0  ➔  0 0 - 4 |10 - 4 1  -4 2 2 | 0 0 1  We now continue the operation by eliminating the coefficients above each of the pivots. This is done by subtracting the appropriate multiple of the third equation from the first two, then subtracting the second equation from the first:

 2  0

0

7 1 1 0 | 2 -1 4 1 1 0 | 3 -1 2 0 -4 |10 -4 1

  2 ➔ 0  

1 -1 0 0 | 2 0 4 1 1 0 | 3 -1 2

0

0 - 4 |10 - 4

1

 ➔ 

Finally, each equation is multiplied by the inverse of the diagonal to get to the identity matrix:

1 0  0

1 0 0 | 4

   -1 4 

-1 8 1 1 0 | 3 -1 2 0

-5 0 1 | 2 1

(45)

The matrix on the right is A-1, the inverse of A. Note that our algorithm would have failed if we had encountered an incurable zero pivot, i.e. if A was singular. In fact, a matrix is invertible if, and only if, it is nonsingular. This process can be repeated in MATLAB using the rref (reduced row echelon form) command. We start by creating the A matrix. We then augment it with the identify matrix using the eye command and use the rref function to create the reduced row echelon form. In the last step we verify that the 3 by 3 matrix on the right is the inverse of A. (The result is presented here as fractions rather than decimals by setting the MATLAB command window display preference to “rational”).

»A = [2 1 1; 4 3 4; -4 2 2]; »B = rref([A eye(3)]) B= 1 0 0

0 1 0

0 0 1

1/4 3 -5/2

0 -1 1

-1/8 1/2 -1/4

»A*B(:,4:6) ans = 1 0 0

0 1 0

0 0 1

A simpler way of obtaining an inverse in MATLAB is to use the inv command:

»Ainv = inv(A) Ainv = 0.2500 3.0000 -2.5000

0.0000 -1.0000 1.0000

-0.1250 0.5000 -0.2500

Here the results are shown in short, rather than rational format. (Note that, in rational format, a * may appear in place of the 0 in the first row due to roundoff error.) Another useful property of inverses is: (AT)-1 = (A-1)T

(46)

We could also verify this for our example in MATLAB:

»inv(A') - inv(A)' ans = 1.0e-15 * -0.3331 0.1110 -0.0555

0.4441 0 0.1110

0 0 0

which is the zero matrix to within machine precision. Vector Spaces and Subspaces Before continuing on it is useful to introduce the concept of a vector space. The most important vector spaces are denoted R1, R2, R3, ....; there is one for every integer. The space Rn consists of all column vectors with n (real) components. The exponent n is said to be the dimension of the space. For instance, the R3 space can be thought of as our “usual” three-dimensional space: each component of the vector is one of the x, y or z axis coordinates. Likewise, R2 defines a planar space, such as the x-y plane. A subspace is a vector space which is contained within another. For instance, a plane within R3 defines a two dimensional subspace of R 3 , and is a vector space in its own right. In general, a subspace of a vector space is a subset of the space where the sum of any two vectors in the subspace is also in the subspace and any scalar multiple of a vector in the subspace is also in the subspace. Linear Independence and Basis Sets In a previous section we defined the rank of a matrix in a purely computational way: the number of non-zero pivots after a matrix had been reduced to echelon form. We are beginning to get the idea that rank and singularity are important concepts when working with systems of equations. To more fully explore these relationships the concept of linear

independence must be introduced. Given a set of vectors v1, v2, ... , vk, if all non-trivial combinations of the vectors are nonzero c1v1 + c2v2 + ... + ckvk ≠ 0 unless c1 = c2 = ... = ck = 0

(47)

then the vectors are linearly independent. Otherwise, at least one of the vectors is a linear combination of the other vectors and they are linearly dependent. It is easy to visualize linear independence. For instance, any two vectors that are not scalar multiples of each other are linearly independent and define a plane. A third vector in that plane must be linearly dependent as it can be expressed as a weighted sum of the other two vectors. Likewise, any 4 vectors in R3 must be linearly dependent. Stated more formally, a set of k vectors in Rm must be linearly dependent if k>m. It is also true that the r nonzero rows of an echelon matrix are linearly independent, and so are the r columns that contain the pivots. A set of vectors w 1, w 2, ... , w k, in R n is said to span the space if every vector v in R n can be expressed as a linear combination of w’s, i.e. v = c1w 1 + c2w 2 + ... + ckw k for some ci.

(48)

Note that for the set of w’s to span R n then lk≥ n. A basis for a vector space is a set of vectors which are linearly independent and span the space. Note that the number of vectors in the basis must be equal to the dimension of the space. Any vector in the space can be specified as one and only one combination of the basis vectors. Any linearly independent set of vectors can be extended to a basis by adding (linearly independent) vectors so that the set spans the space. Likewise, any spanning set of vectors can be reduced to a basis by eliminating linearly dependent vectors. Row Spaces, Column Spaces and Nullspaces Suppose that U is the echelon matrix obtained by performing elimination on the m by n matrix A. Recall that the rank r of A is equal to the number of nonzero rows of U. Likewise, the dimension of the row space of A (the space spanned by the rows of A), denoted R(AT), is equal to r, and the rows of U form a basis for the row space of A, that is, they span the same space. This is true because elementary row operations leave the row space unchanged. The column space of A (the space spanned by the columns of A, also referred to as the range of A), denoted R (A), also has dimension r. This implies that the number of independent columns in A equals the number of independent rows, r. A basis can be constructed for the column space of A by selecting the r columns of U which contain nonzero pivots. The fact that the dimension of the row space is equal to the dimension of the column space is one of the most important theorems in linear algebra. This fact is often expressed simply as “row rank = column rank.” This fact is generally not obvious for a random non-square matrix. For square matrices it implies that, if the rows are linearly independent, so are the columns. Note also that the column space of A is equal to the row space of AT.

The nullspace of A, N(A), is of dimension n - r. N(A) is the space of Rn not spanned by the rows of A. Likewise, the nullspace of AT, N(AT), (also known as the left nullspace of A) has dimension m - r, and is the space of Rm not spanned by the columns of A. Orthogonality of Subspaces Earlier we defined orthogonal vectors. Recall that any two vectors v and w are orthogonal if their inner product is zero. Two subspaces V and W can also be orthogonal, provided that every vector v in V is orthogonal to every vector w in W , i.e. vTw = 0 for all v and w. Given the definition of othogonality of subspaces, we can now state more clearly the properties of the row space, column space and nullspaces defined in the previous section. For any m by n matrix A, the nullspace N (A) and the row space R (A T) are orthogonal subspaces of R n. Likewise, the left nullspace N (A T ) and the column space R (A) are orthogonal subspaces of Rm. The orthogonal complement of a subspace V of Rn is the space of all vectors orthogonal to V and is denoted V⊥ (pronounced V perp). Projections onto Lines The projection of points onto lines is very important in many chemometric methods. In general, the projection problem is, given a vector x and a point defined by the vector y, find the point p along the direction defined by x which is closest to y. This is illustrated graphically in Figure 2. Finding this point p is relatively straightforward once we realize that p must be a scalar multiple of x, i.e. p = b^x, (here we use the the hat symbol ^ to denote estimated variables such as b^), and that the line connecting y to p must be perpendicular to x, i.e. xTy x T (y - b^x) = 0 ➔ xTy = b^xTx ➔ b^ = T x x

(49)

Therefore, the projection p of the point y onto the line spanned by the vector x is given by xTy p = b^x = T x x x

(50)

The (squared) distance from y to p can also be calculated. The result is ||y - p||2 =

(yTy)(xTx)-(xTy)2 xTx

(51)

y x p

Figure 2. The Projection of the vector y onto the vector x. We may also be interested in the projection of a point y onto a subspace. For instance, we may want to project the point onto the plane defined by two vectors, or an n dimensional subspace defined by a collection of vectors, i.e. a matrix. In this case, the vector x in the two equations above will become the matrix X, and the inner product xTx will become XTX. If X is m by n and rank r, the matrix XTX is n by n, symmetric (XTX = (XTX)T) and also of rank r. Thus we see that if X has linearly independent columns, i.e. rank r = n, then X TX is nonsingular and therefore invertible. This is a very important result as we shall see shortly. Projections onto Subspaces and Least Squares Previously we have considered solving Xb = y only for systems where there are the same number of equations as variables (m = n). This system either has a solution or it doesn’t. We now consider the case of m > n, where we have more samples (rows) than variables (columns), with the possibility that the system of equations can be inconsistent, i.e. there is no solution which makes the equality hold. In these cases, one can only solve the system such that the average error in the equations is minimized. Consider for a moment the single variable case xb = y where there is more than 1 equation. The error is nothing more than the length of the vector xb - y, E = ||xb - y||. It is easier, however, to minimize the squared length, E 2 = (xb - y)T(xb - y) = xTxb2 - 2xTyb + yTy

(52)

The minimum of this can be found by taking the derivative with respect to b and setting it to zero: dE2 T T db = 2x xb - 2x y = 0 The least squares solution to the problem is then

(53)

xTy b^ = T x x

(54)

This is the same solution we obtained to the projection problem p = b^x in equation (49). Thus, requiring that the error vector be perpendicular to x gave the same result as the least squares approach. This suggests that we could solve systems in more than one variable by either geometry or calculus. Consider now the system Xb = y, where X is m by n, m > n. Let us solve this system by requiring that Xb^ - y be perpendicular to the column space of X. For this to be true, each vector in the column space of X must be perpendicular to Xb^ - y. Each vector in the column space of X, however, is expressible as a linear combination, say c, of the columns of X, Xc. Thus, for all choices of c, (Xc)T (Xb^ - y) = 0, or cT[X TXb^ - X Ty] = 0

(55)

There is only one way that this can be true for all c: the vector in the brackets must be zero. Therefore, it must be true that X T X b^ = X T y

(56)

These are known as the normal equations for least squares problems. We can now solve for b^ by multiplying through by (XTX)-1 to obtain: b^ = (X T X)-1X T y

(57)

The projection of y onto the column space of X is therefore p = Xb^ = X(X T X)-1X T y

(58)

In chemometrics we commonly refer to b^ as the regression vector. Often X consists of the matrix of m samples with n measured variables. The vector y consists of some property (such as a concentration or quality parameter) for the m samples. We then want to determine a b^ that we can use with new samples, Xnew, to estimate their properties y^new: y^new = Xnewb^

(59)

As an example, suppose we measured two variables, x1 and x2 on four samples, and also had a quality variable y. Our data would then be:  X=  

1 1 2 2

1 2 1 2

 6    and y =  6   7    11  

(60)

We can determine a regression vector b^ that relates the measured variables x1 and x2 to the quality variable y in MATLAB as follows:

»X = [1 1; 1 2; 2 1; 2 2]; y = [6 6 7 11]'; »b = inv(X'*X)*X'*y b= 3.0000 2.0000 Thus, our model for estimating our quality variable is y = 3x1 + 2x2 . In MATLAB, however, we can simplify this calculation by using the “left division” symbol \ instead of the inv operator and get the same result:

»b = X\y b= 3.0000 2.0000 The b^ can now be used to get p, the projection of y into X, and the difference between y and its projection into X, d:

»p = X*b p= 5 7 8 10 »d = y-p d= 1 -1 -1 1 We can also check to see if d is orthogonal to the columns of X. The simplest way to do this is to check the value of XTd:

»X'*d ans = 1.0e-14 * -0.9770 -0.9770 This gives us the inner product of each of the columns of X with d. We see that the values are very close to zero, the only difference from zero being due to machine round off error. Thus, d is orthogonal to the columns of X, as it should be. Ill-conditioned Problems We have now seen how to solve least squares problems in several variables. Once again, we might ask: “when might I expect this procedure to fail?” Based on our previous experience with the solution of m equations in m unknowns and the formation of inverses, we expect that the least squares procedure will fail when X is rank deficient, that is, when rank(X T X) < n. However, what about when X is nearly rank deficient, i.e. when one column of X is nearly a linear combination of the other columns? Let us consider such a system:  X=  

1 2 2 4 3 6 4 8.0001

   and y =     

2 4 6 8

   

(61)

Here the second column of X is almost, but not quite, twice the first column. In other words, X is close to being singular. We can now use MATLAB to estimate a regression vector b^ for this problem:

»X = [1 2; 2 4; 3 6; 4 8.0001]; y = [2 4 6 8]'; »b = X\y b= 2 0 Here we see that the system can be solved exactly by y = 2x1. However, what would happen if the y data changed a very small amount? Consider changing the third element of y from 6 to 5.9999 and 6.0001. The following results:

»X = [1 2; 2 4; 3 6; 4 8.0001]; y = [2 4 6.0001 8]'; b = X\y b= 3.7143 -0.8571 »X = [1 2; 2 4; 3 6; 4 8.0001]; y = [2 4 5.9999 8]'; b = X\y b= 0.2857 0.8571 We see that the regression coefficients changed from b^ = [3.7143 -0.8571] to b^ = [0.2857 0.8571] given this very small change in the data. This is an example of what happens in regression when the X matrix is very close to singular. Methods for dealing with these illconditioned problems will be discussed in the next section. Projection Matrices We saw previously that given the problem Xb = y, we could obtain the projection of y onto the columns of X, p, from: p = X(X T X)-1X T y

(62)

p = Py

(63)

The quantity X(X TX)-1X T, which we shall denote P, is a projection matrix. Projection matrices have two important properties: they are idempotent and symmetric. Thus PP = P 2 = P

(64)

PT = P

(65)

It is easy to see that P must be idempotent: the first application of P projects y into X. A second application won’t change the result since Py will already be in the column space of X. Orthogonal and Orthonormal Bases When we think of a basis for a space, such as the x-y plane, we normally think of an orthogonal basis, i.e. we define our distances in the plane based on two directions that are perpendicular. Orthogonal basis sets, it turns out, are convenient from both a conceptual and mathematical perspective. Basis sets can be made even simpler by normalizing the basis vectors to unit length. The result is an orthonormal basis. An orthonormal basis, v1, v2, ..., vk has the property that:

v iT v j =

for i ≠ j { 10 for i = j

(66)

Of course, we are most familiar with the standard basis, v 1 = [1 0 ... 0]T , v 2 = [0 1 ... 0]T, etc., however, it is possible to construct an infinite number of orthonormal bases for any space of dimension 2 or greater by simply rotating the basis vectors while keeping them at right angles to each other. Suppose for a moment that we are interested in projecting a vector y onto an X that consists of orthonormal columns. In this case, XTX = I, thus P = X(X T X)-1X T = XX T

(67)

A square matrix with orthonormal columns is known as an orthogonal matrix (though orthonormal matrix would have been a better name). An orthogonal matrix Q has the following properties: QTQ = I, QQT = I, and QT = Q-1

(68)

Furthermore, Q will also have orthonormal rows, i.e. if Q is an orthogonal matrix, so is Q T. We will see this again shortly when we consider Principal Components Analysis (PCA) in the next section. Pseudoinverses We saw previously that were not able to solve Xb = y using the “normal equations” in the case where the columns of X were not independent or collinear. But is that the best we can do? Are we simply forced to conclude that there is no solution to the problem? In fact, it is still possible to obtain a least squares solution to the problem. To do so, we must introduce the concept of the pseudoinverse, X +, for a matrix X that is not invertible. The general situation when X is collinear is that there are many solutions to Xb = y. Which one should we choose?: The one with minimum length of b^, ||b^||. To get this we require that the component of b^ in the nullspace of X be zero, which is the same as saying that b^ must lie in the row space of X. Thus, the optimal least squares solution to any Xb = y is b^ such that: X b^ equals the projection of y into the column space of X, and b^ lies in the row space of X. We call the matrix that solves this problem the pseudoinverse X+, where b^ = X + y. Note that when the columns of X are independent, X + = (X T X)-1X T

(69)

This satisfies our requirements stated above. In the case where X is collinear, we must find another way to calculate X+.

The Singular Value Decomposition We now introduce the singular value decomposition, or SVD, which will allow us to obtain a general solution to the problem of finding pseudoinverses. Any m by n matrix X can be factored into X = USV T

(70)

where U is orthogonal and m by m, V is orthogonal and n by n, and S is m by n and diagonal. The non-zero entries of S are referred to as the singular values and decrease monotonically from the upper left to the lower right of S. As an example of the SVD, let us consider the matrix  X=  

1 2 3 2 3 5 3 5 8 4 8 12

   

(71)

Here we have chosen X to be collinear: the third column is the sum of the first two columns. We can now use MATLAB to compute the SVD of X, and verify that the product of the U, S and VT reproduces X.

»X = [1 2 3; 2 3 5; 3 5 8; 4 8 12]; »[U,S,V] = svd(X) U= 0.1935 0.3184 0.5119 0.7740

0.1403 -0.6426 -0.5022 0.5614

-0.9670 0.0341 -0.0341 0.2503

19.3318 0 0 0

0 0.5301 0 0

0 0 0.0000 0

-0.7661 0.6277 -0.1383

0.5774 0.5774 -0.5774

S=

V= 0.2825 0.5221 0.8047

0.0885 0.6961 -0.6961 0.1519

»U*S*V' ans = 1.0000 2.0000 3.0000 4.0000

2.0000 3.0000 5.0000 8.0000

3.0000 5.0000 8.0000 12.0000

Note that S is diagonal, and that the diagonal elements decrease from upper left to lower right. In fact, the last singular value appears to be zero. Recall now that the inverse of a product is the product of the inverses in reverse order. This leads to an explicit formula for the pseudoinverse of X, X+: X + = VS + U T

(72)

Here we must be very careful in formation of S+. In practice, all singular values that are close to zero are set to zero in the inverse (the challenge comes with determining what constitutes “close to zero”). Alternately, one can think of this as truncation of the matrices U, S and V where only the first r columns of U and V are retained, along with the first r rows and columns of S, where r = rank(X). In MATLAB, we see for our example:

»U(:,1:2)*S(1:2,1:2)*V(:,1:2)' ans = 1.0000 2.0000 3.0000 4.0000

2.0000 3.0000 5.0000 8.0000

3.0000 5.0000 8.0000 12.0000

It is apparent that we only need the first two “factors” in the SVD to reproduce the original data. The inverse of X, X+ can then be calculated from:

»Xinv = V(:,1:2)*inv(S(1:2,1:2))*U(:,1:2)' Xinv = -0.2000 0.1714 -0.0286

0.9333 -0.7524 0.1810

0.7333 -0.5810 0.1524

-0.8000 0.6857 -0.1143

Note that this is same result as would be obtained from the MATLAB command pinv, which uses the SVD and sets all singular values that are within machine precision to zero:

»pinv(X) ans = -0.2000 0.1714 -0.0286

0.9333 -0.7524 0.1810

0.7333 -0.5810 0.1524

-0.8000 0.6857 -0.1143

Let us now return to our ill-conditioned regression problem where  X=  

1 2 2 4 3 6 4 8.0001

   and y =     

2 4 6 8

   

(73)

Recall that a small change in the third element of y produced wildly different regression vectors. We can now repeat this problem assuming that the true rank of X is 1, i.e. to within measurement accuracy the second column of X is a multiple of the first column. In MATLAB we would obtain:

»X = [1 2; 2 4; 3 6; 4 8.0001]; y = [2 4 6 8]'; »[U,S,V] = svd(X); »Xinv = V(:,1)*inv(S(1,1))*U(:,1)' Xinv = 0.0067 0.0133

0.0133 0.0267

0.0200 0.0400

0.0267 0.0533

»b = Xinv*y b= 0.4000 0.8000 With this inverse, changing the third element of y would no longer change the result significantly:

»y = [2 4 5.9999 8]'; »b = Xinv*y b= 0.4000 0.8000 »y = [2 4 6.0001 8]'; »b = Xinv*y b= 0.4000 0.8000 Thus we see that the use of a pseudoinverse can enhance the stability of our calculations. Pseudoinverses plan an important role in many chemometric regression methods such as Principal Components Regression (PCR) and Partial Least Squares (PLS), and many different approaches to the defining the pseudoinverse are used. Conclusion This introduction to linear algebra, while brief, should include most of what the reader needs to know in order to tackle the chemometrics literature. Interested readers are encouraged to consult texts on the subject for more detail on the methods, including proofs and additional examples. Acknowledgement This article was adapted from a section of the PLS_Toolbox 2.0 for use with MATLAB (Wise and Gallagher, Eigenvector Research, Inc, 1998) software package.