MATRICES AND MATRIX OPERATIONS: Unit 14

Dr. Wlodzislaw Kostecki

The Papua New Guinea University of Technology (PNGUT)

Department of Electrical and Communication Engineering

Lae, Morobe Province

Papua New Guinea

Copyright ฉ 2000 by Wlodzislaw Kostecki

All rights reserved

-------------------------------------------------------------------

(14) The inverse of a matrix

OBJECTIVES :

• To provide alternative definitions of the inverse or reciprocal of a square matrix and state the necessary condition for its existence.

• To provide alternative methods of matrix inversion with Maple .

• To specify and illustrate some properties of the matrix inverse.

• To show how the multiplicative inverse matrix may be used in Maple for computing the "quotient" of two square non-singular matrices.

> restart : with(linalg, adj, coldim, det, diag, inverse, minor, multiply, rowdim, transpose) :

The inverse or reciprocal of a square matrix is the quotient of the adjoint of the matrix and the determinant of the matrix

Inv [ A ] = Adj [ A ]/ Det [ A ]

The inverse of a matrix exists if and only if its determinant is not . A matrix that has an inverse is called non-singular or invertible .

Consider a ( ื ) matrix [ A ] given as

> A := matrix(2, 2, [a[11], a[12], a[21], a[22]]) : A = matrix(A) ;

The matrix inverse to the matrix [ A ] is a ( ื ) matrix, which can be obtained using any of the following alternative methods.

Method 1 . Using the defining formula for the matrix inverse and the evalm function:

> `inv(A)` := evalm(adj(A)/det(A)) : Inv(A) = matrix(`inv(A)`) ;

Method 2 . Using the inverse function:

> `inv(A)` := inverse(A) : Inv(A) = matrix(`inv(A)`) ;

Method 3 . Using the power and the evalm function:

> `inv(A)` := evalm(A^(-1)) : Inv(A) = matrix(`inv(A)`) ;

N.B. An equivalent form to evalm(A^(-1)) is evalm(1/A) .

As a numerical example, consider a ( ื ) matrix [ B ] given as

> B := matrix(3, 3, [3, -2, 1, -2, 6, 4, 1, 4, 8]) : B = matrix(B) ;

The inverse of the matrix [ B ] is

> `inv(B)` := inverse(B) : Inv(B) = matrix(`inv(B)`) ;

Floating-point evaluation of the inverse of matrix [ B ] gives

> `inv(B)` := evalf(evalm(`inv(B)`)) : Inv(B) = matrix(`inv(B)`) ;

* * *

N.B. Some sources define a matrix inverse using the concept of the cofactor matrix – refer to Unit (13) for the definition of that matrix. According to this approach, the inverse is defined as follows.

The inverse of a square matrix [ A ] is the quotient of the transpose of the cofactor matrix associated with [ A ] and the determinant of [ A ]

Inv [ A ] = Transp(Cofactor [ A ] ) / Det [ A ]

Exemplarily, consider a ( ื ) matrix [ A ] given as

> A := matrix(3, 3, [2, 1, 3, 4, 2, -1, 2, -1, 1]) : A = matrix(A) ;

and obtain the inverse of [ A ].

Step 1 . Compute the cofactor matrix of [ A ] using the double for -loop construct:

> cofactor(A) := matrix(3, 3) :

> for i to rowdim(A) do for j to coldim(A) do cofactor(A)[i,j] := (-1)^(i+j)*det(minor(A, i, j)) : od : od :

> cofactor(A) := matrix(cofactor(A)) : Cofactor(A) = matrix(cofactor(A)) ;

Step 2 . Compute the transpose of the cofactor matrix:

> `transp(cofactor(A))` := transpose(cofactor(A)) :

> Transp(Cofactor(A)) = matrix(`transp(cofactor(A))`) ;

Step 3 . Compute the determinant of [ A ]:

> `det(A)` := det(A) : Det(A) = `det(A)` ;

Step 4 . Compute the inverse of [ A ] using the definition provided earlier:

> `inv(A)` := evalm(`transp(cofactor(A))`/`det(A)`) : Inv(A) = matrix(`inv(A)`) ;

Step 5 . Verify this result using the inverse function:

> `inv(A)` := inverse(A) : Inv(A) = matrix(`inv(A)`) ;

Both matrix inverses are equal.

* * *

N.B. If the commutative property of multiplication holds for non-singular matrices [ A ] and [ B ], then it also holds for the following pairs of matrices:

1. Inv [ A ] Inv [ B ] = Inv [ B ] Inv [ A ]

2. Inv [ A ] [ B ] = [ B ] Inv [ A ]

3. [ A ] Inv [ B ] = Inv [ B ] [ A ]

Exemplarily, consider the commuting matrices [ A ] and [ B ] that are used in Unit (4), i.e.

