Matrices-Unit7.mws

MATRICES AND MATRIX OPERATIONS: Unit 7

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

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

(7) Special types of matrices

OBJECTIVES :

To define the five special types of matrices that can be automatically generated with Maple .

To provide alternative methods of inputting the special matrices.

To introduce the concept of a scalar matrix.

To introduce the concept of the unit matrix or identity matrix .

To introduce the concept of a singular or non-invertible matrix.

To introduce the concept of matrix sparsity .

To specify some properties of the special types of matrices.

> restart : with(linalg, coldim, det, diag, eigenvals, inverse, randmatrix, rowdim, transpose) :

Normally, when defining and inputting a matrix, it is necessary to specify and enter every element of the matrix.

There are some kinds of matrices that may be generated with Maple 'automatically'. All these special types of matrices are included and discussed in this Unit.

A . The diagonal matrix

The diagonal matrix is a square matrix that has non- zero elements only on its principal ( main , leading ) diagonal, which runs from the top left of the matrix to the bottom right.

In Maple , there are two alternative methods of defining and inputting diagonal matrices.

Consider a ( 4 × 4 ) diagonal matrix [ A ].

Method 1 . Using the diag function:

> A := diag(a[11], 0, a[33], a[44]) : A = matrix(A) ;

A = matrix([[a[11], 0, 0, 0], [0, 0, 0, 0], [0, 0, ...

Method 2 . Using the diagonal indexing function at either of the two locations under the array function:

> A := array(diagonal, 1..4, 1..4, [(1,1)=a[11], (2,2)=0, (3,3)=a[33], (4,4)=a[44]]) :

> A := array(1..4, 1..4, [(1,1)=a[11], (2,2)=0, (3,3)=a[33], (4,4)=a[44]], diagonal) : A = matrix(A) ;

A = matrix([[a[11], 0, 0, 0], [0, 0, 0, 0], [0, 0, ...

N.B. Some of the elements in the principal diagonal may also be zeros .

For example, substituting arbitrary numerical values for the non- zero elements of the matrix [ A ] gives

> A := subs(a[11]=3, a[33]=2, a[44]=5, matrix(A)) : A = matrix(A) ;

A = matrix([[3, 0, 0, 0], [0, 0, 0, 0], [0, 0, 2, 0...

* * *

N.B. A diagonal matrix, in which all the diagonal elements are equal is called a scalar matrix.

For example, construct an ( n × n ) scalar matrix [ A ] whose elements are d . Let n = 4 .

This may be done using either of the preceding Methods 1 or 2. Where large diagonal matrices are involved, the method using the double for -loop construct and the diagonal indexing function could be more convenient, viz.

> n := 4 :

> A := array(diagonal, 1..n, 1..n) :

> for i to n do for j to n do if j = i then A[i,j] := d fi : od : od :

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

A = matrix([[d, 0, 0, 0], [0, d, 0, 0], [0, 0, d, 0...

> i := 'i' :

* * *

N.B. A scalar matrix, in which all the diagonal elements are unity is called the unit matrix or identity matrix .

[ For the unit matrix , refer to Unit (9). ]

* * *

N.B. A diagonal matrix is not changed by transposition .

[ For the transpose of a matrix, refer to Unit (10). ]

* * *

N.B. Any square matrix can be reduced to a diagonal matrix having non- zero elements only in the principal diagonal.

* * *

B . The symmetric matrix

The symmetric matrix is a square matrix, in which the ( i , j ) th element is equal to the ( j , i ) th element, so that the pattern of the matrix elements has a symmetry about the principal diagonal.

In Maple , there are two alternative methods of defining and inputting symmetric matrices.

Consider a ( 3 × 3 ) symmetric matrix [ A ].

Method 1 . Using the symmetric indexing function at either of the two locations under the array function:

> A := array(symmetric, 1..3, 1..3, [(1,1)=5, (2,2)=2, (3,3)=7, (1,2)=1, (1,3)=3, (2,3)=-2]) :

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

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

Method 2 . Using any of the methods found in Unit (1), e.g.:

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

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

* * *

N.B. Raising a symmetric matrix to any (positive or negative) integer power yields another symmetric matrix.

For example, compute powers 2 and -3 of the above matrix [ A ].

> `A^2` := evalm(A^2) : `A^(-3)` := evalm(A^(-3)) :

> A^2 = matrix(`A^2`) ; A ^ ` -3` = matrix(`A^(-3)`) ;

A^2 = matrix([[35, 1, 34], [1, 9, -15], [34, -15, 6...

A^` -3` = matrix([[13334/2197, -1798/169, -12989/21...

Both resultant matrices are symmetric.

[ For integer exponentiation of matrices, refer to Unit (15). ]

* * *

Another example of a symmetric matrix is the following ( 3 × 3 ) matrix [ B ]:

> B := matrix(3, 3, [1, a, a^2, a, a^2, 1, a^2, 1, a]) : B = matrix(B) ;

B = matrix([[1, a, a^2], [a, a^2, 1], [a^2, 1, a]])...

* * *

N.B. A symmetric matrix is not changed by transposition . Based on this property, some sources define the symmetric matrix as a matrix, which is equal to its transpose, viz.

[ A ] = Transp [ A ]

Pre-multiplication of either side of this formula by the inverse of [ A ] yields the identity

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

where [ U ] is the unit matrix of the same size as [ A ]. Some sources define the symmetric matrix as a matrix, which satisfies the above identity.

[ For the inverse of a matrix, refer to Unit (14). ]

Notice that the sum [ A ] + Transp [ A ] is also a symmetric matrix.

For example, consider a ( 3 × 3 ) matrix [ A ] given as

> A := matrix(3, 3, [a, b, c, d, e, f, g, h, i]) : A = matrix(A) ;

A = matrix([[a, b, c], [d, e, f], [g, h, i]])

The sum [ A ] + Transp [ A ] is the following ( 3 × 3 ) matrix:

> `A+transp(A)` := evalm(A + transpose(A)) : A + Transp(A) = matrix(`A+transp(A)`) ;

A+Transp(A) = matrix([[2*a, b+d, c+g], [b+d, 2*e, f...

which is symmetric – see also Unit (10) for an interesting use of this property.

* * *

N.B. A symmetric matrix can be reduced to a diagonal matrix.

* * *

C . The skew-symmetric (antisymmetric) matrix

The skew-symmetric (antisymmetric) matrix is a square matrix, in which the ( i , j ) th element is the negative of the ( j , i ) th element, so that all the principal-diagonal elements are necessarily zeros , whilst the pattern of elements has a symmetry about the principal diagonal but with a reversal of sign.

In Maple , there are two alternative methods of defining and inputting skew-symmetric matrices.

Consider a ( 3 × 3 ) skew-symmetric matrix [ A ].

Method 1 . Using the antisymmetric indexing function at either of the two locations under the array function:

> A := array(antisymmetric, 1..3, 1..3, [(1,2)=1, (1,3)=5, (2,3)=-3]) :

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

A = matrix([[0, 1, 5], [-1, 0, -3], [-5, 3, 0]])

Method 2 . Using any of the methods found in Unit (1), e.g.:

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

A = matrix([[0, 1, 5], [-1, 0, -3], [-5, 3, 0]])

* * *

N.B. Raising an antisymmetric matrix to an even integer power yields a symmetric matrix, while an odd integer power yields an antisymmetric matrix.

For example, compute powers 2 and 3 of the above matrix [ A ].

> `A^2` := evalm(A^2) : `A^3` := evalm(A^3) : A^2 = matrix(`A^2`) ; A^3 = matrix(`A^3`) ;

A^2 = matrix([[-26, 15, -3], [15, -10, -5], [-3, -5...

A^3 = matrix([[0, -35, -175], [35, 0, 105], [175, -...

The square of [ A ] is symmetric , but the cube of [ A ] is antisymmetric matrix.

* * *

Another example of a skew-symmetric matrix is the following ( 3 × 3 ) matrix [ B ]:

> B := matrix(3, 3, [0, a, b, -a, 0, c, -b, -c, 0]) : B = matrix(B) ;

B = matrix([[0, a, b], [-a, 0, c], [-b, -c, 0]])

* * *

N.B. An antisymmetric matrix changes sign in transposition but is otherwise unchanged. Based on this property, some sources define the antisymmetric matrix as a matrix, which is equal to its negative transpose, viz.

[ A ] = Transp [ A ]

Notice that the difference [ A ] Transp [ A ] is also an antisymmetric matrix.

For example, consider a ( 3 × 3 ) matrix [ A ] given as

> A := matrix(3, 3, [a, b, c, d, e, f, g, h, i]) : A = matrix(A) ;

A = matrix([[a, b, c], [d, e, f], [g, h, i]])

The difference [ A ] Transp [ A ] is the following ( 3 × 3 ) matrix:

> `A-transp(A)` := evalm(A - transpose(A)) : A - Transp(A) = matrix(`A-transp(A)`) ;

A-Transp(A) = matrix([[0, b-d, c-g], [d-b, 0, f-h],...

which is antisymmetric – see also Unit (10) for an interesting use of this property.

Notice the result of summation of [ A ] + Transp [ A ] and [ A ] Transp [ A ], viz.

> `A+transp(A) + A-transp(A)` := evalm((A + transpose(A)) + (A - transpose(A))) :

> [A + Transp(A)] * ` + ` * [A - Transp(A)] = matrix(`A+transp(A) + A-transp(A)`) ;

[A+Transp(A)]*` + `*[A-Transp(A)] = matrix([[2*a, 2...

which is 2 [ A ] – see Unit (10) for the application of this result.

* * *

N.B. The skew-symmetric matrix is a singular or non-invertible matrix, i.e. its determinant is zero .

For example, the determinant of the above antisymmetrix matrix [ A ] Transp [ A ] is

> `det(A-transp(A))` := det(`A-transp(A)`) : Det(A - Transp(A)) = `det(A-transp(A))` ;

Det(A-Transp(A)) = 0

[ Refer to Unit (11) for the determinant and to Unit (14) for the inverse of a square matrix. ]

* * *

D . The sparse matrix

The sparse matrix is a rectangular matrix, in which most of the elements are zeros .

In Maple , there are two alternative methods of defining and inputting sparse matrices.

Consider a ( 3 × 4 ) sparse matrix [ A ].

Method 1 . Using the sparse indexing function at either of the two locations under the array function:

> A := array(sparse, 1..3, 1..4, [(1,3)=2, (2,1)=-3, (2,4)=4, (3,2)=-1]) :

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

A = matrix([[0, 0, 2, 0], [-3, 0, 0, 4], [0, -1, 0,...

Method 2 . Using any of the methods found in Unit (1), e.g.:

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

A = matrix([[0, 0, 2, 0], [-3, 0, 0, 4], [0, -1, 0,...

Another example of a sparse matrix is the following ( 6 × 6 ) matrix [ Y ]:

> Y := array(sparse, 1..6, 1..6, [(1,1)=15, (1,2)=-10, (1,3)=-5, (2,1)=-10, (2,2)=10, (3,1)=-5, (3,3)=25, (3,4)=-12, (3,5)=-8, (4,3)=-12, (4,4)=12, (5,3)=-8, (5,5)=23, (5,6)=-15, (6,5)=-15, (6,6)=15]) : Y = matrix(Y) ;

Y = matrix([[15, -10, -5, 0, 0, 0], [-10, 10, 0, 0,...

* * *

N.B. The matrix [ Y ] was obtained when analysing electric power flow in a real-life six bus-bar transmission system. Multiplication of every element of [ Y ] by the imaginary unit yields a so-called bus admittance matrix of the system.

It is typical of bus admittance matrices that the vast majority of their elements have the value zero . Such matrices are said to be highly sparse . The notion of matrix sparsity is used to determine how sparse such matrices are.

Since realistic-size systems may have in excess of 100 bus-bars, it is convenient to have a method for fast determination of the sparsity of bus admittance matrices. A suitable method is proposed below, where the above matrix [ Y ] is used as an example.

(a) Find the number of all elements of the matrix:

> dim(Y) := rowdim(Y) * coldim(Y) : all_elements(Y) = dim(Y) ;

all_elements(Y) = 36

(b) Find the number of zero -valued elements of the matrix:

> n := 0 : for i to rowdim(Y) do for j to coldim(Y) do if Y[i,j] = 0 then n := n+1 fi : od : od :

> zeros(Y) := n : zero_elements(Y) = zeros(Y) ;

zero_elements(Y) = 20

(c) Compute the matrix sparsity according to the formula

> sparsity(Y) = zero_elements(Y)/all_elements(Y) ; sparsity(Y) = evalf(zeros(Y)/dim(Y), 3) ;

sparsity(Y) = zero_elements(Y)/all_elements(Y)

sparsity(Y) = .556

* * *

N.B. This particular matrix [ Y ] is singular since its determinant is zero , as evidenced by the following:

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

Det(Y) = 0

This implies that matrix [ Y ] has no inverse .

* * *

E . The unimodular matrix

The unimodular matrix is a square integer matrix that has all elements on its principal diagonal equal to unity . The determinant of a unimodular matrix is unity .

In Maple , a unimodular matrix may be defined and input using any of the methods found in Unit (1) or it can be generated by the program. To this end, the function randmatrix together with the type name unimodular are used.

For example, consider a ( 4 × 4 ) unimodular matrix [ A ] generated by the program.

> A := randmatrix(4, 4, unimodular) : A = matrix(A) ;

A = matrix([[1, -85, -55, -37], [0, 1, -35, 97], [0...

The determinant of [ A ] is

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

Det(A) = 1

* * *

N.B. The inverse of a unimodular matrix is also a unimodular matrix.

The inverse of the above matrix [ A ] is

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

Inv(A) = matrix([[1, 85, 3030, -159708], [0, 1, 35,...

The determinant of the inverse of [ A ] is

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

Det(Inv(A)) = 1

* * *

N.B. Transposition of a unimodular matrix returns a unimodular matrix whose elements are a mirror reflection of the elements of the original matrix about the principal diagonal.

The transpose of the above matrix [ A ] is

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

Transp(A) = matrix([[1, 0, 0, 0], [-85, 1, 0, 0], [...

The determinant of the transpose of [ A ] is

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

Det(Transp(A)) = 1

* * *

N.B. Raising a unimodular matrix to any (positive or negative) integer power yields another unimodular matrix.

For example, compute powers 2 and -3 of the above matrix [ A ].

> `A^2` := evalm(A^2) : `A^(-3)` := evalm(A^(-3)) :

> A^2 = matrix(`A^2`) ; A ^ ` -3` = matrix(`A^(-3)`) ;

A^2 = matrix([[1, -170, 2865, -11069], [0, 1, -70, ...

A^` -3` = matrix([[1, 255, 18015, -1553359], [0, 1,...

Both resultant matrices are unimodular.

* * *

N.B. All the eigenvalues of a unimodular matrix are unity .

The eigenvalues of the above unimodular matrix [ A ] are

> charroots(A) := eigenvals(A) : char_roots(A) = charroots(A) ;

char_roots(A) = (1, 1, 1, 1)

[ For the eigenvalues of a matrix, refer to Unit (21). ]

* * *

Proceed to Unit (8) for " Multiplication of a matrix by a number ".

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