Linear Algebra:
Matrices
B-l
Appendix B: LINEAR ALGEBRA: MATRICES
TABLE OF CONTENTS
Page
§B.l      Matrices....................... B-3
§B.1.1      Square Matrices................ B-3
§B.1.2      Symmetry and Antisymmetry............ B-4
§B.1.3      Special Matrices................ B-4
§B. 1.4      *Are Vectors a Special Case of Matrices?........ B-5
§B.2      Elementary Matrix Operations.............. B-6
§B.2.1      Equality................... B-6
§B.2.2      Transposition................. B-6
§B.2.3      Addition, Subtraction and Scaling .......... B-7
§B.3      Matrix Products ................... B-7
§B.3.1      Matrix by Vector Product ............. B-7
§B.3.2      Matrix by Matrix Product and Matrix Powers ...... B-8
§B.3.3      Matrix Product Properties ............. B-9
§B.4      Bilinear and Quadratic Forms .............. B-10
§B.5      Determinants..................... B-l 1
§B.5.1      Determinant Properties.............. B-ll
§B.5.2      Cramer's Rule and Homogeneous Systems........ B-13
§B.6      Singular Matrices, Rank ................ B-14
§B.6.1      Rank Deficiency ................ B-l5
§B.6.2      Rank of Matrix Sums and Products.......... B-15
§B.6.3      Singular Systems: Particular and Homogeneous Solutions      .   . B-16
§B.6.4      *Rank of Rectangular Matrices........... B-16
§B.7      Matrix Inversion.................... B-17
§B.7.1      Explicit Computation of Inverses .......... B-17
§B.7.2      Some Properties of the Inverse............ B-l 8
§B.8      Characteristic Polynomial and Eigenvalues.......... B-19
§B.8.1      *Similarity Transformations and Invariants ....... B-20
§B.8.2      *Congruential Transformations........... B-20
§B.8.3      *The Cayley-Hamilton Theorem........... B-21
§B.9      *The Inverse of a Sum of Matrices ............ B-21
§B.9.1      *The Sherman-Morrison and Woodbury Formulas ..... B-21
§B.9.2      *Formulas for Modified Determinants......... B-22
§B.9.3      *Inversion of Low Order Symmetric Matrix Sums ..... B-22
§B.10    *Matrix Orthogonality................. B-22
§B.10.1     *Matrix Orthogonalization Via Projectors........ B-23
§B.10.2    *Orthogonal Projector Properties .......... B-23
§B.       Notes and Bibliography ................. B-24
§B.       Exercises...................... B-25
B-2
§B.l MATRICES
§B.l. Matrices
This Appendix introduces the concept of a matrix. Consider a set of scalar quantities arranged in a rectangular array containing m rows and n columns:
	an	■ alj ■	■     a\n ~
a2l		■   <hj ■	■ a2n
	ai2 .	■     "ij ■	■ ain