> A := matrix(2, 2, [6, 8, 4, 6]) : B := matrix(2, 2, [15, 20, 10, 15]) : A = matrix(A) ; B = matrix(B) ;

1(a). The product Inv [ A ] Inv [ B ] is

> `inv(A) inv(B)` := evalm(inverse(A) &* inverse(B)) : Inv(A) * Inv(B) = matrix(`inv(A) inv(B)`) ;

1(b). The product Inv [ B ] Inv [ A ] is

> `inv(B) inv(A)` := evalm(inverse(B) &* inverse(A)) : Inv(B) * `Inv(A)` = matrix(`inv(B) inv(A)`) ;

The resultant matrices of 1(a) and 1(b) are equal.

2(a). The product Inv [ A ] [ B ] is

> `inv(A) B` := evalm(inverse(A) &* B) : Inv(A) * B = matrix(`inv(A) B`) ;

2(b). The product [ B ] Inv [ A ] is

> `B inv(A)` := evalm(B &* inverse(A)) : B * `Inv(A)` = matrix(`B inv(A)`) ;

The resultant matrices of 2(a) and 2(b) are equal.

3(a). The product [ A ] Inv [ B ] is

> `A inv(B)` := evalm(A &* inverse(B)) : A * Inv(B) = matrix(`A inv(B)`) ;

3(b). The product Inv [ B ] [ A ] is

> `inv(B) A` := evalm(inverse(B) &* A) : `Inv(B)` * A = matrix(`inv(B) A`) ;

The resultant matrices of 3(a) and 3(b) are equal.

* * *

N.B. If two matrices satisfy the commutative law of multiplication and their matrix product is the unit matrix, i.e.

[ A ] [ B ] = [ B ] [ A ] = [ U ]

then matrix [ B ] is the inverse of matrix [ A ].

Exemplarily, consider ( ื ) matrices [ A ] and [ B ] given as

> A := matrix(2, 2, [3, 4, 2, 3]) : B := matrix(2, 2, [3, -4, -2, 3]) : A = matrix(A) ; B = matrix(B) ;

(a) The product [ A ] [ B ] is the ( ื ) unit matrix:

> `AB` := multiply(A, B) : `A B` = matrix(`AB`) ;

(b) The product [ B ] [ A ] is the ( ื ) unit matrix:

> `BA` := multiply(B, A) : `B A` = matrix(`BA`) ;

(c) The inverse of matrix [ A ] is the following matrix:

> `inv(A)` := inverse(A) : Inv(A) = matrix(`inv(A)`) ;

which is equal to matrix [ B ].

* * *

N.B. If inversion of a matrix [ A ] does not change the matrix, then

( [ U ] – [ A ] ) ( [ U ] + [ A ] ) = [ 0 ]

Exemplarily, consider a non-singular ( ื ) matrix [ A ] given as

> A := matrix(3, 3, [1, 0, 0, 0, 0, 1, 0, 1, 0]) : A = matrix(A) ;

(a) The inverse of [ A ] is the matrix

> `inv(A)` := inverse(A) : Inv(A) = matrix(`inv(A)`) ;

which is equal to the matrix [ A ].

(b) With the appriopriately sized unit matrix [ U ], i.e.

> U := diag(1, 1, 1) : U = matrix(U) ;

the product ( [ U ] – [ A ] ) ( [ U ] + [ A ] ) is

> `(U-A) (U+A)` := evalm((U-A) &* (U+A)) : '(U - A) * (U + A)' = matrix(`(U-A) (U+A)`) ;

which is a ( ื ) matrix.

* * *

N.B. The product of the inverse of a matrix and the matrix itself obeys the commutative law . The product matrix is the unit matrix

(Inv [ A ] ) [ A ] = [ A ] (Inv [ A ] ) = [ U ]

Exemplarily, consider a ( ื ) invertible matrix [ A ] given as

> A := matrix(3, 3, [2, 1, 0, 1, -1, 3, 0, 6, 2]) : A = matrix(A) ;

(a) The product (Inv [ A ] ) [ A ] is

> `inv(A) A` := multiply(inverse(A), A) : Inv(A)*A = matrix(`inv(A) A`) ;

(b) The product [ A ] (Inv [ A ] ) is

> `A inv(A)` := multiply(A, inverse(A)) : A * `Inv(A)` = matrix(`A inv(A)`) ;

Both product matrices are identical and equal to the ( ื ) unit matrix [ U ]

> U = matrix(U) ;

* * *

N.B. The product of the inverse and transpose of a symmetric matrix is the unit matrix

(Inv [ A ] ) (Transp [ A ] ) = [ U ]

Exemplarily, consider a ( ื ) symmetric matrix [ A ] given as

> A := matrix(3, 3, [5, 1, 3, 1, 2, -2, 3, -2, 7]) : A = matrix(A) ;

