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):
cA =
4 10 6 12
14 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 10 6 12
14 6 4 2
10 4 0 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
a =
2
5
1
and b =
4
3
5
(11)
the inner product is
a
T
b =
[] 2 5 1
4
3
5
= [2*4 + 5*3 + 1*5] = 28 (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|| =
√
a
T
a (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
a =
2
5
1
and b =
4
3
5
7
9
(14)
The outer product of a and b is then
ab
T
=
2
5
1
[] 4 3 5 7 9
=
2*4 2*3 2*5 2*7 2*9
5*4 5*3 5*5 5*7 5*9
1*4 1*3 1*5 1*7 1*9
=
8 6 10 14 18
20 15 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 6 10 14 18
20 15 25 35 45
4 3 5 7 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 v
i
, with i = 1, 2, n, then
v
i
T
v
j
=
{
0 for i ≠ j
1 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 A
mxn
by B
nxp
and the result will be
C
mxp
. The ij
th
element of C consists of the inner product of the i
th
row of A with the j
th
column of B. As an example, if
A =
[]
2 5 1
4 5 3
and B =
4 3 5 7
9 5 3 4
5 3 6 7
(17)
then
AB =
[]
2*4+5*9+1*5 2*3+5*5+1*3 2*5+5*3+1*6 2*7+5*4+1*7
4*4+5*9+3*5 4*3+5*5+3*3 4*5+5*3+3*6 4*7+5*4+3*7
=
[]
58 34 31 41
76 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 34 31 41
76 46 53 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
(A
T
)
T
= A (22)
(A + B)
T
= A
T
+ B
T
(23)
(AB)
T
= B
T
A
T
(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
D =
4 0 0 0
0 3 0 0
0 0 7 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
I
4x4
=
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 A
mxn
and identity matrix I of
appropriate size
A
mxn
I
nxn
= I
mxm
A
mxn
= A
mxn
(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 = A
T
(28)
Another way to say this is that A
ij
= A
ji
. 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
2b
1
+b
2
+b
3
=1
4b
1
+3b
2
+4b
3
= 2 (29)
-4b
1
+2b
2
+2b
3
=-6
This system can also be written using matrix multiplication as follows:
2 1 1
4 3 4
-4 2 2
b
1
b
2
b
3
=
1
2
-6
(30)
or simply as
Xb = y (31)
where
X =
2 1 1
4 3 4
-4 2 2
, b =
b
1
b
2
b
3
, and y =
1
2
-6
(32)
The problem is to find the values of b
1
, b
2
and b
3
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 b
1
is eliminated from
them. We see here that subtracting 2 times the first equation from the second equation
should eliminate b
1
. Similarly, subtracting -2 times the first equation from the third
equation will again eliminate b
1
. The result after these first two operations is:
2 1 1
0 1 2
0 4 4
b
1
b
2
b
3
=
1
0
-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 b
2
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 1 1
0 1 2
0 0 -4
b
1
b
2
b
3
=
1
0
-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 b
3
= 1. Substituting this into the second equation yields b
2
= -2.
Substituting b
3
= 1 and b
2
= -2 into the first equation we find that b
1
= 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 1
2 4 5
3 6 1
b
1
b
2
b
3
=
1
3
7
(35)
Elementary row operations can be used to produce the system:
1 1 1
0 0 3
0 0 -2
b
1
b
2
b
3
=
1
1
4
(36)
This system has no solution as the second equation requires that b
3
= 1/3 while the third
equation requires that b
3
= -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 1
2 4 5
3 6 1
b
1
b
2
b
3
=
1
-4
7
(37)
This would be reduced by elementary row operations to:
1 2 1
0 0 3
0 0 -2
b
1
b
2
b
3
=
1
-6
4
(38)
Note that the final two equations have the same solution, b
3
= -2. Substitution of this result
into the first equation yields b
1
+ 2b
2
= 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. b
3
= -2. From the
infinite number of possible solutions for b
1
and b
2
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 1
2 4 5
3 6 1
➔
1 2 1
0 0 3
0 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 3 2
0 0 3
0 0 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
-1
A = 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
-1
A
-1
(42)
This can be extended to hold with three or more matrices, for example:
(ABC)
-1
= C
-1
B
-1
A
-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
A =
2 1 1
4 3 4
-4 2 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
4 3 4 | 0 1 0
-4 2 2 | 0 0 1
➔
2 1 1 | 1 0 0
0 1 2 |-2 1 0
0 0 -4 |10 -4 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 1 0 |
7
2
-1
1
4
0 1 0 | 3 -1
1
2
0 0 -4 |10 -4 1
➔
2 0 0 |
1
2
0
-1
4
0 1 0 | 3 -1
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
4
0
-1
8
0 1 0 | 3 -1
1
2
0 0 1 |
-5
2
1
-1
4
(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 1/4 0 -1/8
0 1 0 3 -1 1/2
0 0 1 -5/2 1 -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:
Không có nhận xét nào:
Đăng nhận xét