-aml	am2 ■	■     amj ■	amn -
This will be called a rectangular matrix of order mbyn, or, briefly, an m x n matrix. Not every rectangular array is a matrix; to qualify as such it must obey the operational rules discussed below.
The quantities <2,7 are called the entries or components of the matrix.1 As in the case of vectors, the term "matrix element" will be avoided to lessen the chance of confusion with finite elements. The two subscripts identify the row and column, respectively.
Matrices are conventionally identified by bold uppercase letters such as A, B, etc. The entries of matrix A may be denoted as A,7 or <2,7, according to the intended use. Occassionally we shall use the short-hand component notation
A = [aul (B.2)
Example B.l. The following is a 2 x 3 numerical matrix:
(B.3)
2 6 3 4   9 1
This matrix has 2 rows and 3 columns. The first row is (2, 6, 3), the second row is (4, 9, 1), the first column is (2, 4), and so on.
In some contexts it is convenient or useful to display the number of rows and columns. If this is so we will write them underneath the matrix symbol.2 For the example matrix (B.3) we would show
B (B.4)
2x3
Remark B.l. Matrices should not be confused with determinants.3 A determinant is a number associated with square matrices (m = n), defined according to the rules stated in §B.5
As in the case of vectors, the components of a matrix may be real or complex. If they are real numbers, the matrix is called real, and complex otherwise. For the present exposition all matrices will be assumed real unless otherwise stated.
Of these two terms, "entry" is preferred when talking about the computer implementation, whereas "component" is preferred when connecting matrices to physical entities.
2 A convention introduced in Berkeley courses by Ray Clough. It is particularly convenient in blackboard expositions.
3 This confusion is apparent in the literature of the period 1860-1920.
B-3
Appendix B: LINEAR ALGEBRA: MATRICES
§B.1.1. Square Matrices
The case m = n is important in practical applications. Such matrices are called square matrices of order n. Matrices for which m ^ n are called non-square. The term "rectangular" is also used in this context, but this is slightly fuzzier because squares are special cases of rectangles.
Square matrices enjoy properties that are not shared by non-square matrices, such as the symmetry and antisymmetry conditions defined below. Furthermore many operations, such as taking determinants and computing eigenvalues, are only defined for square matrices.
Example B.2.
is a square matrix of order 3.
12 6 8 24 2 5
3' 7 11
(B.5)
Consider a square matrix A = [a^] of order n x n. Its n components an form the main diagonal, which runs from top left to bottom right. The cross diagonal or antidiagonal runs from bottom left to upper right. The main diagonal of the example matrix (B.5) is {12, 24, 11} whereas its cross diagonal is {2, 24, 3}.
Entries that run parallel to and above (below) the main diagonal form superdiagonals (subdiagonals). For example, {6, 7} is the first superdiagonal of the matrix (B.5).
The trace of an x n square matrix A, denoted trace(A), is the sum of its diagonal coefficients:
trace(A) = %•
(B.6)
§B.1.2. Symmetry and Antisymmetry
Square matrices for which a^- = a,-, are called symmetric about the main diagonal or simply symmetric. Square matrices for which atj = —ajt are called antisymmetric4. The diagonal entries of an antisymmetric matrix must be zero.
Example B.3. The following is a symmetric matrix of order 3:
(B.7)
	"11	6	1
S =	6	3	-1
	1	-1	-6
The following is an antisymmetric matrix of order 4:
W
0 -3
1 5
3 0 -7 2
-1 7 0 0
-5 -2 0 0
(B.8)
The alternative names skew-symmetric and antimetric are also used in the literature.
B-4
5B.1 MATRICES
§B.1.3. Special Matrices
The null matrix, written 0, is the matrix all of whose components are zero. Example B.4. The null matrix of order 2 x 3 is
TO 0 0" Lo 0 0.
(B.9)
The identity matrix, written I, is a square matrix all of which entries are zero except those on the main diagonal, which are ones.
Example B.5. The identity matrix of order 4 is
I =
10 0 0 0 10 0 0 0 10 0   0   0 1
(B.10)
A diagonal matrix is a square matrix all of which entries are zero except for those on the main diagonal, which may be arbitrary.
Example B.6. The following matrix of order 4 is diagonal:
D
14	0	0	0
0	-6	0	0
0	0	0	0
0	0	0	3
(B.ll)
A short hand notation which lists only the diagonal entries is sometimes used for diagonal matrices to save writing space. This notation is illustrated for the above matrix:
D = diag [ 14   - 6   0   3 ].
(B.12)
An upper triangular matrix is a square matrix in which all entries underneath the main diagonal vanish. A lower triangular matrix is a square matrix in which all entries above the main diagonal vanish. These forms are important in the solution of linear systems of equations.
Example B.7. Here are examples of each kind:
U =
6 4 2 1
0 6 4 2
0 0 6 4
0 0 0 6
L =
5	0	0	0
10	4	0	0
-3	21	6	0
-15	-2	18	7
(B.13)
An upper triangular matrix with all diagonal entries equal to 1 is called a unit upper triangular matrix. Similarly, a lower triangular matrix with all diagonal entries equal to 1 is called a unit upper triangular matrix.
B-5
Appendix B: LINEAR ALGEBRA: MATRICES
§B.1.4. * Are Vectors a Special Case of Matrices?
Consider the 3-vector x and a 3 x 1 matrix X with the same components:
	1 "		' xn ~
X	X2	X =	
	.x3.		-X3l _
(B.14)
in which x\ — xu, x2 — x22 and x3 = x33. Are x and X the same thing? If so we could treat column vectors as one-column matrices and dispense with the distinction. Indeed in many contexts a column vector of order n may be treated as a matrix with a single column, i.e., as a matrix of order nxl. Similarly, a row vector of order m may be sometimes treated as a matrix with a single row, i.e., as a matrix of order 1 x m.
There are operations, however, for which the analogy does not carry over, and one must consider vectors as different from matrices. The dichotomy is reflected in our notational conventions of lower versus upper case.
Another practical distinction between matrices and vectors can be cited: although we speak of "matrix algebra" as embodying vectors as special cases of matrices, in practice the quantities of primary interest to the engineer are vectors rather than matrices. For example, a structural engineer may be interested in displacement vectors, force vectors, vibration eigenvectors, buckling eigenvectors, etc. In finite element analysis even stresses and strains are often arranged as vectors5 although they are actually tensors.
On the other hand, matrices are rarely the quantities of primary interest: they work silently in the background where they are normally engaged in operating on vectors.
§B.2. Elementary Matrix Operations §B.2.1. Equality
Two matrices A and B of same order m x n are said to be equal if and only if all of their components are equal: a{- = b-, for all i = 1, ... m, j = 1, ... n. We then write A = B. If the inequality test fails for at least one entry, the matrices are said to be unequal and we write A^B.
Two matrices of different order cannot be compared for equality or inequality.
There are specialized tests for greater-than or less-than, but these are not considered here.
§B.2.2. Transposition
The transpose of am x n real matrix A = [a,7] is an x m matrix denoted by Ar, which satifies
(B.15)
The rows of Ar are the columns of A, and the rows of A are the columns of Ar. Obviously the transpose of Ar is again A, that is, (Ar)r = A.
Example B.8.
A- P   7 °" Ll   0 4.
AT
(B.16)
The transpose of a square matrix is also a square matrix. The transpose of a symmetric matrix A is equal to the original matrix, that is A = Ar. The negated transpose of an antisymmetric matrix matrix A is equal to the original matrix, that is A = —Ar.
The equivalent operation for a complex matrix involves conjugation and is not considered here.
5 Through the "casting" operation described in Chapter 11.
B-6
Example B.9.
A =
'4   7 0' 7   1 2 0   2 3
= AT,
W =
0 7 -7 0 0 2
5B.2.3. Addition, Subtraction and Scaling
5B.3 MATRIX PRODUCTS
= -W
(B.17)
The simplest operation acting on two matrices is addition. The sum of two matrices of the same order, A and B, is written A + B and defined to be the matrix
def r
B = [a,:
(B.18)
Like vector addition, matrix addition is commutative: A+B = B+A, and associative: A+(B + C) = (A + B) + C. For n = 1 or m = 1 the operation reduces to the addition of two column or row vectors, respectively.
For matrix subtraction, replace + by — in the definition (B.18).
Example B.10. The sum of
-3 0 4 2-1
[I
and
is   A +
»=[il I
Multiplication of a matrix A by a scalar c is defined by means of the relation
*   def r n
cA = [cau]
(B.19)
(B.20)
That is, each entry of the matrix is multiplied by c. This operation is called scaling of a matrix. If c = 0, the result is the null matrix. Division of a matrix by a nonzero scalar c is equivalent to multiplication by (1/c).
Example B.ll.
0"
-1
3A
-[
3
12
-9 6
0"
-3
(B.21)
§B.3. Matrix Products §B.3.1. Matrix by Vector Product
Before describing the general matrix product of two matrices, let us consider the particular case in which the second matrix is a column vector. This so-called matrix-vector product merits special attention because it occurs frequently in the applications. Let A = [a-] be an m x n matrix, x = {jc .} a column vector of order n, and y = {y{} a column vector of order m. The matrix-vector product is symbolically written
y = Ax, (B.22)
to mean the linear transformation
n
yi = ^2 aijxj = aijxj>       i = 1, ... ,m. (B.23)
SC
The symbol X may be omitted if using the summation convention, as ilustrated above after = .
B-7
Appendix B: LINEAR ALGEBRA: MATRICES
Example B.12. The product of a 2 x 3 matrix and a 3-vector is a 2-vector:
1 -3 0" 4 2-1.
This product definition is not arbitrary but emanates from the analytical and geometric properties of entities represented by matrices and vectors.
For the product definition to make sense, the column dimension of the matrix A (called the pre-multiplicand) must equal the dimension of the vector x (called the post-multiplicand). For example, the reverse product x A would not make sense unless m = n = 1.
If the row dimension m of A is one, the matrix formally reduces to a row vector, and the matrix-vector product reduces to the inner product defined by equation (A. 15) of Appendix A. The result of this operation is a one-dimensional vector or scalar. Thus we see that the present definition properly embodies previous cases.
The associative and commutative properties of the matrix-vector product fall under the rules of the more general matrix-matrix product discussed next.
§B.3.2. Matrix by Matrix Product and Matrix Powers
We now pass to the most general matrix-by-matrix product, and consider the operations involved in computing the product C of two matrices A and B:
C = AB. (B.25)
Here A = [a^] is a matrix of order m x n, B = [bjk] is a matrix of order n x p, and C = [cjk] is a matrix of order m x p. The entries of the result matrix C are defined by the formula
n
cik = J2 aUbJk = aijbjk>    i = l, ... ,m,    k=\, ... ,p. (B.26)
7 = 1
We see that the (i, k)th entry of C is computed by taking the inner product of the ith row of A with the kth column of B. For this definition to work and the product be possible, the column dimension of A. must be the same as the row dimension of~B. Matrices that satisfy this rule are said to be product-conforming, or conforming for short. If the two matrices do not conform, their product is undefined. The following mnemonic notation often helps in remembering this rule:
C = A    B (B.27)
mxp      mxn nxp
For the matrix-by-vector case treated in the preceding subsection, p = 1.
Matrix A is called the pre-multiplicand and is said to premultiply B. Matrix B is called the post-multiplicand and is said to postmultiply A. This careful distinction on which matrix comes first is a consequence of the absence of commutativity: even if B A exists (it only does if m = n), it is not generally the same as AB.
If A = B, the product A A is called the square of A and is denoted by A2. Note that for this definition to make sense, A must be a square matrix — else the factors would not be conforming. Similarly, A3 = AAA = A2A = AA2. Other positive-integer powers can be analogously defined.
[1x1 + (-3) x2 + 0x31_r-51 U x 1 +2 x 2+ (-1) x 3_| ~ L 5J
(B.24)
B-8
5B.3 MATRIX PRODUCTS
The above definition does not encompass negative powers. For inverse, A-1 denotes the inverse of A, which is studied later. The general power Am, where m can be a real or complex scalar, can be defined with the help of the matrix spectral form — an advanced topic covered in Appendix E.
A square matrix A that satisfies A = A2 is called idempotent. We shall see later that this condition characterizes the so-called projector matrices. A square matrix A whose pth power is the null matrix is called p-nilpotent.
Remark B.2. For hand computations, the matrix product is most conveniently organized by the so-called Falk's scheme:
nk
(B.28)
Each entry in row i of A is multiplied by the corresponding entry in column k of B (note the arrows), and the products are summed and stored in the (i, k)th entry of C.
Example B.13. To illustrate Falk's scheme, let us form the product C = AB of the following matrices
'2   1 0-5'
"3
0 -1
B
4 3 0 1
(B.29)
The matrices are conforming because the column dimension of A and the row dimension of B are the same (3). We arrange the computations as shown below:
2 1 0-5"
4 3-10 =B
.0 1 -7     4 J                  . (B.30)
T3     0   21    T6 5 -14   -71 _C_AR
A"   |_4   -1   5 J    U 6 -34     oj CAB Here 3x2 + 0x4 + 2x0 = 6 and so on.
§B.3.3. Matrix Product Properties
Associativity. The associative law is verified:
A(BC) = (AB)C. (B.31)
Hence we may delete the parentheses and simply write ABC.
Distributivity. The distributive law also holds: If B and C are matrices of the same order, then
A(B + C) = AB + AC,       and      (B + C) A = BA + C A. (B.32)
Commutativity. The commutativity law of scalar multiplication does not generally hold. If A and B are square matrices of the same order, then the products A B and B A are both possible but in general AB^BA. If A B = B A, the matrices A and B are said to commute. One important case is when A and B are diagonal. In general A and B commute if they share the same eigensystem.
B-9
Appendix B: LINEAR ALGEBRA: MATRICES
Example B.14. Matrices
a - ß b
[l   cl »=['
(B.33)
b      c - p.
commute for any a, b, c, p. More generally, A and B = A — p I commute for any square matrix A.
Transpose of a Product. The transpose of a matrix product is equal to the product of the transposes of the operands taken in reverse order:
(AB)r=BrAr. (B.34) The general transposition formula for an arbitrary product sequence is
(ABC...MN)r =NrMr...CrBrAr. (B.35)
Congruent Transformation. If B is a symmetric matrix of order m and A is an arbitrary m x n matrix, then
S = ArB A. (B.36)
is a symmetric matrix of order n. Such an operation is called a congruent or congruential transformation. It is often found in finite element analysis when changing coordinate bases because such a transformation preserves certain key properties (such as symmetry).
Loss of Symmetry. The product of two symmetric matrices is not generally symmetric.
Null Matrices may have Non-null Divisors. The matrix product AB can be zero although A ^ 0 and B^O. Likewise it is possible that A ^ 0, A2 ^ 0, ..., but Ap = 0.
§B.4. Bilinear and Quadratic Forms
Let x and y be two column vectors of order n, and A a real square n x n matrix. Then the following triple product produces a scalar result:
s = yT   A   x (B.37)
Ixn nxn nxl
This is called a bilinear form. Matrix A is called the kernel of the form.
Transposing both sides of (B.37) and noting that the transpose of a scalar does not change, we obtain the result
s = xrAry = yr Ax. (B.38)
If A is symmetric and vectors x and y coalesce, i.e. Ar = A and x = y, the bilinear form becomes a quadratic form
s = xTAx. (B.39) Transposing both sides of a quadratic form reproduces the same equation.
B-10
§B.5 DETERMINANTS
Example B.15. The kinetic energy of a dynamic system consisting of three point masses ml, m2, m3 moving in one dimension with velocities v\, v2 and w3, respectively, is
T = )^{mxv\ + in > r.' + m3vf).
This can be expressed as the quadratic form
in which
		T	1 T = -V 2 v
	ml	0	0 "
M =	0	m2	0
	. 0	0	m3 _
V =
(B.40)
(B.41)
(B.42)
Here M denotes the system mass matrix whereas v is the system velocity vector.
The following material discusses more specialized properties of matrices, such as determinants, inverses, rank and orthogonality. These apply only to square matrices unless othewise stated.
§B.5. Determinants
The determinant of a square matrix A = [a^] is a number denoted by |A| or det(A), through which important properties such as matrix singularity or its eigenvalue spectrum can be compactly characterized. This number is defined as the following function of the matrix elements:
|A| = det(A) = ± Y[ aljx a.
2j2 • • • anj„>
(B.43)
where the column indices j\, J2, ... jn are taken from the set {1,2,... n}, with no repetitions allowed. The plus (minus) sign is taken if the permutation (ji J2 ... jn) is even (odd).
Example B.16. If the order n does not exceed 3, closed forms of the determinant are practical. For a 2 x 2 matrix,
a2l a22
Ml "22 "12"21-
(B.44)
For a 3 x 3 matrix,
a\\ Cl\2 «13 a2\ a22 a23 Ö31     Ö32 Ö33
— ^11^22^33   ^~ ^12^23^31   ^~ ^13^21^32     ^13^22^31      ^12^21^33 Cl^^d^^d^^.
(B.45)
Remark B.3. The concept of determinant is not applicable to non-square matrices or to vectors. Thus the notation |x| for a vector x can be reserved for its magnitude (as in Appendix A) without risk of confusion.
Remark B.4. Inasmuch as the product (B.43) contains n\ terms, the calculation of |A| from the definition is impractical for general matrices whose order exceeds 3 or 4. For example, if n = 10, the product (B.43) contains 10! = 3, 628, 800 terms, each involving 9 multiplications, so over 30 million floating-point operations would be required to evaluate |A| according to that definition. A more practical method based on matrix decomposition is described in Remark B.3.
B-ll
Appendix B: LINEAR ALGEBRA: MATRICES §B.5.1. Determinant Properties
Some useful rules associated with the calculus of determinants are listed next.
I.      Rows and columns can be interchanged without affecting the value of a determinant. Consequently
|A| = |Ar|. (B.46)
II.     If two rows, or two columns, are interchanged the sign of the determinant is reversed. For example:
1 A 1 O
(B.47)
3 4		1	-2
1 -2		3	4
III.    If a row (or column) is changed by adding to or subtracting from its elements the corresponding elements of any other row (or column) the determinant remains unaltered. For example:
3 4		3+1 4-2		4	2
1 -2		1 -2		1	-2
-10.
(B.48)
IV.    If the elements in any row (or column) have a common factor a then the determinant equals the determinant of the corresponding matrix in which a = 1, multiplied by a. For example:
6 8	= 2	3	4
1 -2		1	-2
2 x (-10)
-20.
(B.49)
V. When at least one row (or column) of a matrix is a linear combination of the other rows (or columns) the determinant is zero. Conversely, if the determinant is zero, then at least one row and one column are linearly dependent on the other rows and columns, respectively. For example, consider
3     2 1
1      2   -1 . (B.50)
2-1 3
This determinant is zero because the first column is a linear combination of the second and third columns:
column 1 = column 2 + column 3. (B.51) Similarly, there is a linear dependence between the rows, which is given by the relation
row 1 = | row 2 + | row 3. (B.52)
VI.    The determinant of an upper triangular or lower triangular matrix is the product of the main diagonal entries. For example,
3	2	1		
0	2	-1	= 3 x 2 x 4 = 24.	(B.53)
0	0	4		
B-12
§B.5 DETERMINANTS
This rule is easily verified from the definition (B.43) because all terms vanish except j\ = 1, j2 = 2, ... jn = n, which is the product of the main diagonal entries. Diagonal matrices are a particular case of this rule.
VII.  The determinant of the product of two square matrices is the product of the individual determinants:
|AB| = |A||B|. (B.54)
The proof requires the concept of triangular decomposition, which is covered in the Remark below. This rule can be generalized to any number of factors. One immediate application is to matrix powers: |A2| = |A||A| = |A|2, and more generally |A™| = |A|™ for integers.
VIII. The determinant of the transpose of a matrix is the same as that of the original matrix:
IAI.
(B.55)
This rule can be directly verified from the definition of determinant, and also as direct consequence of Rule I.
Remark B.5. Rules VI and VII are the key to the practical evaluation of determinants. Any square nonsingular matrix A (where the qualifier "nonsingular" is explained in §B.3) can be decomposed as the product of two triangular factors
A = LU, (B.56)
in which L is unit lower triangular and U is upper triangular. This is called a LU triangularization, LU factorization or LU decomposition. It can be carried out in 0(n3) floating point operations. According to rule VII:
|A| = |L| |U|. (B.57)
According to rule VI, |L| = 1 and |U| = unu22 ... unn. The last operation requires only 0(n) operations. Thus the evaluation of |A| is dominated by the effort involved in computing the factorization (B.56). For n = 10, that effort is approximately 103 = 1000 floating-point operations, compared to approximately 3 x 107 from the naive application of the definition (B.43), as noted in Remark B.4. Thus the LU-based method is roughly 30, 000 times faster for that modest matrix order, and the ratio increases exponentially for large n.
§B.5.2. Cramer's Rule and Homogeneous Systems
Cramer's rule provides a recipe for solving linear algebraic equations directly in terms of determinants. Let the system of equations be — as usual — denoted by
Ax = y,
(B.58)
in which A is a given n x n matrix, y is a given hxI vector, and x is the «xl vector of unknowns. The explicit form of (B.58) is Equation (A.l) of Appendix A, with n = m.
A closed form solution for the components x\, X2 ■ ■., xn of x in terms of determinants is
yi an «13 yi «23
yn     ßfi2 ßfi3
an   yi an
Ü2\ yi «23 fl-nl     yn CLrii
O-ln O-ln
(B.59)
B-13
Appendix B: LINEAR ALGEBRA: MATRICES
The rule can be remembered as follows: in the numerator of the quotient for x ■, replace the jth column of A by the right-hand side y.
This method of solving simultaneous equations is known as Cramer's rule. Because the explicit computation of determinants is impractical for n > 3 as explained in Remark C.3, direct use of the rule has practical value only for n = 2 and n = 3 (it is marginal for n = 4). But such small-order systems arise often in finite element calculations at the Gauss point level; consequently implementors should be aware of this rule for such applications. In addition the rule may be advantageous for symbolic systems; see Example below.
Example B.17. Solve the 3x3 linear system
by Cramer's rule:
Xi =
8 2 1 5 2 0 3   0 2
5 2 1 3 2 0 1   0 2
6
*2
~5	2	1		X\		8
3	2	0			—	5
_1	0	2		x%		3
5 8 1 3 5 0 1   3 2
5 2 1 3 2 0 1   0 2
-5-1.
6
x3
5 2 8 3 2 5 1   0 3
5 2 1 3 2 0 1   0 2
(B.60)
6
(B.61)
Example B.18. Solve the 2x2 linear algebraic system
L -P   i + p\Yx2\ LoJ
by Cramer's rule:
	5 -ß 0   l + ß		5 + 5ß	x2 = -		2 + ß 5 -ß 0	5ß	
2 + ß -ß -ß l+ß			2 + 3ß'		2 + ß -ß -ß     l + ß			2 + 3ß
(B.63)
Remark B.6. Creamer's rule importance has grown in symbolic computations carried out by computer algebra systems. This happens when the entries of A and y are algebraic expressions. For example the example system (B.62). In such cases Cramer's rule may be competitive with factorization methods for up to moderate matrix orders, for example n < 20. The reason is that determinantal products may be simplified on the fly.
One immediate consequence of Cramer's rule is what happens if
yi = y2 = .. . = yn = 0. (B.64)
A linear equation system with a null right hand side
Ax = 0, (B.65)
is called a homogeneous system. From the rule (B.59) we see that if |A| is nonzero, all solution components are zero, and consequently the only possible solution is the trivial one x = 0. The case in which |A| vanishes is discussed in the next section.
B-14
§B.6 SINGULAR MATRICES, RANK
§B.6. Singular Matrices, Rank
If the determinant | A| of a n x n square matrix A = An is zero, then the matrix is said to be singular. This means that at least one row and one column are linearly dependent on the others. If this row and column are removed, we are left with another matrix, say A„_i, to which we can apply the same criterion. If the determinant |A„_i | is zero, we can remove another row and column from it to get A„_2, and so on. Suppose that we eventually arrive at an r x r matrix Ar whose determinant is nonzero. Then matrix A is said to have rank r, and we write rank(A) = r.
If the determinant of A is nonzero, then A is said to be nonsingular. The rank of a nonsingular n x n matrix is equal to n.
Obviously the rank of Ar is the same as that of A since it is only necessary to transpose "row" and "column" in the definition.
The notion of rank can be extended to rectangular matrices; as covered in the textbooks cited in Notes and Bibliography. That extension, however, is not important for the material covered here.
Example B.19. The 3x3 matrix
has rank r = 3 because I Al
Example B.20. The matrix
-3 / 0.
	"3	2	2"
A =	1	2	-1
	_2	-1	3_
	"3	2	1"
A =	1	2	-1
	_2	-1	3_
(B.66)
(B.67)
already used as example in §B.5.1, is singular because its first row and column may be expressed as linear combinations of the others as shown in (B.51) and (B.52). Removing the first row and column we are left with a 2 x 2 matrix whose determinant is 2 x 3 — (—1) x (— 1) = 5 ^ 0. Thus (B.67) has rank r = 2.
§B.6.1. Rank Deficiency
If the square matrix A is supposed to be of rank r but in fact has a smaller rank f < r, the matrix is said to be rank deficient. The number r — f > 0 is called the rank deficiency.
Example B.21. Suppose that the unconstrained master stiffness matrix K of a finite element has order n, and that the element possesses b independent rigid body modes. Then the expected rank of K is r = n — b. If the actual rank is less than r, the finite element model is said to be rank-deficient. This is usually undesirable.
Example B.22. An an illustration of the foregoing rule, consider the two-node, 4-DOF, Bernoulli-Euler plane beam element stiffness derived in Chapter 12.
(B.68)
in which EI and L are positive scalars. It may be verified (for example, through an eigenvalue analysis) that this 4x4 matrix has rank 2. The number of rigid body modes is 2, and the expected rank is r = 4 — 2 = 2. Consequently this FEM model is rank sufficient.
	12 6L	-12	6L
EI	4L2	-61	2L2
U		12	-6L
	symm		4L2
B-15
Appendix B: LINEAR ALGEBRA: MATRICES §B.6.2. Rank of Matrix Sums and Products
In finite element analysis matrices are often built through sum and product combinations of simpler matrices. Two important rules apply to "rank propagation" through those combinations.
The rank of the product of two square matrices A and B cannot exceed the smallest rank of the multiplicand matrices. That is, if the rank of A is ra and the rank of B is rb,
rank(AB) < min(ra, rb). (B.69)
In particular, if either matrix is singular, so is the product. Regarding sums: the rank of a matrix sum cannot exceed the sum of ranks of the summand matrices. That is, if the rank of A is ra and the rank of B is rb,
rank(A + B) < ra + rb. (B.70)
Rule (B.70) is the basis for studying the effect of the number of Gauss points on the rank of an element stiffness matrix being built through numerical integration, as discussed in Chapter 19.
§B.6.3. Singular Systems: Particular and Homogeneous Solutions
Having introduced the notion of rank we can now discuss what happens to the linear system (B.58) when the determinant of A vanishes, meaning that its rank is less than n. If so, (B.58) has either no solution or an infinite number of solutions. Cramer's rule is of limited or no help in this situation.
To discuss this case further we note that if |A| =0 and the rank of A is r = n — d, where d > 1 is the rank deficiency, then there exist d nonzero independent vectors z,, i = 1, ... d such that
Az,=0. (B.71)
These d vectors, suitably orthonormalized, are called null eigenvectors of A, and form a basis for its null space.
Let Z denote the n x d matrix obtained by collecting the z, as columns. If y in (B.58) is in the range of A, that is, there exists an nonzero xp such that y = Axp, its general solution is
x = Xp + xA = xp + Zw, (B.72)
where w is an arbitrary d x 1 weighting vector. This statement can be easily verified by substituting this solution into Ax = y and noting that AZ vanishes.
The components xp and Xk are called the particular and homogeneous portions respectively, of the total solution x. (The terminology: particular solution and particular solution, are often used.) If y = 0 only the homogeneous portion remains.
If y is not in the range of A, system (B.58) does not generally have a solution in the conventional sense, although least-square solutions can usually be constructed. The reader is referred to the many textbooks in linear algebra for further details.
B-16
§B.7 MATRIX INVERSION
§B.6.4. *Rank of Rectangular Matrices
The notion of rank can be extended to rectangular matrices, real or complex, as follows. Let A be m x n. Its column range space 1Z(A) is the subspace spanned by Ax where x is the set of all complex n-vectors. Mathematically: 1Z(A) = {Ax : x e C"}. The rank r of A is the dimension of 1Z(A).
The null space J\f (A) of A is the set of n-vectors z such that A z = 0. The dimension of J\f (A) is n — r.
Using these definitions, the product and sum rules (B.69) and (B.70) generalize to the case of rectangular (but conforming) A and B. So does the treatment of linear equation systems A x = y in which A is rectangular. Such systems often arise in the fitting of observation and measurement data.
In finite element methods, rectangular matrices appear in change of basis through congruential transformations, and in the treatment of multifreedom constraints.
§B.7. Matrix Inversion
The inverse of a square nonsingular matrix A is represented by the symbol A-1 and is defined by the relation
AA '=1. (B.73)
The most important application of the concept of inverse is the solution of linear systems. Suppose that, in the usual notation, we have
Ax = y. (B.74) Premultiplying both sides by A-1 we get the inverse relationship
x = A~1y. (B.75)
More generally, consider the matrix equation for multiple (m) right-hand sides:
A   X = Y , (B.76)
nxn nxm nxm
X = A_1Y. (B.77)
which reduces to (B.74) for m = 1. The inverse relation that gives X as function of Y is In particular, the solution of
AX = I, (B.78)
is X = A-1. Practical methods for computing inverses are based on directly solving this equation; see Remark below.
§B.7.1. Explicit Computation of Inverses
The explicit calculation of matrix inverses is seldom needed in large matrix computations. But occasionally the need arises for the explicit inverse of small matrices that appear in element level computations. For example, the inversion of Jacobian matrices at Gauss points, or of constitutive matrices.
A general formula for elements of the inverse can be obtained by specializing Cramer's rule to
"» = M- <B'79)
(B.78). LetB = [by] = A-1. Then
B-17
Appendix B: LINEAR ALGEBRA: MATRICES
in which Aji denotes the so-called adjoint of entry aij of A. The adjoint A7, is defined as the
th
determinant of the submatrix of order (n — 1) x (n — 1) obtained by deleting the j   row and i column of A, multiplied by (—1)'+7.
This direct inversion procedure is useful only for small matrix orders, say 2 or 3. In the examples below the explicit inversion formulas for second and third order matrices are listed.
Example B.23. For order n — 2:
la2i
in which | A| is given by (B.44). Example B.24. For order n = 3:
"12
a2i
where
	an	an	a13
A =	a2\	a21	a23
	_a3l	a32	a33-
bn =	ail	a23	h\
	an	ÍZ33	
bn —
bn —
a.21 Ü23
a3\ a33
«21 a22
ÍZ31 ÍZ32
in which IAI is given by (B.45) Example B.25.
A =
3 1 2 3 -1 2
£22 =
b23 =
0 5
_ _L ra
~ |A| l-t
22 a2i
A"1 =
1
TAT
Cl\2 ÍZ13
a32 a33
an a3l
au a31
a33
an
032
A"1 = -
1
36
"21 b 31
hi
h2 Z733
15 -10 -1
-an
a22
-3 14 5
bn	hn	
b22	b23	
b32	b33-	
	a\2	an
	a22	a23
	a\\	an
	a2i	a23
	an	an
	a.21	a22
(B.80)
(B.81)
(B.82)
(B.83)
If the order exceeds 3, the general inversion formula based on Cramer's rule becomes rapidly useless because it displays combinatorial complexity as noted in a previous Remark. For numerical work it is preferable to solve (B.78) after A is factored. Those techniques are described in detail in linear algebra books; see also Remark below.
§B.7.2. Some Properties of the Inverse
I. If A is nonsingular, the determinant of A-1 is the inverse of the determinant of A
II. If A-1 exists, the inverse of its transpose is equal to the transpose of the inverse:
(Ar)_1 = (A_1)r, (B.84)
because (AA"1) = (AA_1)r = (A_1)r Ar = I.
III. The inverse of a symmetric matrix is also symmetric. Because of the previous rule, (Ar) ~1 = A-1 = (A_1)r, hence A-1 is also symmetric.
B-18
§B.8 CHARACTERISTIC POLYNOMIAL AND EIGENVALUES
IV. The inverse of a matrix product is the reverse product of the inverses of the factors:
(AB)_1 =B_1A_1. (B.85)
which is easily proven by substitution. This rule generalizes to an arbitrary number of factors.
V. For a diagonal matrix D in which all diagonal entries are nonzero, D 1 is again a diagonal matrix with entries 1 /dn. The verification is straightforward.
VI. If S is a block diagonal matrix:
	rSn	0	0	o -
	0		0	0
s =	0	0	$33 •	0
	_ 0	0	0	•   &nn -
= diag[SH]
(B.86)
in which all S77 are nonsingular, the inverse matrix is also block diagonal and is given by
0
0      S^1 0
0      0 S331
L 0      0 0
0 0 0
= diag[S//1]
(B.87)
VII. The inverse of an upper triangular matrix is also an upper triangular matrix. The inverse of a lower triangular matrix is also a lower triangular matrix. Both inverses can be computed in 0(n2) floating-point operations.
VIII. An antisymmetric matrix is always singular.
Remark B.7. The practical numerical calculation of inverses is based on triangular factorization.6 Given
a nonsingular n x n matrix A, calculate its LU factorization A = LU, which can be obtained in 0(n3) operations. Then solve the linear triangular systems:
UY = I,
LX = Y,
(B.88)
and the computed inverse A~ appears in X. One can overwrite I with Y and Y with X. The whole process can be completed in 0(n3) floating-point operations. For symmetric matrices the alternative decomposition A = LDLT, where L is unit lower triangular and D is diagonal, is generally preferred to save computing time and storage.
See References at the end of the Appendix.
B-19
Appendix B: LINEAR ALGEBRA: MATRICES
§B.8.  Characteristic Polynomial and Eigenvalues
Let A be a square n x n matrix, and X a (real or complex) scalar variable. The determinant of A — XI is a polynomial in X of order n:
P(X) = det(A - X I) = C0 + CiX + ... + C„ Xn, (B.89)
which is called the charactristic polynomial of A. The C, are the characteristic coefficients. It is easy shown that Co = det(A) and C„_i = (—1)"_1 trace(A). The n roots of P(X) = 0 are the eigenvalues of A. These may be complex even if A is real. This topic is developed in Appendix E.
Example B.26. The characteristic polynomial of the matrix A in (B.83) is
P(X) = det(A -Xh)
3 - X      1 -1 2      3 — X 0 -1       2 5-X
36 - 36A + llA2-A3. (B.90)
The roots of P(X) — 0, in ascending order, are X\ — 2, X2 — 3 and X3 — 6. These are the eigenvalues of A. It may be verified that the eigenvalues of A-1 are 1/2, 1/3 and 1/6. This illustrates a general rule: the eigenvalues of a nonsingular A and its inverse A-1 are reciprocal, as further discussed in Appendix E.
§B.8.1. *Similarity Transformations and Invariants
Let again A be a square n x n matrix, whereas T denotes a nonsingular n x n matrix. The transformation
B = T 'AT, (B.91)
is called a similarity transformation. A and B are said to be similar matrices. The determinants of A and B are identical since |B| = |T| | A| (1/|T|) — see Rule VII in §B.5.1. Likewise the determinants of A — XI and B — XI stay the same. This implies that the characteristic polynomials of A and B are identical, whence their eigenvalues coincide.
A matrix invariant or simply invariant is any expression that does not change under a similarity transformation. The characteristic coefficients furnish an obvious example since P(X) is unchanged. In particular Co = trace(A) and C„ = det(A) are invariants. The n coefficients C„ of P(X) are said to form an invariant basis. Any expression that is only a function of those coefficients, such as the eigenvalues, is also an invariant.
§B.8.2. * Congruential Transformations
A matrix transformation that often appears in FEM developments is the congruential transformation (also called congruent):
B = TTAT. (B.92)
Matrices B and A related through (B.92) are said to be congruent. This operation becomes a similarity transformation if TT = T-1, or TTT = I, which for real A characterizes an orthogonal matrix.
Aside from the special case of T being orthogonal, a congruential transformation does not preserve eigenvalues (or generally, any invariants). But (B.92) can be generalized to product-conforming rectangular T matrices, whereas (B .91) cannot. Another key difference: if A is real symmetric, a congruential transformation preserves symmetry, whereas a similarity transformation with a non-orthogonal matrix destroys it.
B-20
§B.9 *THE INVERSE OF A SUM OF MATRICES
§B.8.3. *The Cayley-Hamilton Theorem
The Cayley-Hamilton theorem is of great importance in advanced matrix theory. It states that a square matrix satisfies its own characteristic equation. Symbolically: P(a) = 0.
Example B.27. To see the theorem in action, consider the matrix a given in (b.83). Its characteristic polynomial was computed in (b.90) to be P{X)36 - 36 X + 11 X2 - X3. Replacing Xm by am (if m = 0, by i3) we can check that
,0   0 0"
P(a) = 36i3 - 36 a + 11 a2 - a3
0 0 0 0   0 0
0.
(b.93)
§B.9. *The Inverse of a Sum of Matrices
The formula for the inverse of a matrix product: (ab)-1 = b 1 a-1 is not too different from its scalar counterpart: (ab)~l = (l/a)(l/b) = (l/b)(l/a), except that factor order matters. On the other hand, formulas for matrix sum inverses in terms of the summands are considerably more involved, and there are many variants. We consider here the expression of (a + b)-1 where both a and b are square and a is nonsingular. We begin from the identity introduced by Henderson and Searle in their review article [385]:
(i+p)-1 = (i + pr1 (i + p - p) = i - (i + pr1 p, (b.94)
in which p is a square matrix. Using (b.94) we may develop (a + b)-1 as follows:
(a + Br1 = (a(i + a"1 Br1) 1 = (i + a"1 Br1 a *b
= (i - (i + a"1 Br1 a"1) a"1 = a"1 - (i + a"1 Br1 a"1 b a"
(b.95)
Here b may be singular. If b = 0, it reduces to a = a . The check b = f> a also works. The last expression in (b.95) may be further transformed by matrix manipulations as
(a + Br1 = a"1 - (i + a"1 Br1 a"1 b a"1 = a"1 - a '(i + ba ') 'ba 1 = a"1 - a *b (i + a^b)-1 a"1
= a^1 - a-'ba-'ci + ba-1)-1.
In all of these forms b may be singular (or null). If b is also invertible, the third expresion in (b.96) may be transformed to the a commonly used variant
(a + Br1 = a"1 - a"1 (a-1 + b-'r1 a"1. (b.97)
The case of singular a may be handled using the notion of generalized inverses. This is a topic beyond the scope of this course, which may be studied, e.g., in the textbooks [83,129,686]. The special case of b being of low rank merges with the Sherman-Morrison and Woodbury formulas, covered below.
The Sherman-Morrison formula discussed below gives the inverse of a matrix modified by a rank-one matrix. The Woodbury formula extends it to a modification of arbitrary rank. In structural analysis these formulas are of interest for problems of structural modifications, in which a finite-element (or, in general, a discrete model) is changed by an amount expressable as a low-rank correction to the original model.
B-21
Appendix B: LINEAR ALGEBRA: MATRICES
§B.9.1. *The Sherman-Morrison and Woodbury Formulas
Let A be a square n x n invertible matrix, whereas u and v are two n-vectors and /3 an arbitrary scalar. Assume that a = 1 + p vTA 'u / 0. Then
(A + ySuv^r1 = A"1 - - A 'uvT A-1. (B.98)
When /3 = 1 this is called the Sherman-Morrison formula after [740]. (For a history of this remarkable expression and its extensions, which are important in many applications such as statistics and probability, see the review paper by Henderson and Searle cited previously.) Since any rank-one correction to A can be written as /3uvT, (B.98) gives the rank-one change to its inverse. The proof can be done by direct multiplication; see Exercise B.15.
For practical computation of the change one solves the linear systems Aa = u and Ab = v for a and b, using the known A-1. Compute a — 1 + /3vTa. If a ^ 0, the change to A-1 is the dyadic — {fi/a)ahT.
Consider next the case where U and V are two n x k matrices with k < n and /3 an arbitrary scalar. Assume that the k x k matrix £ — Ik + /3 VTA~'U, in which 1^ denotes the k x k identity matrix, is invertible. Then
(A + y6UVr)_1 = A"1 - /? \~lVITl\TArl. (B.99)
This is called the Woodbury formula, after [890]. It reduces to the Sherman-Morrison formula (B.98) if k = 1, in which case U = a is a scalar.
§B.9.2. * Formulas for Modified Determinants
Let A denote the adjoint of A. Taking the determinants from both sides of A + /3uvT one obtains
| A + /3uvT| = | A| +/3vTAu. (B.100) If A is invertible, replacing A = |A| A-1 this becomes
|A + /3uvT| = |A| (1 + /3v7A~1u). (B.101) Similarly, one can show that if A is invertible, and U and V are n x k matrices,
| A + /3UV7'| = |A| |It+ ^VTA-1U|. (B.102)
§B.9.3. * Inversion of Low Order Symmetric Matrix Sums
The inversion formulas given here are important in recent developments in continuum mechanics. Let A be an x n positive definite symmetric matrix with n = 2, 3,1 the identity matrix of order n and c a scalar. We want to express A + cl)~' directly in terms of A, I, c and the invariants of A provided by its characteristic coefficients. The derivations are carried out in [402], and make essential use of the Cayley-Hamilton theorem.
Two Dimensional Case. Here n = 2. The invariants are IA = trace(A) and IIa = det(A). Then
(A + cD-' =-- \    tt  (A-(c+/A)I). (B.103)
c (c + IA) + IIA
Three Dimensional Case. Here n — 3. The invariants are IA — trace(A), IIA — (trace(A2) — IA)/2, and IIIA = det(A). Then
(A + cir1 = 1 (A2 - (c + IA) A + [c(c + IA) + IIa\ D) ■ (B.104)
c[c(c + IA) + IIA] + HIa v '
B-22
§B.10 *MATRIX ORTHOGONALITY
§B.10. *Matrix Orthogonality
Let A and B be two product-conforming real matrices. For example, A is k x m whereas B is m x n. If their product is the null matrix
C = AB = 0, (B.105)
the matrices are said to be orthogonal. This is the generalization of the notions of vector orthogonality discussed in the previous Appendix.
§B.10.1. *Matrix Orthogonalization Via Projectors
The matrix orthogonalization problem can be stated as follows. Product conforming matrices A and B are given but their product is not zero. How can A be orthogonalized with respect to B so that (B.105) is verified? Suppose that B is m x n with m >n and that BTB is nonsingular (equivalently, B has full rank).7 Then form the m x m orthogonal projector matrix, or simply projector
1 TlT
PB = I-B(B' B) 'B
(B.106)
in which I is the m x m identity matrix. Since PB = PB, the projector is square symmetric. Note that
PBB = B-B(BTB) 1 (BTB) = B - B = 0. (B.107) It follows that PB projects B onto its null space. Likewise BT PB = 0. Postmultiplying A by PB yields
A = APB = A - AB (BTB) 1 BT.
(B.108)
Matrix A is called the projection of A onto the null space of B.8 It is easily verified that A and B are orthogonal:
AB = AB - AB (BTB) 1 (BTB) = AB - AB = 0. (B.109)
Consequently, forming A via (B.106) and (B.108) solves the orthogonalization problem.
If B is square and nonsingular, A = 0, as may be expected. If B has more columns than rows, that is m < n, the projector (B.106) cannot be constructed since BBT is necessarily singular. A similar difficulty arises if m > n but BTB is singular. Such cases require treatment using generalized inverses, which is a topic beyond the scope of this Appendix.9
In some applications, notably FEM, matrix A is square symmetric and it is desirable to preserve symmetry in A. That can be done by pre-and postmultiplying by the projector:
A = P B A P B.
(B.110)
Since     = PB, the operation (B. 110) is a congruential transformation, which preserves symmetry.
7 If you are not sure what "singular", "nonsingular" and "rank" mean or what (.) 1 stands for, please read §B.7.
8 In contexts such as control and signal processing,     is called a filter and the operation (B.108) is called filtering.
9 See e.g., the textbooks [83,686].
B-23
Appendix B: LINEAR ALGEBRA: MATRICES §B.10.2. *Orthogonal Projector Properties
The following properties of the projector (B.106) are useful when checking out computations. Forming its square as
p2 = pB pB = i _ 2B (BTB) lBT + B (BTB) lBT B (BTB) lBT
= I-2B(BTB) 'B^ + B^B) 'BT = I-B(BtB) 'Bt = Pb, (B.lll)
shows that the projector matrix is idempotent. Repeating the process one sees that p'^ = PB, in which n is an arbitrary nonnegative integer.
If B is m xn with m >n and full rank n, PB has m — n unit eigenvalues and n zero eigenvalues. This is shown in the paper [258], in which various applications of orthogonal projectors and orthogonalization to multilevel FEM computations are covered in detail.
Notes and Bibliography
Much of the material summarized here is available in expanded form in linear algebra textbooks. Highly reputable ones include Bellman [73], Golub and Van Loan [335], Householder [410], and Strang [773]. As books of historical importance we may cite Aitken [9], Muir [555], Turnbull [833], and Wilkinson and Reinsch [876]. For references focused on the eigenvalue problem, see Appendix E.
For inverses of matrix sums, there are two SIAM Review articles: [359,385]. For an historical account of the topic and its close relation to the Schur complement, see the bibliography in Appendix P.
B-24
Exercises
Homework Exercises for Appendix B: Matrices
In all ensuing exercises matrices and vectors are assumed real unless otherwise stated. EXERCISE B.l [A: 10] Given the three matrices
2 4 -1 2 2 5
1 0' 3 1 -1 2
	2	-2"
B =	1	0
	4	1
	-3	2
ri -3
1.2 0
(EB.l)
compute the product D = ABC by hand using Falk's scheme. Hint: do BC first, then premultiply that by A.
EXERCISE B.2 [A: 10] Given the square matrices
verify by direct computation that A B ^ B A. EXERCISE B.3 [A: 10] Given the matrices
3
-4 2
(EB.2)
	" 1	0"	
A =	-1	2	B =
	2	0	
-1 4' 2 0 0 0
(EB.3)
(note that B is symmetric) compute S = A B A, and verify that S is symmetric. EXERCISE B.4 [A: 10] Given the square matrices
A =
-1 0
-2
B =
-6 -14 2
(EB.4)
verify that A B = 0 although A / 0 and B / 0. Is B A also null? EXERCISE B.5 [A: 10] Given the square matrix
A =
0 a b 0 0c 0   0 0
(EB.5)
in which a,b,c are real or complex scalares, show by direct computation that A2    0 but A3 = 0. EXERCISE B.6 [A: 10] Can a diagonal matrix be antisymmetric?
EXERCISE B.7 [A:20] Prove the matrix product transposition rule (B.34). Hint: call C = (AB)T, D = BTAT, and use the matrix product definition (B.26) to show that the generic entries of C and D agree.
EXERCISE B.8 [A:20] If A is an arbitrary m x n real matrix, show: (a) both products ATA and AAT exist, and (b) both are square and symmetric. Hint: for (b) use the symmetry condition S = ST and (?).
EXERCISE B.9 [A: 15] Show that A2 only exists if and only if A is square.
B-25
Appendix B: LINEAR ALGEBRA: MATRICES
EXERCISE B.10 [A:20] If A is square and antisymmetric, show that A2 is symmetric. Hint: start from A = — AT and apply the results of Exercise B.8.
EXERCISE B.ll [A:15] If A is a square matrix of order n and c a scalar, show that det(cA) = c" det A.
EXERCISE B.12 [A:20] Let u and v denote real n-vectors normalized to unit length, so that uT u = 1 and yTy = 1, and let I denote the n x n identity matrix. Show that
det(I - uvT) = 1 - vTu (EB.6)
EXERCISE B.13 [A:25] Let u denote a real n-vector normalized to unit length so that uT u = 1, while I denotes the n x n identity matrix. Show that
H = I - 2uuT (EB.7)
is orthogonal: HTH = I, and idempotent: H2 = H. This matrix is called a elementary Hermitian, a Householder matrix, or a reflector. It is a fundamental ingredient of many linear algebra algorithms; for example the QR algorithm for finding eigenvalues.
EXERCISE B.14 [A:20] The trace of a n x n square matrix A, denoted trace(A) is defined by (B.6). Show that if the entries of A are real,
n n
trace(ATA) = ^ ^ a2 (EB.8)
i=\ j=i
This is the same as the Frobenius norm 11 A| | F (incorrectly called Euclidean norm by some authors).
EXERCISE B.15 [A:25] Prove the Sherman-Morrison formula (B.98) by direct matrix multiplication.
EXERCISE B.16 [A:30] Prove the Woodbury formula (B.99) by considering the following block bordered system
in which Ik and I„ denote the identity matrices of orders k and n, respectively. Solve (EB.9) two ways: eliminating first B and then C, and eliminating first C and then B. Equate the results for B.
EXERCISE B.17 [A:20] Show that the eigenvalues of a real symmetric square matrix are real, and that the eigenvectors are real vectors.
EXERCISE B.18 [A:25] Let the n real eigenvalues A, of a real n x n symmetric matrix A be classified into two subsets: r eigenvalues are nonzero whereas n — r are zero. Show that A has rank r.
EXERCISE B.19 [A: 15] Show that if A is positive definite, Ax = 0 implies that x = 0.
EXERCISE B.20 [A:20] Show that for any real m x n matrix A, ATA is nonnegative. (Existence is proven in a previous Exercise).
EXERCISE B.21 [A: 15] Show that a triangular matrix is normal if and only if it is diagonal.
EXERCISE B.22 [A:20] Let A be a real orthogonal matrix. Show that all of its eigenvalues A,, which are generally complex, have unit modulus.
EXERCISE B.23 [A: 15] Let A and B be two real symmetric matrices. If the product A B is symmetric, show that A and B commute.
B-26
Exercises
EXERCISE B.24 [A:30] Let A and T be real n x n matrices, with T nonsingular. Show that T 'AT and A have the same eigenvalues. (This is called a similarity transformation in linear algebra).
EXERCISE B.25 [A:40] (Research level) Let Abemxn and Bbenxm. Show that the nonzero eigenvalues of AB are the same as those of BA (Kahan).
EXERCISE B.26 [A:25] Let A be real skew-symmetric, that is, A = — AT. Show that all eigenvalues of A are either purely imaginary or zero.
EXERCISE B.27 [A:30] Let A be real skew-symmetric, that is, A = -AT. ShowthatU= (1 +A)_1(I —A), called a Cayley transformation, is orthogonal.
EXERCISE B.28 [A:30] Let P be a real square matrix that satisfies
P2 = P. (EB.10)
Such matrices are called idempotent, and also orthogonal projectors. Show that (a) all eigenvalues of P are either zero or one; (b) I — P is also a projector; (c) PT is also a projector.
EXERCISE B.29 [A:35] The necessary and sufficient condition for two square matrices to commute is that they have the same eigenvectors.
EXERCISE B.30 [A:30] A matrix whose elements are equal on any line parallel to the main diagonal is called a Toeplitz matrix. (They arise in finite difference or finite element discretizations of regular one-dimensional grids.) Show that if Ti and t2 are any two Toeplitz matrices, they commute: t1t2 = t2t1. Hint: do a Fourier transform to show that the eigenvectors of any Toeplitz matrix are of the form [e'mih}\ then apply the previous Exercise.
B-27