The product (Inv [ A ] ) (Transp [ A ] ) is the ( ื ) unit matrix

> `inv(A) transp(A)` := evalm(inverse(A) &* transpose(A)) :

> Inv(A) * Transp(A) = matrix(`inv(A) transp(A)`) ;

* * *

N.B. The inverse of the matrix inverse is equal to the original matrix

Inv(Inv [ A ] ) = [ A ]

Exemplarily, consider a ( ื ) matrix [ A ] given as

> A := matrix(2, 2, [a[11], a[12], a[21], a[22]]) : A = matrix(A) ;

The inverse of the inverse of the matrix, Inv(Inv [ A ] ) , is the following ( ื ) matrix:

> `inv(inv(A))` := inverse(inverse(A)) : Inv(Inv(A)) = matrix(`inv(inv(A))`) ;

This operation may be displayed in "like-in-a-book" form, viz.

> Inv(Inv(matrix(A))) = matrix(`inv(inv(A))`) ;

* * *

N.B. The transpose of the matrix inverse is equal to the inverse of the matrix transpose

Transp(Inv [ A ] ) = Inv(Transp [ A ] )

Exemplarily, consider a ( ื ) matrix [ A ] given as

> A := matrix(2, 2, [a[11], a[12], a[21], a[22]]) : A = matrix(A) ;

(a) The transpose of the matrix inverse, Transp(Inv [ A ] ) , is the following ( ื ) matrix:

> `transp(inv(A))` := transpose(inverse(A)) : Transp(Inv(A)) = matrix(`transp(inv(A))`) ;

(b) The inverse of the matrix transpose, Inv(Transp [ A ] ) , is the following ( ื ) matrix:

> `inv(transp(A))` := inverse(transpose(A)) : Inv(Transp(A)) = matrix(`inv(transp(A))`) ;

This operation may be displayed in "like-in-a-book" form, viz.

> Transp(Inv(matrix(A))) = Inv(Transp(matrix(A))) ;

* * *

N.B. The inverse of the product of two matrices is the product of the two inversed matrices but in the reverse order

Inv( [ A ] [ B ] ) = Inv [ B ] Inv [ A ]

Exemplarily, consider ( ื ) matrices [ A ] and [ B ] given as

> A := matrix(2, 2, [a[11], a[12], a[21], a[22]]) : B := matrix(2, 2, [b[11], b[12], b[21], b[22]]) :

> A = matrix(A) ; B = matrix(B) ;

(a) The inverse of the product of the two matrices, Inv( [ A ] [ B ] ) , is the following ( ื ) matrix:

> `inv(AB)` := inverse(multiply(A, B)) : Inv(A * B) = matrix(`inv(AB)`) ;

(b) The product of the two inversed matrices multiplied in the reverse order, Inv [ B ] Inv [ A ], is the following ( ื ) matrix:

> `inv(B) inv(A)` := multiply(inverse(B), inverse(A)) : Inv(B) * Inv(A) = matrix(`inv(B) inv(A)`) ;

To obtain the display of this product matrix in a compact form, the map function is used to access the matrix elements and apply the function normal with the argument expanded to each element. It is the arrow-type procedure, ( x-> ) , in which x is a local variable . In this way, the procedure is applied to each operand (matrix element) of the second argument (the matrix) of the map function.

> `inv(B) inv(A)` := map(x->normal(x, expanded), `inv(B) inv(A)`) :

> Inv(B) * Inv(A) = matrix(`inv(B) inv(A)`) ;

The product matrix Inv [ B ] Inv [ A ] displayed in the compact form is clearly equal to the matrix Inv( [ A ] [ B ] ) of (a).

* * *

N.B. The determinant of a non-singular square matrix is the reciprocal of the determinant of the inversed matrix

Det [ A ] = / Det(Inv [ A ] )

Exemplarily, consider a (2 ื 2) matrix [ A ] given as

> A := matrix(2, 2, [1, 3, 2, 4]) : A = matrix(A) ;

The inverse of [ A ] is the following (2 ื 2) matrix:

> `inv(A)` := inverse(A) : Inv(A) = matrix(`inv(A)`) ;

(a) The determinant of the matrix, Det [ A ], is

> `det(A)` := det(A) : Det(A) = `det(A)` ;

(b) The reciprocal of the determinant of the inversed matrix, / Det(Inv [ A ] ) , is

> `1/det(inv(A))` := 1/det(`inv(A)`) : 1/Det(Inv(A)) = `1/det(inv(A))` ;

N.B. This property of the determinant is used in solving systems of linear algebraic inhomogeneous equations by means of Cramer’ s rule – refer to Unit (19).

* * *

N.B. The determinant of the inverse of a non-singular square matrix is the reciprocal of the determinant of the matrix

Det(Inv [ A ] ) = / Det [ A ]

