U

Fred E. Szabo PhD , in The Linear Algebra Survival Guide, 2015

Manipulation

Exploring upper-triangular matrices

Manipulate [MatrixForm [UpperTriangularize [{{1, 2 a, 3}, {4, 5, 6 b}, {7, 8c, 9}}]], {a, −   2, 2, 1}, {b, −   3, 3, 1}, {c, −   5, 5, 1}]

We use Manipulate, MatrixForm, and UpperTriangularize to construct and explore upper-triangular matrices. If we let a  =   1, b  =   2, and c  =   3, then the UpperTriangualize function converts the matrix

MatrixForm[{{1, 2, 3}, {4, 5, 12}, {7, 24, 9}}]

1 2 3 4 5 12 7 24 9

to the upper-triangular matrix

1 2 3 0 5 12 0 0 9

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B978012409520550028X

Determinants and Eigenvalues

Stephen Andrilli , David Hecker , in Elementary Linear Algebra (Fourth Edition), 2010

Highlights

The determinant of an upper (or lower) triangular matrix is the product of the main diagonal entries.

A row operation of type (I) involving multiplication by c multiplies the determinant by c.

A row operation of type (II) has no effect on the determinant.

A row operation of type (III) negates the determinant.

If an n × n matrix A is multiplied by c to produce B, then |B| = cn |A|.

The determinant of a matrix can be found by row reducing the matrix to upper triangular form and keeping track of the row operations performed and their effects on the determinant.

An n × n matrix A is nonsingular iff |A| ≠ 0 iff rank(A) = n.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780123747518000226

Algorithms

William Ford , in Numerical Linear Algebra with Applications, 2015

9.3 The Solution to Upper and Lower Triangular Systems

This section presents algorithms for solving upper- and lower-triangular systems of equations. In addition to providing additional algorithms for study, we will need to use both these algorithms throughout this book.

An upper-triangular matrix is an n  × n matrix whose only nonzero entries are below the main diagonal; in other words

a ij = 0 , j < i , 1 i , j n .

If U is an n  × n upper-triangular matrix, we know how to solve the linear system Ux  = b using back substitution. In fact, this is the final step in the Gaussian elimination algorithm that we discussed in Chapter 2. Compute the value of xn   = bn/unn , and then insert this value into equation (n    1) to solve for x n    1. Continue until you have found x 1. Algorithm 9.4 presents back substitution in pseudocode.

Algorithm 9.4

Solving an Upper Triangular System

function BACKSOLVE(U,b)

%   Find the solution to Ux  = b,   where U is an n  × n upper-triangular matrix.

xn   = bn/unn

for i   =   n-1:-1:1 do

  sum   =   0.0

for j   =   i+1:n do

sum  = sum  + uijxj

end for

x  (i)   =   (b  (i)   sum)   /uii

end for

return x

end function

NLALIB: The function backsolve implements Algorithm 9.4.

A lower-triangular matrix is a matrix all of whose elements above the main diagonal are 0; in other words

a ij = 0 , j > i , 1 i , j n .

A lower-triangular system is one with a lower-triangular coefficient matrix.

a 11 0 0 0 a 21 a 22 0 0 a 31 a 32 a 33 0 0 a n 1 a n 2 a n , n 1 a nn x 1 x 2 x 3 x n = b 1 b 2 b 3 b n

The solution to a lower-triangular system is just the reverse of the algorithm for solving an upper-triangular system—use forward substitution. Solve the first equation for x 1 = b 1 a 11 , and insert this value into the second equation to find x 2, and so forth.

Example 9.7

Solve

2 0 0 3 1 0 1 4 5 x 1 x 2 x 3 = 2 1 8

x 1 = 2 / 2 = 1 3 1 + x 2 = 1 , x 2 = 4 1 1 + 4 4 + 5 x 3 = 8 , x 3 = 23 / 5

SOLUTION: x = 1 4 23 / 5 T .

Algorithm 9.5

Solving a Lower Triangular System

function FORSOLVE(L,b)

%   Find the solution to the system Lx  = b,   where L is an n  × n lower-triangular matrix.

x 1  = b 1/l 11

for i   =   2:n do

  sum   =   0.0

for j   =   1:i-1 do

sum  = sum  + lijxj

end for

x  (i)   =   (b  (i)   sum)   /lii

end for

return x

end function

NLALIB: The function forsolve implements Algorithm 9.5.

Example 9.8

Solve the systems

1 1 3 0 2 9 0 0 1 x 1 x 2 x 3 = 1 9 2

and

1 0 0 1 2 0 3 4 5 y 1 y 2 y 3 = 1 9 2 .

