Matrix Operations
A major part of numerical methods deals with the mechanics of matrix manipulations. Code is given but it is impossible to appreciate without a reasonable understanding of the field and basic concepts.
A good text is: Matrix Computations by Gene Golub and Charles van Loan, Johns Hopkins University Press. I will supplement the material in the chapter with material from this reference. Another is Numerical Recipes in C++, 2nd Edition, Press, and others, Cambridge Press.
The basic problem is to solve Ax = b where A is a matrix of coefficients, x is a column vector of unknowns, and b is a column vector of RHS (right hand sides) of equations.
To solve the matrix equation A x = b where
the first row of A isA11 A12 ....A1n and the last row is An1 An2 ... Ann
x is a column vector whose elements are x1 x2 ...xn
b is a column vector whose elements are b1 b2 ...bn
The brute force way to solve the equation is to find and use the matrix inverse.
A-1Ax = A-1b, which results in x = b, the solution
A-1 can be found from the relation AA-1 = I where I is the identity matrix (1s along the diagonal,
0s elsewhere)
The text mentions that matrices follow some basic laws of arithmetic and not others. These laws, for numbers, are provided below. Which apply to matrices?
Arithmetic Laws
Commutative: a + b = b + a, a - b = b - a, ab = ba. Why not true for matrices?
Associative: (a + b) + c = a + (b + c). True for matrices?
Distributive: a * (b + C) = ab + ac. True for matrices?
Matrix Operations
Addition of matrices: add the corresponding elements
Subtraction of matrices: subtract the corresponding elements
Multiplication of matrices: rotate and multiply - note dimension requirements. The number of columns of A must equal the number of rows of B for the product AB to be defined.
Division of matrices: not defined but the inverse matrix is defined so that A*A-1= I
where I is the identity matrix. This holds for square matrices.
Identity matrix:
the main diagonal elements are 1 and the rest are 0.
Upper triangular matrix:
all elements below the main diagonal are zero.
Lower triangular matrix:
all elements above the main diagonal are zero.
Banded matrix:
all zero elements except for a band centered on the main diagonal.
Sparse matrix:
most of the elements are zero. How could this occur?
Inverse of a matrix A:
Matrix that when multiplied by A results in the Identity matrix.It is written as A-1
Transpose of a matrix A
It is a matrix with aij becoming aji.
Transpose of A
Designated AT. Some transpose relationships:
(AT)T = A
(A + B)T = AT + BT
(AB)T = BTAT
If A-1 exists, then (A-1)T = (AT)-1
Determinant of a matrix.
The determinant of a 2x2 matrix A is a11xa22-a21xa12
If the determinant is zero then no solution exists or many solutions exist. In general it is found as follows.
If A is an n by n matrix then the minor Mij is the determinant of the (n-1)x(n-1) submatrix of A obtained by deleting the ith row and jth foclumn of hte matrix A.
Then the determinant of A is given by
det A = |
|
(-1)i+j aijMij for any i = 1,2,...n |
The basic tasks (problems) in computational linear algebra usually fall into one of the following categories. In each case, the matrix equation Ax = b is being addressed. The complete treatment of any of these is a lengthy task. In many (maybe most) case the approach is to use a series of good routines. Many are available commercially. One of the best known free versions is LINPAK, developed at Argonne National laboratories.
Solve the equation for an unknown vector x where b is a known vector and the matrix
of coefficients A is known.
Solution of more than one matrix equation for a set of vectors x, each corresponding to
a vector b. In this case the A's are held constant while the b's vary.
Singular value decomposition of matrix A. This approach is used in the case where
M (number of equations) is less than N (number of unknowns)
If there are as many equations as unknowns (M = N) then three is a good chance of solving for a unique solution x. Analytically, there can fail to be a unique solution, however, if two or more of the M equations is a linear combination of the others - row degeneracy - or if all equations contain certain variables only in the same linear combination - column degeneracy.
Equations that are degenerate are called singular. Numerically, two things (at least) can go wrong.
Even though not exact linear combinations of each other, some of the equations may be
so close to being linearly dependent that round off errors render them so.
Accumulated round off errors can cause problems - especially for large N.
Most of the software packages and approaches (like LINPAK mentioned above) are devoted to addressing one of these two issues.
Direct Methods
A direct method is one that gives the exact solution to the system, if it is assumed that all of the calculations can be performed without the effects of round-off errors. This assumption, however, is idealized. The effects of these ever present errors must always be considered. A major goal in developing all NM algorithms is to arrange the calculations so that these effects are minimized.
Iterative Methods
Iterative methods will not return the exact solution even if all the calculations could be performed using exact arithmetic. In many instances, however, they are more effective than the direct methods. This is true because the require less computational effort and round off error is reduced. This is especially true for sparce matrices.
An added difficulty, however, is that the question of convergence must be addressed. The objective, after all, is to obtain an approximation that is wihtin a certin tolerance of the exact solution. The calculation of norms is used for this purpose - not addressed here.
Gauss-Jordan Elimination
Variables are eliminated from the equations one equation involves only one variable, a second inovlves that variable and one other, etc. The solution is found by solving for hte variable in the single equation, using htis to reduce the second equation until all variables are found.
Three operations are permitted:
An equation can be multiplied by any non-zero constant and this equation used in
place of the original equation.
An equaiton can be multiplied by any constant and added to another equation, with
the resulting equation used in place of the original.
Two equations can be transposed in order
The approach is to form an "augmented" matrix formed by placing b as a new
final column of A. When a situation is reached when there is one constant in column
1, two in column 3 etc a solution can be found. The process is called Backward
Substitution.
Example
2x + y + z = 8
1/2 y + 1/2 z = 1
2y + z = 5
Eliminate y from the first equation by adding -2 times the second equation to the first
Eliminate y from the third equation by adding -4 times the second equation to the third
2x - 2z = 6
1/2 y + 1/2 z = 1
-z = 1
Eliminate z from the first equation by adding -2 times the third equation to the first
and then eliminate z from the second equation by adding 0.5 itmes the third equation to the second
2x = 4
1/2y = 3/2
-z = 1
The version used in the code is similar except that the objective is to obtain 1s along the diagonal of A and 0s elsewhere
Matrix Factorization - LU Decomposition
This involves using Gaussian elimination to factor a matrix into a product of
matrices that are easier to manipulate.
The factorization is particularly useful when it has the form A = LU, whrere L is
lower triangular and U is upper triangular.
Now Ax = b can be written as LUx = b.
We then use a two-step process and let y = Ux. We solve the system Ly = b for y.
Now we solve Y = Ux for x.
It can be shown that this greatly reduces the number of steps required.
Special Matrices
Solution techniques are facilitated for certain types of matrices.
Strictly Diagonally Dominant
A matrix A is strictly diagonally dominant when abs (aii) > the summation of from
j = 1 to n (j does not equal i) of the absolute values of aij
A matrix with the following row vectors is strictly diagonally dominant.
(7 2 0), (3 5 -1), (0, 5 -6)
A strictly diagonally dominant matrix A is non-singular. In addition, in this case
Gaussian elimination can be performed on any linear system of the form Ax = b to
obtain its unique solution without row or column interchanges, and the computations
are stable with respect to the growth of round-off errors.