Exemplarily, consider the same (2 ื 2) matrix [ A ] as above.

(a) The determinant of the matrix inverse, Det(Inv [ A ] ) , is

> `det(inv(A))` := det(inverse(A)) : Det(Inv(A)) = `det(inv(A))` ;

(b) The reciprocal of the determinant of the matrix, / Det [ A ], is

> `1/det(A)` := 1/det(A) : 1/Det(A) = `1/det(A)` ;

* * *

N.B. If [ A ] and [ B ] are non-singular matrices of the same order, then

Det(Inv [ A ] [ B ] [ A ] ) = Det [ B ]

Exemplarily, consider ( ื ) non-singular matrices [ A ] and [ B ] given as

> A := matrix(3, 3, [2, 3, 1, 3, -1, -2, -1, 2, 1]) : B := matrix(3, 3, [12, 8, -1, 6, 4, 0, -7, -2, 1]) :

> A = matrix(A) ; B = matrix(B) ;

(a) The determinant Det(Inv [ A ] [ B ] [ A ] ) is

> `det(inv(A) B A)` := det(evalm(inverse(A) &* B &* A)) : Det(Inv(A)*B*A) = `det(inv(A) B A)` ;

(b) The determinant of [ B ] is

> `det(B)` := det(B) : Det(B) = `det(B)` ;

* * *

N.B. Although the operation of division is not defined for matrices, the multiplicative inverse matrix is used in computing the "quotient" of two square non-singular matrices. In this sense, Maple recognises "division" of two matrices and performs it according to the relationship

[ A ]/[ B ] = [ A ] Inv [ B ]

Exemplarily, consider non-singular ( ื ) matrices [ A ] and [ B ] given as

> A := matrix(3, 3, [2, 1, 0, 1, -1, 3, 0, 6, 2]) : B := matrix(3, 3, [3, -2, 1, 0, 4, -3, 2, 5, -1]) :

> A = matrix(A) ; B = matrix(B) ;

(a) The "quotient" of the matrices, [ A ]/[ B ], is returned by Maple as the following ( ื ) matrix:

> `A/B` := evalm(A/B) : A/B = matrix(`A/B`) ;

(b) The product matrix [ A ] Inv [ B ] is the following ( ื ) matrix:

> `A inv(B)` := multiply(A, inverse(B)) : A * Inv(B) = matrix(`A inv(B)`) ;

Both matrices of (a) and (b) are equal.

* * *

N.B. It follows from the relationships

[ A ]/[ B ] = [ A ] Inv [ B ]

and

Inv(Inv [ A ] ) = [ A ]

that for square non-singular matrices

Inv [ A ]/ Inv [ K ] = Inv [ A ] {Inv(Inv [ K ] )} = (Inv [ A ] ) [ K ] = [ X ]

where the matrix [ X ] satisfies the equation

[ A ] [ X ] = [ K ]

To verify this conclusion, consider the following ( ื ) matrices:

> A := matrix(3, 3, [2, 1, 0, 1, -1, 3, 0, 6, 2]) : X := matrix(3, 3, [3, -2, 1, 0, 4, -3, 2, 5, -1]) :

> K := matrix(3, 3, [6, 0, -1, 9, 9, 1, 4, 34, -20]) : A = matrix(A) ; X = matrix(X) ; K = matrix(K) ;

(a) The "quotient" of the matrices Inv [ A ]/ Inv [ K ] is the following ( ื ) matrix

> `inv(A)/inv(K)` := evalm(inverse(A)/inverse(K)) : Inv(A)/Inv(K) = matrix(`inv(A)/inv(K)`) ;

which is equal to the matrix [ X ].

(b) The product matrix [ A ] [ X ] is the following ( ื ) matrix

> `AX` := multiply(A, X) : A * X = matrix(`AX`) ;

which is equal to the matrix [ K ].

N.B. If [ A ] is a square non-singular matrix and [ K ] is a non- column matrix with the number of rows equal to that in matrix [ A ], the relationship [ A ] [ X ] = [ K ] is used for representing a system of linear algebraic inhomogeneous equations in the matrix form – refer to Unit (19).

* * *

N.B. It follows from the relationships

[ A ]/[ B ] = [ A ] Inv [ B ]

and

[ A ] (Inv [ A ] ) = [ U ]

that for a square non-singular matrix [ A ]

[ A ]/[ A ] = [ U ]

To verify this conclusion, consider the above ( ื ) matrix [ A ] and compute the result of [ A ]/[ A ].

> `A/A` := multiply(A, A^(-1)) : `A/A` = matrix(`A/A`) ;

or, using the evalm function,

> `A/A` := evalm(A &* evalm(A^(-1))) : `A/A` = matrix(`A/A`) ;

* * *

Proceed to Unit (15) for " Integer exponentiation of matrices ".

-------------------------------------------------------------------