>>   U   =   [1   -1 3;0 2 9;0 0 1];

>>   L   =   [1 0 0;-1 2 0;3 4 5];

>>   b   =   [1 9   -2]';

>>   x   =   backsolve(U,b)

x   =

20.5000

13.5000

-2.0000

>>   U\b

ans  =

20.5000

13.5000

-2.0000

>>   y   =   forsolve(L,b)

y   =

1

5

-5

>>   L\b

ans  =

1

5

-5

9.3.1 Efficiency Analysis

Algorithm 9.4 executes 1 division and then begins an outer loop having n    1 iterations. The inner loop executes n    (i  +   1)   +   1   = n  i times, and each loop iteration performs 1 addition and 1 multiplication, for a total of 2 (n  i) flops. After the inner loop finishes, 1 subtraction and 1 division execute. The total number of flops required is

1 + i = 1 n 1 2 n i + 2 = 1 + 2 n 1 + 2 i = 1 n 1 n i = 1 + 2 n 1 + 2 n 1 + n 2 + + 1 = 1 + 2 n 1 + 2 n n 1 2 = n 2 + n 1

Thus, back substitution is an O (n 2)(quadratic) algorithm. It is left as an exercise to show that Algorithm 9.5 has exactly the same flop count.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780123944351000090

STABILITY, INERTIA, AND ROBUST STABILITY

BISWA NATH DATTA , in Numerical Methods for Linear Control Systems, 2004

Computing the Inertia of a Symmetric Matrix

If A is symmetric, then Sylvester's law of inertia provides an inexpensive and numerically effective method for computing its inertia.

A symmetric matrix A admits a triangular factorization:

where U is a product of elementary unit upper triangular and permutation matrices, and D is a symmetric block diagonal with blocks of order 1 or 2. This is known as diagonal pivoting factorization. Thus, by Sylvester's law of inertia In(A) = In(D)). Once this diagonal pivoting factorization is obtained, the inertia of the symmetric matrix A can be obtained from the entries of D as follows:

Let D have p blocks of order 1 and q blocks of order 2, with p + 2q = n. Assume that none of the 2 × 2 blocks of D is singular. Suppose that out of p blocks of order 1, p′ of them are positive, p″ of them are negative, and p″ of them are zero (i.e., p′ + p″ + p″ = p). Then,

π ( A ) = p + q , v ( A ) = p + q , δ ( A ) = p .

The diagonal pivoting factorization can be achieved in a numerically stable way. It requires only n 3/3 flops. For details of the diagonal pivoting factorization, see Bunch (1971), Bunch and Parlett (1971), and Bunch and Kaufman (1977).

LAPACK implementation: The diagonal pivoting method has been implemented in the LAPACK routine SSYTRF.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780122035906500112

The Inverse

Richard Bronson , Gabriel B. Costa , in Matrix Methods (Third Edition), 2009

3.2 Calculating Inverses

In Section 2.3, we developed a method for transforming any matrix into row-reduced form using elementary row operations. If we now restrict our attention to square matrices, we may say that the resulting row-reduced matrices are upper triangular matrices having either a unity or zero element in each entry on the main diagonal. This provides a simple test for determining which matrices have inverses.

Theorem 1 A square matrix has an inverse if and only if reduction to row-reduced form by elementary row operations results in a matrix having all unity elements on the main diagonal.

We shall prove this theorem in the Final Comments to this chapter as

Theorem 2 An n × n matrix has an inverse if and only if it has rank n.

Theorem 1 not only provides a test for determining when a matrix is invertible, but it also suggests a technique for obtaining the inverse when it exists. Once a matrix has been transformed to a row-reduced matrix with unity elements on the main diagonal, it is a simple matter to reduce it still further to the identity matrix. This is done by applying elementary row operation (E3)–adding to one row of a matrix a scalar times another row of the same matrix–to each column of the matrix, beginning with the last column and moving sequentially toward the first column, placing zeros in all positions above the diagonal elements.

Example 1 Use elementary row operations to transform the upper triangular matrix

A = [ 1 2 1 0 1 3 0 0 1 ]

to the identity matrix.

Solution

[ 1 2 1 0 1 3 0 0 1 ] [ 1 2 1 0 1 0 0 0 1 ] { by addding to the second row (−3) times the third row [ 1 2 0 0 1 0 0 0 1 ] { by addding to the first row (−1) times the third row [ 1 0 0 0 1 0 0 0 1 ] { by addding to the first row (−2) times the second row

To summarize, we now know that a square matrix A has an inverse if and only if it can be transformed into the identity matrix by elementary row operations. Moreover, it follows from the previous section that each elementary row operation is represented by an elementary matrix E that generates the row operation under the multiplication EA. Therefore, A has an inverse if and only if there exist a sequence of elementary matrices. E 1, E 2,…, E k such that

But, if we denote the product of these elementary matrices as B, we then have BA = I, which implies that B = A −1. That is, the inverse of a square matrix A of full rank is the product of those elementary matrices that reduce A to the identity matrix! Thus, to calculate the inverse of A, we need only keep a record of the elementary row operations, or equivalently the elementary matrices, that were used to reduce A to I. This is accomplished by simultaneously applying the same elementary row operations to both A and an identity matrix of the same order, because if

then

(E k E k 1 E 3 E 2 E 1 )I=E k E k 1 E 3 E 2 E 1 =A 1 .

We have, therefore, the following procedure for calculating inverses when they exist. Let A be the n × n matrix we wish to invert. Place next to it another n × n matrix B which is initially the identity. Using elementary row operations on A, transform it into the identity. Each time an operation is performed on A, repeat the exact same operation on B. After A is transformed into the identity, the matrix obtained from transforming B will be A −1.

If A cannot be transformed into an indentity matrix, which is equivalent to saying that its row-reduced from contains at least one zero row, then A does not have an inverse.

Example 2 Invert

A = [ 1 2 3 4 ]

Solution

[ 1 2 1 0 3 4 0 1 ] [ 1 2 1 0 0 2 3 1 ] { by adding to the second row(–3) times the first row [ 1 2 1 0 0 1 3 2 1 2 ] { by multiplying the second row by ( 1 2 )

A has been transformed into row-reduced form with a main diagonal of only unity elements; it has an inverse. Continuing with transformation process, we get

[ 1 0 2 1 0 1 3 2 1 2 ] { by adding to the first row (–2) times the second row

Thus,

A 1 = [ 2 1 3 2 1 2 ] .

Example 3 Find the inverse of

A = [ 5 8 1 0 2 1 4 3 1 ] .

Solution

[ 5 8 1 1 0 0 0 2 1 0 1 0 4 3 1 0 0 1 ] [ 1 1.6 0.2 0.2 0 0 0 2 1 0 1 0 4 3 1 0 0 1 ] { by multiplying the first row by ( 0.2 ) [ 1 1.6 0.2 0.2 0 0 0 2 1 0 1 0 0 3.4 1.8 0.8 0 1 ] { by adding to the third row ( 4 ) times the first row [ 1 1.6 0.2 0.2 0 0 0 1 0.5 0 0.5 0 0 3.4 1.8 0.8 0 1 ] { by multiplying the second row by ( 0.5 ) [ 1 1.6 0.2 0.2 0 0 0 1 0.5 0 0.5 0 0 0 0.1 0.8 1.7 1 ] { by adding to the third row (3 .4) times the second row [ 1 1.6 0.2 0.2 0 0 0 1 0.5 0 0.5 0 0 0 1 8 17 10 ] { by multiplying the third row by ( 0.1 )

A has been transformed into row-reduced form with a main diagonal of only unity elements; it has an inverse. Continuing with the transformation process, we get

[ 1 1.6 0.2 0.2 0 0 0 1 0 4 9 5 0 0 1 8 17 10 ] { by adding to the second row (−0 .5) times the third row [ 1 1.6 0 1.4 3.4 2 0 1 0 4 9 5 0 0 1 8 17 10 ] { by adding to the first row (−0 .2) times the third row [ 1 0 0 5 11 6 0 1 0 4 9 5 0 0 1 8 17 10 ] . { by adding to the first row (−1 .6) times the second row

Thus,

A 1 = [ 5 11 6 4 9 5 8 17 10 ] .

Example 4 Find the inverse of

A = [ 0 1 1 1 1 1 1 1 3 ] .

Solution

[ 0 1 1 1 0 0 1 1 1 0 1 0 1 1 3 0 0 1 ] [ 1 1 1 0 1 0 0 1 1 1 0 0 1 1 3 0 0 1 ] { by interchanging the first and second rows [ 1 1 1 0 1 0 0 1 1 1 0 0 0 0 2 0 1 1 ] { by adding to the the third row (−1) times the first row [ 1 1 1 0 1 0 0 1 1 1 0 0 0 0 1 0 1 2 1 2 ] { by multiplying the third row by ( 1 2 ) [ 1 1 1 0 1 0 0 1 0 1 1 2 1 2 0 0 1 0 1 2 1 2 ] { by adding to the second row(−1) times the third row [ 1 1 0 0 3 2 1 2 0 1 0 1 1 2 1 2 0 0 1 0 1 2 1 2 ] { by adding to the first row(−1) times the third row [ 1 0 0 1 1 0 0 1 0 1 1 2 1 2 0 0 1 0 1 2 1 2 ] { by adding to the first row(−1) times the second row

Thus,

A 1 = [ 1 1 0 1 1 2 1 2 0 1 2 1 2 ] .

Example 5 Invert

A = [ 1 2 2 4 ] .

Solution

[ 1 2 1 0 2 4 0 1 ] [ 1 2 1 0 0 0 2 1 ] . { by adding to the second row (−2) tiems the first row

A has been transformed into row-reduced form. Since the main diagonal contains a zero element, here in the 2−2 position, the matrix A does not have an inverse. It is singular.

Problems 3.2

In Problems 1−20, find the inverses of the given matrices, if they exist.

1.

[ 1 1 3 4 ] ,

2.

[ 2 1 1 2 ] ,

3.

[ 4 4 4 4 ] ,

4.

[ 2 1 3 4 ] ,

5.

[ 8 3 5 2 ] ,

6.

[ 1 1 2 1 2 1 3 ] ,

7.

[ 1 1 0 1 0 1 0 1 1 ] ,

8.

[ 0 0 1 1 0 0 0 1 0 ] ,

9.

[ 2 0 1 0 1 2 3 1 1 ] ,

10.

[ 1 2 3 4 5 6 7 8 9 ] ,

11.

[ 2 0 0 5 1 0 4 1 1 ] ,

12.

[ 2 1 5 0 3 1 0 0 2 ] ,

13.

[ 3 2 1 4 0 1 3 9 2 ] ,

14.

[ 1 2 1 2 0 1 1 1 3 ] ,

15.

[ 1 2 1 3 2 4 2 3 1 ] ,

16.

[ 2 4 3 3 4 4 5 0 1 ] ,

17.

[ 5 0 1 2 1 2 2 3 1 ] ,

18.

[ 3 1 1 1 3 1 2 3 1 ] ,

19.

[ 1 1 1 2 0 1 1 1 0 0 2 3 0 0 0 2 ] ,

20.

[ 1 0 0 0 2 1 0 0 4 6 2 0 3 2 4 1 ] .

21.

Use the results of Problems 11 and 20 to deduce a theorem involving inverses of lower triangular matrices.

22.

Use the results of Problems 12 and 19 to deduce a theorem involving the inverses of upper triangular matrices.

23.

Matrix inversion can be used to encode and decode sensitive messages for transmission. Initially, each letter in the alphabet is assigned a unique positive integer, with the simplest correspondence being

A B C D E F G H I J K L M N O P Q R S T U V W X Y Z 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

Zeros are used to separate words. Thus, the message

SHE IS A SEER

is encoded

19 8 5 0 9 19 0 1 0 19 5 5 18 0

This scheme is too easy to decipher, however, so a scrambling effect is added prior to transmission. One scheme is to package the coded string as a set of 2-tuples, multiply each 2-tuple by a 2 × 2 invertible matrix, and then transmit the new string. For example, using the matrix

A = [ 1 2 2 3 ] ,

the coded message above would be scrambled into

[ 1 2 2 3 ] [ 19 8 ] = [ 35 62 ] , [ 1 2 2 3 ] [ 5 0 ] = [ 5 10 ] , [ 1 2 2 3 ] [ 9 19 ] = [ 47 75 ] , etc

and the scrambled message becomes

35 62 5 10 47 75

Note an immediate benefit from the scrambling: the letter S, which was originally always coded as 19 in each of its three occurrences, is now coded as a 35 the first time and as 75 the second time. Continue with the scrambling, and determine the final code for transmitting the above message.

24.

Scramble the message SHE IS A SEER using, matrix

A = [ 2 3 4 5 ] .

25.

Scramble the message AARON IS A NAME using the matrix and steps described in Problem 23.

26.

Transmitted messages are unscrambled by again packaging the received message into 2-tuples and multiplying each vector by the inverse of A. To decode the scrambled message

18 31 44 72

using the encoding scheme described in Problem 23, we first calculate

A 1 = [ 3 2 2 1 ]

and then

[ 3 2 2 1 ] [ 18 31 ] = [ 8 5 ] , [ 3 2 2 1 ] [ 44 72 ] = [ 12 16 ] .

The unscrambled message is

8 5 12 16

which, according to the letter-integer correspondence given in Problem 23, translates to HELP. Using the same procedure, decode the scrambled message

26 43 40 60 18 31 28 51 .

27.

Use the decoding procedure described in Problem 26, but with the matrix A given in Problem 24, to decipher the transmitted message

16 120 39 131 27 45 38 76 51 129 28 56 .

28.

Scramble the message SHE IS A SEER by packaging the coded letters into 3-tuples and then multiplying by the 3 × 3 invertible matrix

A = [ 1 0 1 0 1 1 1 1 0 ] .

Add as many zeros as necessary to the end of the message to generate complete 3-tuples.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780080922256500097

SOME FUNDAMENTAL TOOLS AND CONCEPTS FROM NUMERICAL LINEAR ALGEBRA

BISWA NATH DATTA , in Numerical Methods for Linear Control Systems, 2004

LU Factorization from Gaussian Elimination with Partial Pivoting

Since the interchange of two rows of a matrix is equivalent to premultiplying the matrix by a permutation matrix, the matrix A (k) is related to A (k − 1) by the following relation:

A ( k ) = M k P k A ( k 1 ) , k = 1 , 2 , , n 1 ,

where P k is the permutation matrix obtained by interchanging the rows k and r k of the identity matrix, and M k is an elementary lower triangular matrix resulting from the elimination process. So,

U = A ( n 1 ) = M n 1 P n 1 A ( n 2 ) = M n 1 P n 1 M n 2 P n 2 A ( n 3 ) = = M n 1 P n 1 M n 2 p n 2 M 2 P 2 M 1 P 1 A .

Setting M = M n-1 P n-1 M n-2 P n-2M 2 P 2 M 1 P 1, we have the following factorization of A:

The above factorization can be written in the form: PA = LU, where P = P n-1 P n-2P 2 P 1, U = A (n-1), and the matrix L is a unit lower triangular matrix formed out of the multipliers. For details, see Golub and Van Loan (1996, pp. 99).

For n = 4, the reduction of A to the upper triangular matrix U can be schematically described as follows:

1.

A P 1 P 1 A M 1 M 1 P 1 A = ( × × × × 0 × × × 0 × × × 0 × × × ) = A ( 1 ) .

2.

A ( 1 ) P 2 P 2 A ( 1 ) M 1 M 2 P 2 A ( 1 ) = M 2 P 2 M 1 P 1 A = ( × × × × 0 × × × 0 0 × × 0 0 × × ) = A ( 2 ) .

3.

A ( 2 ) P 3 P 3 A ( 2 ) M 3 M 3 P 3 A ( 2 ) = M 3 P 3 M 2 P 2 M 1 P 1 A = ( × × × × 0 × × × 0 0 × × 0 0 0 × ) = A ( 3 ) = U .

The only difference between L here and the matrix L from Gaussian elimination without pivoting is that the multipliers in the kth column are now permuted according to the permutation matrix P ˜ k = P n 1 P n 2 P k + 1 .

Thus, to construct L, again no explicit products or matrix inversions are needed. We illustrate this below.

Consider the case n = 4, and suppose P 2 interchanges rows 2 and 3, and P 3 interchanges rows 3 and 4.

The matrix L is then given by:

L = ( 1 0 0 0 m 31 1 0 0 m 21 m 42 1 0 m 41 m 32 m 34 1 ) .

Example 3.4.1.

k = 1

1.

The pivot entry is 7: r 1 = 3.

2.

Interchange rows 3 and 1.

P 1 = ( 0 0 1 0 1 0 1 0 0 ) , P 1 A = ( 7 8 9 4 5 6 1 2 4 ) .

3.

Form the multipliers: a 21 m 21 = 4 7 , a 31 m 31 = 1 7 .

4.

A ( 1 ) = M 1 P 1 A = ( 1 0 0 4 7 1 0 1 7 0 1 ) ( 7 8 9 4 5 6 1 2 4 ) ( 7 8 9 0 3 7 6 7 0 6 7 19 7 ) . .

k = 2

1.

The pivot entry is 6 7 , r 2 = 3 .

2.

Interchange rows 2 and 3.

P 2 = ( 1 0 0 0 0 1 0 1 0 ) , P 2 A ( 1 ) = ( 7 8 9 0 6 7 19 7 0 3 7 6 7 ) .

3.

Form the multiplier: m 32 = 1 2

A ( 2 ) = M 2 P 2 A ( 1 ) = ( 1 0 0 0 1 0 0 1 2 1 ) ( 7 8 9 0 6 7 19 7 0 3 7 6 7 ) = ( 7 8 9 0 6 7 19 7 0 0 1 2 ) .

Form L = ( 1 0 0 m 31 1 0 m 21 m 32 1 ) = ( 1 0 0 1 7 1 0 4 7 1 2 1 ) .

Verify. P A = ( 7 8 9 1 2 4 4 5 6 ) = L U . .

Flop-count. Gaussian elimination with partial pivoting requires only 2 3 n 3 flops. Furthermore, the process with partial pivoting requires at most O(n 2) comparisons for identifying the pivots.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780122035906500070

LINEAR STATE-SPACE MODELS AND SOLUTIONS OF THE STATE EQUATIONS

BISWA NATH DATTA , in Numerical Methods for Linear Control Systems, 2004

Algorithm 5.3.2.

The Schur Algorithm for e A .

Input. A ∈ ℝ nxn

Output. e A .

Step 1.

Transform A to R, an upper triangular matrix, using the QR iteration algorithm (Chapter 4):

P T A P = R .

(Note that when the eigenvalues of A are all real, the RSF is upper triangular.)

Step 2.

Compute e R = G = (g ij ):

For i = 1,…,n do

g ii = e r ii

End

For k = 1, 2,…,n − 1 do

For i = 1, 2,…,nk do

Set j = i + k

g i j = 1 ( r i i r j j ) [ r i j ( g i i g j j ) + p = i + 1 j 1 ( g i p r p j r i p g p j ) ] .

End

End

Step 3.

Compute e A = Pe R P T.

Flop-count: Computation of e R in Step 2 requires about (2n 3/3) flops.

MATCONTROL note: Algorithm 5.3.2 has been implemented in MATCONTROL function expmschr.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780122035906500094

Linear Equations and Eigensystems

George Lindfield , John Penny , in Numerical Methods (Fourth Edition), 2019

2.9 QR Decomposition

We have seen how a square matrix can be decomposed or factorized into the product of a lower and an upper triangular matrix by the use of elementary row operations. An alternative decomposition of A is into an upper triangular matrix and an orthogonal matrix if A is real, or into an upper triangular matrix and a unitary matrix if A is complex. This is called QR decomposition. Thus

A = Q R

where R is the upper triangular matrix and Q is the orthogonal, or the unitary matrix. If Q is orthogonal, Q 1 = Q T , and if Q is unitary, Q 1 = Q H . The preceding properties are very useful.

There are several procedures which provide QR decomposition; here we present Householder's method. To decompose a real matrix, Householder's method begins by defining a matrix P thus:

(2.21) P = I 2 ww T

where w is a column vector and P is symmetrical matrix. Provided w T w = 1 , P is also orthogonal. The orthogonality can easily be verified by expanding the product P T P = PP as follows:

PP = ( I 2 ww T ) ( I 2 ww T ) = I 4 ww T + 4 ww T ( ww T ) = I 4 ww T + 4 w ( w T w ) w T = I

since w T w = 1 .

To decompose A into QR, we begin by forming the vector w 1 from the coefficients of the first column of A as follows:

w 1 T = μ 1 [ ( a 11 s 1 ) a 21 a 31 a n 1 ]

where

μ 1 = 1 2 s 1 ( s 1 a 11 )  and s 1 = ± ( j = 1 n a j 1 2 ) 1 / 2

By substituting for μ 1 and s 1 in w 1 it can be verified that the necessary orthogonality condition, w 1 T w 1 = 1 , is satisfied. Substituting w 1 into (2.21) we generate an orthogonal matrix P ( 1 ) .

The matrix A ( 1 ) is now created from the product P ( 1 ) A . It can easily be verified that all elements in the first column of A ( 1 ) are zero except for the element on the leading diagonal which is equal to s 1 . Thus,

A ( 1 ) = P ( 1 ) A = [ s 1 + + 0 + + 0 + + 0 + + ]

In the matrix A ( 1 ) , + indicates a non-zero element.

We now begin the second stage of the orthogonalization process by forming w 2 from the coefficients of the second column of A ( 1 ) thus:

w 2 T = μ 2 [ 0 ( a 22 ( 1 ) s 2 ) a 32 ( 1 ) a 42 ( 1 ) a n 2 ( 1 ) ]

where a i j are the coefficients of A and

μ 2 = 1 2 s 2 ( s 2 a 22 ( 1 ) )  and s 2 = ± ( j = 2 n ( a j 2 ( 1 ) ) 2 ) 1 / 2

Then the orthogonal matrix P ( 2 ) is generated from

P ( 2 ) = I 2 w 2 w 2 T

The matrix A ( 2 ) is then created from the product P ( 2 ) A ( 1 ) as follows:

A ( 2 ) = P ( 2 ) A ( 1 ) = P ( 2 ) P ( 1 ) A = [ s 1 + + 0 s 2 + 0 0 + 0 0 + ]

Note that A ( 2 ) has zero elements in its first two columns except for the elements on and above the leading diagonal. We can continue this process n 1 times until we obtain an upper triangular matrix R. Thus,

(2.22) R = P ( n 1 ) P ( 2 ) P ( 1 ) A

Note that since P ( i ) is orthogonal, the product P ( n 1 ) ... P ( 2 ) P ( 1 ) is also orthogonal.

We wish to determine the orthogonal matrix Q such that A = QR . Thus R = Q 1 A or R = Q T A . Hence, from (2.22),

Q T = P ( n 1 ) P ( 2 ) P ( 1 )

Apart from the signs associated with the columns of Q and the rows of R, the decomposition is unique. These signs are dependent on whether the positive or negative square root is taken in determining s 1 , s 2 , etc. Complete decomposition of the matrix requires 2 n 3 / 3 multiplications and n square roots. To illustrate this procedure consider the decomposition of the matrix

A = [ 4 2 7 6 2 3 3 4 4 ]

Thus,

s 1 = 4 2 + 6 2 + 3 2 = 7.8102 μ 1 = 1 2 × 7.8102 × ( 7.8102 4 ) = 0.1296 w 1 T = 0.1296 [ ( 4 7.8102 ) 6 3 ] = [ 0.4939 0.7777 0.3889 ]

Using (2.21) we generate P ( 1 ) and hence A ( 1 ) thus:

P ( 1 ) = [ 0.5121 0.7682 0.3841 0.7682 0.2097 0.6049 0.3841 0.6049 0.6976 ]

A ( 1 ) = P ( 1 ) A = [ 7.8102 2.0486 2.8168 0 4.3753 3.5873 0 0.8123 7.2936 ]

Note that we have reduced the elements of the first column of A ( 1 ) below the leading diagonal to zero. We continue with the second stage thus:

s 2 = 4.3753 2 + 0.8123 2 = 4.4501 μ 2 = 1 2 × 4.4501 × ( 4.4501 + 4.3753 ) = 0.1128 w 2 T = 0.1128 [ 0 ( 4.3753 4.4501 ) 0.8123 ] = [ 0 0.9958 0.0917 ]

P ( 2 ) = [ 1 0 0 0 0.9832 0.1825 0 0.1825 0.9832 ]

R = A ( 2 ) = P ( 2 ) A ( 1 ) = [ 7.8102 2.0486 2.8168 0 4.4501 2.1956 0 0 7.8259 ]

Note that we have now reduced the first two columns of A ( 2 ) below the leading diagonal to zero. This completes the process to determine the upper triangular matrix R. Finally, we determine the orthogonal matrix Q as follows:

Q = ( P ( 2 ) P ( 1 ) ) T = [ 0.5121 0.6852 0.5179 0.7682 0.0958 0.6330 0.3841 0.7220 0.5754 ]

It is not necessary for the reader to carry out the preceding calculations since Matlab provides the function qr to carry out this decomposition. For example,

>> A = [4 -2 7;6 2 -3;3 4 4]

A =

4   -2   7

6   2   -3

3   4   4

>> [Q R] = qr(A)

Q =

-0.5121   0.6852   0.5179

-0.7682   -0.0958   -0.6330

-0.3841   -0.7220   0.5754

R =

-7.8102   -2.0486   -2.8168

0   -4.4501   2.1956

0   0   7.8259

Note that the matrices Q and R in the Matlab output are the negative of the hand calculations of Q and R above. This is not significant since their product is equal to A, and in the multiplication, the signs cancel.

One advantage of QR decomposition is that it can be applied to non-square matrices, decomposing an m × n matrix into an m × m orthogonal matrix and an m × n upper triangular matrix. Note that if m > n , the decomposition is not unique.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128122563000117

Determinants and Eigenvalues

Stephen Andrilli , David Hecker , in Elementary Linear Algebra (Fifth Edition), 2016

Calculating the Determinant by Row Reduction

We will now illustrate how to use row operations to calculate the determinant of a given matrix A by finding an upper triangular matrix B that is row equivalent to A.

Example 4

Let

A = 0 14 8 1 3 2 2 0 6 .

We row reduce A to upper triangular form, as follows, keeping track of the effect on the determinant at each step:

Because the last matrix B is in upper triangular form, we stop. (Notice that we do not target the entries above the main diagonal, as in reduced row echelon form.) From Theorem 3.2, B = ( 1 ) ( 1 ) 46 7 = 46 7 . Since |B| = + 1 14 | A | , we see that | A | = 14 | B | = 14 ( 46 7 ) = 92 .

A more convenient method of calculating A is to create a variable P (for "product") with initial value 1, and update P appropriately as each row operation is performed. That is, we replace the current value of P by

P × c for Type (I) row operations P × ( 1 ) for Type (III) row operations .

Of course, Type (II) row operations do not affect the determinant. Then, using the final value of P, we can solve for |A| using A = ( 1 / P ) B , where B is the upper triangular result of the row reduction process. This method is illustrated in the next example.

Example 5

Let us redo the calculation for A in Example 4. We create a variable P and initialize P to 1. Listed below are the row operations used in that example to convert A into upper triangular form B, with B = 46 7 . After each operation, we update the value of P accordingly.

Row Operation Effect P
(III): 1 2 Multiply P by − 1 − 1
(II): 3 2 1 + 3 No change − 1
(I): 2 1 14 2 Multiply P by 1 14 1 14
(II): 3 6 2 + 3 No change 1 14
Then A equals the reciprocal of the final value of P times B ; that is, A = ( 1 / P ) | B | = 14 × 46 7 = 92 .

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128008539000037

Additional Applications

Stephen Andrilli , David Hecker , in Elementary Linear Algebra (Fifth Edition), 2016

Exercises for Section 8.10

1.

In each part of this exercise, a quadratic form Q: R n R is given. Find an upper triangular matrix C and a symmetric matrix A such that, for every x R n , Q(x) =x T Cx =x T Ax.

⋆(a)

Q([x,y]) = 8x 2 − 9y 2 + 12xy

(b)

Q([x,y]) = 7x 2 + 11y 2 − 17xy

⋆(c)

Q ( [ x 1 , x 2 , x 3 ] ) = 5 x 1 2 2 x 2 2 + 4 x 1 x 2 3 x 1 x 3 + 5 x 2 x 3

2.

In each part of this exercise, use the Quadratic Form Method to diagonalize the given quadratic form Q: R n R . Your answers should include the matrices A, P, and D defined in that method, as well as the orthonormal basis B. Finally, calculate Q(x) for the given vector x in the following two different ways: first, using the given formula for Q, and second, calculating Q = [ x ] B T D [ x ] B where [x] B =P −1 x and D =P −1 AP.

⋆(a)

Q([x,y]) = 43x 2 + 57y 2 − 48xy;x = [1,−8]

(b)

Q ( [ x 1 , x 2 , x 3 ] ) = 5 x 1 2 + 37 x 2 2 + 49 x 3 2 + 32 x 1 x 2 + 80 x 1 x 3 + 32 x 2 x 3 ; x = [7,−2,1]

⋆(c)

Q ( [ x 1 , x 2 , x 3 ] ) = 18 x 1 2 68 x 2 2 + x 3 2 + 96 x 1 x 2 60 x 1 x 3 + 36 x 2 x 3 ; x = [4,−3,6]

(d)

Q ( [ x 1 , x 2 , x 3 , x 4 ] ) = x 1 2 + 5 x 2 2 + 864 x 3 2 + 864 x 4 2 24 x 1 x 3 + 24 x 1 x 4 + 120 x 2 x 3 + 120 x 2 x 4 + 1152 x 3 x 4 ; x = [5,9,−3,−2]

3.

Let Q: R n R be a quadratic form, and let A and B be symmetric matrices such that Q(x) =x T Ax =x T Bx. Prove that A = B (the uniqueness assertion from Theorem 8.14). (Hint: Use x =e i to show that a ii = b ii . Then use x =e i +e j to prove that a ij = b ij when ij.)

⋆4.

Let Q: R n R be a quadratic form. Is the upper triangular representation for Q necessarily unique? That is, if C 1 and C 2 are upper triangular n × n matrices with Q(x) =x T C 1 x =x T C 2 x, for all x R n , must C 1 =C 2? Prove your answer.

5.

A quadratic form Q(x) on R n is positive definite if and only if both of the following conditions hold:

(i)

Q(x) ≥ 0, for all x R n .

(ii)

Q(x) = 0 if and only if x = 0.

A quadratic form having only property (i) is said to be positive semidefinite.

Let Q be a quadratic form on R n , and let A be the symmetric matrix such that Q(x) =x T Ax.

(a)

Prove that Q is positive definite if and only if every eigenvalue of A is positive.

(b)

Prove that Q is positive semidefinite if and only if every eigenvalue of A is nonnegative.

⋆6.

True or False:

(a)

If Q(x) =x T Cx is a quadratic form, and A = 1 2 ( C + C T ) , then Q(x) =x T Ax.

(b)

Q(x,y) = xy is not a quadratic form because it has no x 2 or y 2 terms.

(c)

If x T Ax =x T Bx for every x R n , then A = B.

(d)

Every quadratic form can be diagonalized.

(e)

If A is a symmetric matrix and Q(x) =x T Ax is a quadratic form that diagonalizes to Q ( x ) = [ x ] B T D [ x ] B , then the main diagonal entries of D are the eigenvalues of A.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128008539000086