Egwald's popular web pages are provided without cost to users. Please show your support by joining Egwald Web Services as a Facebook Fan:
Follow Elmer Wiens on Twitter:

A matrix A = [a_{i, j }] is a rectangular m by n array of numbers (or other mathematical objects) expressed as:

A =

a_{1,1} a_{1,2} . . . . . a_{1,n-1} a_{1,n}

a_{2,1} a_{2,2} . . . . . a_{2,n-1} a_{2,n}

. . . . . . . . . . .

a_{m,1} a_{m,2} . . . a_{m,n-1} a_{m,n}

such as:

A =

10 5 7 9

22 7 5 6

3 15 4 4

A matrix of dimension m by 1 is called a column vector; a matrix of dimension 1 by n is called a row vector.

The Zero Matrix

The zero matrixO = [o_{i, j}], with o_{i, j} = 0 for all i and j.

The Identity Matrix.

A square matrix is an identity matrixI = [i_{i, j}], if i_{i, j} = 1 where i = j, and i_{i, j} = 0 otherwise. Thus, an identity matrix is a diagonal matrix with 1's along its diagonal and zeros elsewhere. A 3 by 3 identity matrix appears below.

I =

1 0 0

0 1 0

0 0 1

Rules of Operations on Matrices.

Scalar Multiplication of a Matrix.

The product of a matrix A by the scalar (real or complex number) s is given by:

s * A =

s * a_{1,1} s * a_{1,2} . . . . . s * a_{1,n}

s * a_{2,1} s * a_{2,2} . . . . . s * a_{2,n}

. . . . . . . . . . .

s * a_{m,1} s * a_{m,2} . . . s * a_{m,n}

or s * A = [s * a_{i, j}], with each entry multiplied by the number s.

Addition of Matrices.

Given matrices A = [a_{i, j}] and B = [b_{i, j }] of dimension m by n, their sum A + B is another matrix C = [c_{i, j }] of dimension m by n, where:

c_{i, j } = a_{i, j} + b_{i, j} , for i = 1, . . . m; j = 1, . . . n .

For example,

5 4 10

9 3 7

-8 -3 10

+

-5 6 11

-4; 7 5

8 -3 2

=

0 10 21

5 10 12

0 -6 12

Matrix addition is commutative, since A + B = B + A.

Subtraction of Matrices.

Given matrices A and B of dimension m by n, their difference A - B is another matrix C of dimension m by n, where:

C = A - B = A + (-1) * B.

Equality of Matrices

Matrices A and B with the same dimensions are equal if A - B equals the zero matrix O.

Multiplication of two Matrices.

Given matrices A = [a_{i, j}] of dimension m by r and B = [b_{i, j }] of dimension r by n, their product A * B is another matrix C = [c_{i, j }] of dimension m by n, where:

The outer product of a column vector x = [x_{i}] (n by 1) and a row vector y^{T} = [y_{i}] (1
by n) is a matrix Z = [z_{i, j}] = x * y^{T}, where z_{i, j} = x_{i} * y_{j}.

For example,

3

-6

5

*

2

4

-3

=

6

12

-9

-12

-24

18

10

20

-15

Norm of a Vector.

The Euclidean norm ||x|| of a column vector x = [x_{i}] is the square root of x^{T} * x, ie ||x|| = (x^{T} * x)^{1/2}.

For example, if v^{T} = (3, -6, 5), ||v|| = sqrt(9 + 36 + 25) = 8.3666.

The different types of vector norms and their properties are described on the vectors web page.

Properties of Matrices

A Square Matrix.

A square matrix has the same number of columns and rows.

A Symmetric Matrix.

A square matrix A is symmetric if A = A^{T}.

Given matrices A and A^{T}defined above, then the matrices:

A * A^{T} =

255 344 169

344 594 215

169 215 266

A^{T} * A =

593 249 192 234

249 299 130 147

192 130 90 109

234 147 109 133

are symmetric matrices.

A is skew-symmetric if A = -A^{T}.

A Triangular Matrix.

A square matrix A is triangular if all entries above or below its diagonal are zero.

A square matrix A is upper triangular if all entries below its diagonal are zero.

A square matrix A is lower triangular if all entries above its diagonal are zero. The matrix below is lower triangular.

1.0 0 0

0.74 1.0 0

0.49 0.79 1.0

Diagonal Matrix.

A square matrix A is a diagonal matrix if all entries are zero, except possibly the entries on A's diagonal.

The Trace of a Matrix.

The trace of a square matrix A = [a_{i, j}] is tr(A) = a_{1, 1} + a_{2, 2} . . . . . + a_{n, n}, the sum of the diagonal elements of A.

The Rank of a Matrix.

The rank of an m by n matrix A, rank(A), equals the maximum number of A's linearly independent columns, which equals the maximum number of A's linearly independent rows. See the vectors web page for a definition of linearly dependent vectors.

The Quadratic Form of a Matrix.

Given a symmetric matrix A of dimension n, the function f(x) from R^{n} to R is called a quadratic form if it is defined by:

f(x) = x^{T} * A * x, for all vectors x in R^{n}.

Positive Definite Matrix.

A real, symmetric matrix A is called a positive definite matrix if its quadratic form is positive for all nonzero x in R^{n}. ie

f(x) = x^{T} * A * x > 0, for all nonzero vectors x in R^{n}.

Example.

The identity matrix I is positive definite because x^{T} *I * x = x_{1}^{2} + x_{2}^{2} + . . . + x_{n}^{2} > 0 for all nonzero x.

Matrix Norms.

A matrix norm ||A|| of a square matrix A of dimension n is a measure of the "distance" of A from the zero matrix, O, with properties similar to those of a vector norm.

Induced (Natural) Matrix Norms.

Many important matrix norms are derived from vector norms. If ||x|| is a vector norm for vectors of R^{n}, the matrix norm induced by ||x|| is defined by:

||A|| = max{ ||A * x|| | all x in R^{n} with ||x|| = 1 }

Properties of Matrix Norms.

For a given matrix norm || . ||, any scalar s, and square matrices A and B of dimension n:

ie ||A||_{F} equals the sum of the absolute values of the entries of A.

Special Matrices

Permutation Matrix.

A permutation matrix P is formed by rearranging the columns of an identity matrix. The permutation matrix below will switch rows 1 and 3 of a matrix A if multiplied as P * A.

P =

0 0 1

0 1 0

1 0 0

Rotation Matrices.

Rotation matrices play an important role in 3-D computer graphics. When a vector is multiplied by a rotation matrix, the vector is rotated through a given angle Ø.

Q_{Ø} =

cos(Ø)

-sin(Ø)

sin(Ø)

cos(Ø)

The matrix Q_{Ø} with Ø = 0.524 radians is:

Q_{0.524} =

0.866

-0.5

0.5

0.866

The next diagram shows the vectors v and Q_{Ø} * v with v = (5, 3) and an angle of rotation of Ø = 0.524 radians.

You can change the vector, v, and the rotation angle, Ø in radians, in the form above to see how v and Q_{Ø} * v change. For example, try v = (-5, -3), Ø = 5 radians and click "Rotate Vector".

Properties of Rotation Matrices.

Q_{Ø} * Q_{-Ø} = Q_{0}

Q_{Ø} * Q_{Ø} = Q_{2*Ø}

Q_{Ø1} * Q_{Ø2} = Q_{(Ø1+ Ø2)}

||Q_{Ø}|| = 1

Q_{Ø} * Q_{Ø}^{T} = I, the identity matrix.

||Q_{Ø} * x|| = ||x||, ie Q_{Ø} is length preserving.

Projection Matrices.

The concept of projecting one vector onto another vector was used with vector dot products, and the Gram-Schmidt Process. When a vector is multiplied by the following matrix, the vector is projected onto the line that passes through the origin making an angle of Ø radians with the x-axis.

P_{Ø} =

cos(Ø)*cos(Ø)

cos(Ø)*sin(Ø)

cos(Ø)*sin(Ø)

sin(Ø)*sin(Ø)

The matrix P_{Ø} with Ø = 0.524 radians is:

Q_{0.524} =

0.75

0.433

0.433

0.25

The next diagram shows the vectors v and P_{Ø} * v with v = (2, 4) and the projection line with an angle of Ø = 0.524 radians.

You can change the vector, v, and the projection line angle, Ø in radians, in the form above to see how v and P_{Ø} * v change. For example, try v = (-4, 4), Ø = 5 radians and click "Project Vector".

Properties of Projection Matrices.

P_{Ø} * P_{Ø} = P_{Ø}.

P_{Ø} * x = x, if x lies on the Ø line.

P_{Ø} = P_{Ø}^{T}, ie P_{Ø} is symmetric.

||P_{Ø} * x|| <= ||x||.

Reflection Matrices.

Reflection matrices act as a mirror with respect to the Ø line through the origin. When a vector is multiplied by a reflection matrix, the vector is reflected through the Ø line an equal distant to the other side of the line.

H_{Ø} =

2*cos(Ø)*cos(Ø) - 1

2*cos(Ø)*sin(Ø)

2*cos(Ø)*sin(Ø)

2*sin(Ø)*sin(Ø) - 1

The matrix H_{Ø} with Ø = 0.785 radians is:

H_{0.785} =

0

1

1

-0

The next diagram shows the vectors v and H_{Ø} * v with v = (2, 4) and the reflection line with an angle of Ø = 0.785 radians.

You can change the vector, v, and the reflection line angle, Ø in radians, in the form above to see how v and H_{Ø} * v change. For example, try v = (-4, 4), Ø = 2 radians and click "Reflect Vector".

Properties of Reflection Matrices.

H_{Ø} = 2*P_{Ø} - I.

H_{Ø} * H_{Ø} = I.

H_{Ø} = H_{Ø}^{T}.

||H_{Ø}|| = 1.

||H_{Ø} * x|| = ||x||, ie H_{Ø} is length preserving.

Orthogonal and Orthonormal Matrices.

A square matrix A is orthogonal if for each column a_{i} of A, a_{i}^{T} * a_{j} = 0 for any other column a_{j} of A. If each column a_{i} of A has a norm of 1, then A is orthonormal. If A is orthonormal, then:

A^{T} * A = I, the identity matrix.

Example of an Orthonormal Matrix.

Q =

-0.4472 -0.8944 0

0 0 1

0.8944 -0.4472 0

Permutation, rotation, and reflection matrices are orthonormal matrices.

Orthogonal / orthonormal matrices play an important role in the transformation and decomposition of matrices into special forms.

Inverse Matrix.

A square matrix A has an inverse if a matrix B exists with B * A = I. A's inverse is usually written as A^{-1}.

Properties of A^{-1}.

If A is a square matrix of dimension n, then A^{-1} obeys:

A^{-1} * A = A * A^{-1} = I

A^{-1} is unique if it exists

(A^{-1})^{-1} = A

A^{T} = A^{-1} if A is orthonormal

A^{-1} exists if and only if rank(A) = n.

(A^{T})^{-1} = (A^{-1})^{T}

Example of a Matrix Inverse.

A =

1 1 0

2 2 1

0 1 1

A^{-1} =

-1 1 -1

2 -1 1

-2 1 0

Computing the Inverse of a Matrix.

The Gauss-Jordan method permits one to calculate the inverse of a matrix. Begin by writing the matrix A and the identity matrix I side by side as shown in Tableau^{0} below.

Using the following row operations simultaneously on both matrices in the tableau, convert the left hand matrix to the identity matrix. When this is done, the right hand matrix is the inverse matrix.

1. exchange two rows.

2. divide or multiply a row by a constant, called the pivot.

3. add a multiple of one row to another row.

Begin in the upper left hand corner of the left hand matrix. Proceed along the diagonal converting each column into the appropriate column of the identity matrix as shown in the sequence of tableaux. Only exchange rows if the diagonal entry is zero. Moreover, if the diagonal entry is zero, exchange rows with a lower row that has a nonzero entry in that column. If no suitable row exists, the matrix does not have an inverse.

Tableau^{0}

A

I

R^{0}_{1}

R^{0}_{2}

R^{0}_{3}

1

1

0

2

2

1

0

1

1

1

0

0

0

1

0

0

0

1

Tableau^{1}

R^{1}_{1} = R^{0}_{1} / 1

R^{1}_{2} = R^{0}_{2} - (2) * R^{1}_{1}

R^{1}_{3} = R^{0}_{3} - (0) * R^{1}_{1}

1

1

0

0

0

1

0

1

1

1

0

0

-2

1

0

0

0

1

Tableau^{2}

Switch R3 and R2

R^{2}_{1} = R^{1}_{1} - (1) * R^{2}_{2}

R^{2}_{2} = R^{1}_{2} / 1

R^{2}_{3} = R^{1}_{3} - (0) * R^{2}_{2}

1

0

-1

0

1

1

0

0

1

1

0

-1

0

0

1

-2

1

0

Tableau^{3}

I

A^{-1}

R^{3}_{1} = R^{2}_{1} - (-1) * R^{3}_{3}

R^{3}_{2} = R^{2}_{2} - (1) * R^{3}_{3}

R^{3}_{3} = R^{2}_{3} / 1

1

0

0

0

1

0

0

0

1

-1

1

-1

2

-1

1

-2

1

0

Use the form below to compute the inverse of some other 3 by 3 A matrix.

Matrix Determinant.

The determinant of a square matrix A, det(A), is a scalar calculated as follows:

3. If A is an n by n matrix, the determinant M_{i, j} of the matrix formed by deleting A's ith row and jth column is called the minor determinant for a_{i, j}.

4. The cofactor A_{i, j} = (-1)^{(i+j)} * M_{i, j}.

5. The determinant det(A) = a_{i, 1} * A_{i, 1} + a_{i, 2} * A_{i, 2} + . . . . . + a_{i, n } * A_{i, n}, for any i = 1, . . . n. (Expanding by any row of A).

5. The determinant det(A) = a_{1, j} * A_{1, j} + a_{2, j} * A_{2, j} + . . . . . + a_{n, j } * A_{n, j}, for any j = 1, . . . n. (Expanding by any column of A).

For example, compute det(A) by expanding by the 2nd row:

A =

1 1 2

2 3 0

0 1 1

det(A) =

(-1)^{(2+1)} * 2 * det(

1 2

1 1

) + (-1)^{(2+2)} * 3 * det(

1 2

0 1

) + (-1)^{(2+3)} * 0 * det(

1 1

0 1

) = 5

Properties of Determinants.

If A is square n by n matrix:

1. det(A) = det(A^{T}).

2. det(A * B) = det(A) * det(B), where B is any square n by n matrix.

3. det(A) != 0 if and only if A has an inverse. If det(A) = 0, A is called a singular matrix.

4. If A^{-1} exists, det(A) = 1 / det(A^{-1}).

5. det(A) is unchanged if a multiple of one row is added to another row; det(A) is unchanged if a multiple of one column is added to another column.

6. det(A) changes sign if two rows are switched; det(A) changes sign if two columns are switched.

7. Multiplying a column or row by a constant scalar multiplies the determinant by the constant scalar.

Eigenvalues and Eigenvectors.

The scalar and vector pair, (µ, v), is called an eigenvalue and eigenvector pair of the square matrix A if:

A *v = µ * v, or

(µ * I - A ) * v = 0.

Characteristic Equation.

The eigenvalues of A can be found by finding the roots of the polynomial in µ:

f(µ) = det(µ * I - A) = 0.

Computing Eigenpairs: Introduction.

For example, compute the eigenvalues and eigenvectors for the matrix:

The polynomial f(µ) has roots µ_{1} = 4, and µ_{2} = 2.

Eigenvectors.

Each eigenvalue µ_{k} provides a system of equations (µ_{k} * I - A ) * v = 0 that will yield the eigenvector(s) v_{k} associated with eigenvalue µ_{k} .

For the eigenvalue µ_{1} = 4:

(µ_{1} * I - A ) * v = (4 * I - A ) * v =

-1

1

-3

3

*

v_{1}

v_{2}

=

- v_{1} + v_{2} = 0

- 3*v_{1} + 3*v_{2} = 0

e_{1}

e_{2}

Equations e_{1} and e_{2} are the same, after dividing equation e_{2} by 3. Solving for v_{2} in terms of v_{1} yields:

v_{2} = v_{1}.

This means that the vector v_{1}^{T} = (1, 1), and any scalar multiple s * v_{1}, are eigenvectors associated with the eigenvalue µ_{1} = 4.

The set of eigenvectors associated with the eigenvalue µ_{1} = 4 is the kernel of the matrix A_{1}, ker(A_{1}), where A_{1} = (4 * I - A ).

For the eigenvalue µ_{2} = 2:

(µ_{2} * I - A ) * v = (2 * I - A ) * v =

-3

1

-3

1

*

v_{1}

v_{2}

=

- 3*v_{1} + v_{2} = 0

- 3*v_{1} + v_{2} = 0

e_{1}

e_{2}

Equations e_{1} and e_{2} are the same. Solving for v_{2} in terms of v_{1} yields:

v_{2} = 3 * v_{1}.

This means that the vector v_{2}^{T} = (1, 3), and any scalar multiple s * v_{2}, are eigenvectors associated with the eigenvalue µ_{2} = 2.

The set of eigenvectors associated with the eigenvalue µ_{2} = 2 is the kernel of the matrix A_{2}, ker(A_{2}), where A_{2} = (2 * I - A ).

The next diagram shows v_{1}, A * v_{1}, v_{2}, A * v_{2}, u = (1.5, -2), and A * u. See what happens if you set the vector u equal to a scalar multiple of one of the eigenvectors.

You can change the u vector in the form above to see how u and A * u change. For example, try u = (-1.5, 2) and click "Eigenvector".

Eigenpair Decomposition of A.

The matrix A has two unique eigenvalues and two linearly independent eigenvectors. These eigenvalues and eigenvectors can be used to decompose (factor) the matrix A as follows. Let S be the matrix whose first column is the eigenvector v_{1}, and whose second column is the eigenvector v_{2}. Let D be the diagonal matrix formed from the eigenvalues of A in the order of µ_{1}, µ_{2}.

S =

1

1

1

3

D =

4

0

0

2

By construction of S and D, A * S = S * D. Since the columns of S are linearly independent, S has an inverse matrix, S^{(-1)}. Multiplying A * S = S * D on the right by S^{(-1)} yields the eigenpairs decomposition of the matrix A as A = S * D * S^{(-1)}.

A = S * D * S^{(-1)} =

1

1

1

3

*

4

0

0

2

*

1.5

-0.5

-0.5

0.5

=

5

-1

3

1

The Determinant and Trace of the matrix A.

The determinant of A, det(A), equals the product of the eigenvalues of A.

det(A) = µ_{1} * µ_{2} = 4 * 2 = 8.

Consequently, the determinant of A is nonzero if and only if A has no zero eigenvalue. Furthermore, the matrix A is invertible if and only if A has no zero eigenvalue.

The trace of A, trace(A), equals the sum of the eigenvalues of A.

trace(A) = µ_{1} + µ_{2} = 4 + 2 = 6.

The general 2 by 2 matrix:

A =

a

b

c

d

has trace(A) = a + d, and det(A) = a*d - b*c. Moreover, its characteristic equations is:

det(

µ - a

b

c

µ - d

)

= µ^{2} - trace(A) * µ + det(A)

Factoring this quadratic equation yields the pair of eigenvalues, µ_{1} and µ_{2}.

Jordan Normal Form.

A matrix with repeated eigenvalues has an eigenpair decomposition if it has a complete set of linearly independent eigenvectors. Otherwise, the Jordan Normal Form augmented with generalized eigenvectors is used to decompose the matrix.

with eigenvalues µ_{1} = 2, and µ_{2} = 2. However, the only eigenvector is v^{T}_{1} = (1, -1) and its multiples.

Recall that the eigenvector v^{T}_{1} = (1, -1) is computed by solving the system of equations (2 * I - A ) * v = 0:

(µ_{1} * I - A ) * v = (2 * I - A ) * v =

-1

-1

1

1

*

v_{1}

v_{2}

=

- 1*v_{1} - 1*v_{2} = 0

v_{1} + v_{2} = 0

e_{1}

e_{2}

Equations e_{1} and e_{2} are the same. Solving for v_{2} in terms of v_{1} yields:

v_{2} = -v_{1}.

This means that the vector v_{1}^{T} = (1, -1) (and its scalar multiples) is the only eigenvector associated with the eigenvalues µ_{1,2} = 2.

Computing Generalized Eigenvectors: Introduction.

One can find a generalized eigenvalue v_{2} for the eigenvalues µ_{1,2} = 2 by solving the system of equations:

(A - 2 * I)*v_{2} = v_{1}, or

(2 * I - A )*v_{2} = - v_{1}.

Expanding:

(2 * I - A )* v = -v_{1}

-1

-1

1

1

*

-v_{1}

-v_{2}

=

- 1*v_{1} - 1*v_{2} = -1

v_{1} + v_{2} = 1

e_{1}

e_{2}

Adding equation e_{1} to equation e_{2}, and multiplying e_{1} by -1, yields the echelon form:

1*v_{1} + 1*v_{2} = 1

0 * v_{1} + 0 *v_{2} = 0

e_{1}

e_{2}

The particular solution v^{T}_{p} = (1, 0) = v_{2} is the generalized eigenvalue of the matrix A associated with the pair of equal eigenvalues µ_{1,2} = 2, while the homogeneous solution v_{h} is a multiple of the eigenvector v_{1}:

v

=

v_{p}

+ s *

v_{h}

v_{1}

v_{2}

=

1

0

+ s *

-1

1

I obtained this result by using the echelon algorithm form on the linear equations web page to solve the system of equations.

As an eigenvector v_{1} satisfies:

(A - 2 * I)*v_{1} = o, or A * v_{1} = 2 * v_{1}.

As the generalized eigenvector constructed from v_{1}, v_{2} satisfies:

(A - 2 * I)*v_{2} = v_{1}, or A * v_{2} = 2 * v_{2} + 1 * v_{1}.

Jordan Decomposition.

Let V be the matrix whose first column is the eigenvector v_{1}, and whose second column is the generaralized eigenvector v_{2}. Let J be the upper triangular matrix whose diagonal consists of eigenvalues of A, µ_{1,2} = 2, and with a 1 in the upper right hand corner.

V =

1

1

-1

0

J =

2

1

0

2

By construction of v_{1} and v_{2}, A * V = V * J. Because v_{1} and v_{2} are linearly independent, the matrix V has an inverse. The Jordan Decomposition is:

The eigenpair and Jordan decompositions of the matrices above are just two of the ways that matrices can be factored. Matrices can be decomposed into the product of other matrices, such as triangular, diagonal and/or orthogonal matrices. The decomposed matrices can be used to solve systems of linear equations.

The methods available to decompose a matrix depends on the dimensions of the matrix and on the rank of the matrix.
A general matrix that is square and of full rank, ie its rows and columns are linearly independent, can be factored using Gaussian or Gram-Schmidt Decomposition. Matrices that are positive definite, symmetric, or triangular can be factored by more specific methods.

An m by n matrix with m > n can be factored with Householder or Jacobi transformations, or through singular value decomposition (svd).

An m by n matrix with m < n can be factored through singular value decomposition (svd).

Factor A Square Matrix of Full Rank

Gauss Matrix Decomposition.

The objective of Gaussian Decomposition is to write a matrix A as the product of a lower triangular matrix L and an upper triangular matrix U.

The Gaussian method is similar to the Gauss-Jordan method for computing a matrix's inverse. However, in each tableau the Gaussian method searches for the pivot with the largest absolute value and switches rows in both matrices accordingly. Moreover, it reduces the LHS matrix of the tableau to an upper triangular matrix U with the pivots along its diagonal, and the RHS matrix to a lower triangular matrix L with zeros along its diagonal, storing the multipliers used to zero out the U matrix below its diagonal. Furthermore, it keeps track of the row switches to construct a permutation matrix P. The actual L matrix is obtained from the RHS matrix of the final tableau by placing ones along its diagonal.

Gauss matrix decomposition yields matrices P, L, and U such that:

P * A = L * U, or

A = P^{T} * L * U.

Factor a matrix A so that: P * A = L *U.

Tableau^{0}

A

L

R^{0}_{1}

R^{0}_{2}

R^{0}_{3}

R^{0}_{4}

R^{0}_{5}

55

56

59

65

75

56

58

62

69

80

59

62

68

77

90

65

69

77

89

105

75

80

90

105

125

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

Tableau^{1}

U

L

Switch R1 and R5

R^{1}_{1}

R^{1}_{2} = R^{0}_{2} - (0.7467) * R^{1}_{1}

R^{1}_{3} = R^{0}_{3} - (0.7867) * R^{1}_{1}

R^{1}_{4} = R^{0}_{4} - (0.8667) * R^{1}_{1}

R^{1}_{5} = R^{0}_{5} - (0.7333) * R^{1}_{1}

75

80

90

105

125

0

-1.7333

-5.2

-9.4

-13.3333

0

-0.9333

-2.8

-5.6

-8.3333

0

-0.3333

-1

-2

-3.3333

0

-2.6667

-7

-12

-16.6667

0

0

0

0

0

0.7467

0

0

0

0

0.7867

0

0

0

0

0.8667

0

0

0

0

0.7333

0

0

0

0

Tableau^{2}

U

L

Switch R2 and R5

R^{2}_{1}

R^{2}_{2}

R^{2}_{3} = R^{1}_{3} - (0.35) * R^{2}_{2}

R^{2}_{4} = R^{1}_{4} - (0.125) * R^{2}_{2}

R^{2}_{5} = R^{1}_{5} - (0.65) * R^{2}_{2}

75

80

90

105

125

0

-2.6667

-7

-12

-16.6667

0

0

-0.35

-1.4

-2.5

0

0

-0.125

-0.5

-1.25

0

0

-0.65

-1.6

-2.5

0

0

0

0

0

0.7333

0

0

0

0

0.7867

0.35

0

0

0

0.8667

0.125

0

0

0

0.7467

0.65

0

0

0

Tableau^{3}

U

L

Switch R3 and R5

R^{3}_{1}

R^{3}_{2}

R^{3}_{3}

R^{3}_{4} = R^{2}_{4} - (0.1923) * R^{3}_{3}

R^{3}_{5} = R^{2}_{5} - (0.5385) * R^{3}_{3}

75

80

90

105

125

0

-2.6667

-7

-12

-16.6667

0

0

-0.65

-1.6

-2.5

0

0

0

-0.1923

-0.7692

0

0

0

-0.5385

-1.1538

0

0

0

0

0

0.7333

0

0

0

0

0.7467

0.65

0

0

0

0.8667

0.125

0.1923

0

0

0.7867

0.35

0.5385

0

0

Tableau^{4}

U

L

Switch R4 and R5

R^{4}_{1}

R^{4}_{2}

R^{4}_{3}

R^{4}_{4}

R^{4}_{5} = R^{3}_{5} - (0.3571) * R^{4}_{4}

75

80

90

105

125

0

-2.6667

-7

-12

-16.6667

0

0

-0.65

-1.6

-2.5

0

0

0

-0.5385

-1.1538

0

0

0

0

-0.3571

0

0

0

0

0

0.7333

0

0

0

0

0.7467

0.65

0

0

0

0.7867

0.35

0.5385

0

0

0.8667

0.125

0.1923

0.3571

0

P

A

0

0

0

0

1

1

0

0

0

0

0

1

0

0

0

0

0

1

0

0

0

0

0

1

0

55

56

59

65

75

56

58

62

69

80

59

62

68

77

90

65

69

77

89

105

75

80

90

105

125

L

U

1

0

0

0

0

0.7333

1

0

0

0

0.7467

0.65

1

0

0

0.7867

0.35

0.5385

1

0

0.8667

0.125

0.1923

0.3571

1

75

80

90

105

125

0

-2.6667

-7

-12

-16.6667

0

0

-0.65

-1.6

-2.5

0

0

0

-0.5385

-1.1538

0

0

0

0

-0.3571

Notice that in each tableau the relevant information in the L and U matrices fit onto one 5 by 5 matrix, which is what computer packages do to save storage space.

Compute the Determinant of A with Gauss Decomposition.

With the P * A = L * U decomposition the determinant of A, det(A), is the product of the diagonal elements of the matrix U, times +1 or -1 if an even or an odd number of row switches were used in the decomposition.

Use the form below to factor some other square matrix A, up to 5 by 5.

If your matrix is smaller, specify its dimension and fill in / change the numbers in the upper left hand corner, and then click "Factor".

Gaussian Decomposition.

Gram-Schmidt Matrix Decomposition.

Given a square matrix A of dimension 4, the objective is to decompose A into two matrices Q and R, where Q is an orthonormal matrix, and R is an upper triangular matrix, such that:

A = Q * R.

Write A as consisting of 4 column vectors.

A =

a_{1}

a_{2}

a_{3}

a_{4}

=

1

2

1

3

2

1

3

2

1

2

3

2

4

1

7

4

The Gram-Schmidt process is based on the method used on the vectors web page. The objective is to construct a set of 4 orthonormal vectors {
q_{1}, q_{2}, q_{3}, q_{4}}, from the set of 4 linearly independent vectors {
a_{1}, a_{2}, a_{3}, a_{4}}.

Begin by writing q_{1} = a_{1}. Consider a_{2}. The projection of a_{2} on q_{1} is the vector:

is orthogonal to q_{1}, and {q_{1}, q_{2}} span the same vector subspace as {a_{1}, a_{2}}.

Repeat this process with the remaining {
a_{3}, a_{4}} vectors. At stage k for a general matrix A, subtract the projections of a_{k} on the orthogonal vectors {q_{1}, ... q_{k-1}} from a_{k} to construct the next orthogonal vector q_{k}.

Do this at each stage k, by constructing new projection vectors, p_{j}, obtained by projecting a_{k} on each q_{j} as given by:

The required vector q_{k} orthogonal to {q_{1}, ... q_{k-1}} is:

q_{k} = a_{k} - p_{1} - p_{2} - . . . - p_{k-1}

The calculations for each step of the Gram-Schmidt decomposition for the example matrix are displayed below. The columns of the orthogonal matrix are divided by their norms to produce an orthonormal matrix.

q_{1}

=

1

2

1

4

q_{2}

=

2

1

2

1

- (0.4545 ) *

1

2

1

4

=

1.54545

0.09091

1.54545

-0.81818

q_{3}

=

1

3

3

7

- (1.7273 ) *

1

2

1

4

- (0.1333 ) *

1.54545

0.09091

1.54545

-0.81818

=

-0.93333

-0.46667

1.06667

0.2

q_{4}

=

3

2

2

4

- (1.1364 ) *

1

2

1

4

- (0.85 ) *

1.54545

0.09091

1.54545

-0.81818

- (-0.3529 ) *

-0.93333

-0.46667

1.06667

0.2

=

0.22059

-0.51471

-0.07353

0.22059

Orthogonal Matrix Q

1

1.5455

-0.9333

0.2206

2

0.0909

-0.4667

-0.5147

1

1.5455

1.0667

-0.0735

4

-0.8182

0.2

0.2206

Orthonormal Matrix Q

0.2132

0.6617

-0.6199

0.3638

0.4264

0.0389

-0.31

-0.8489

0.2132

0.6617

0.7085

-0.1213

0.8528

-0.3503

0.1328

0.3638

Gram-Schmidt Decomposition A = Q * R

Orthonormal matrices have the property that Q^{T} * Q = I, the identity matrix. With A and Q known matrices, compute R as follows:

Q^{T} * A = Q^{T} * Q * R = R.

Expanding this last equation yields:

R =

(q_{1}^{T} * a_{1})

(q_{1}^{T} * a_{2})

(q_{1}^{T} * a_{3})

(q_{1}^{T} * a_{4})

(q_{2}^{T} * a_{1})

(q_{2}^{T} * a_{2})

(q_{2}^{T} * a_{3})

(q_{2}^{T} * a_{4})

(q_{3}^{T} * a_{1})

(q_{3}^{T} * a_{2})

(q_{3}^{T} * a_{3})

(q_{3}^{T} * a_{4})

(q_{4}^{T} * a_{1})

(q_{4}^{T} * a_{2})

(q_{4}^{T} * a_{3})

(q_{4}^{T} * a_{4})

=

(q_{1}^{T} * a_{1})

(q_{1}^{T} * a_{2})

(q_{1}^{T} * a_{3})

(q_{1}^{T} * a_{4})

(q_{2}^{T} * a_{2})

(q_{2}^{T} * a_{3})

(q_{2}^{T} * a_{4})

(q_{3}^{T} * a_{3})

(q_{3}^{T} * a_{4})

(q_{4}^{T} * a_{4})

The entries below the diagonal of R are zero, because a_{j} is orthogonal to q_{k} for j < k.

The A = Q * R factorization for the example matrix is:

Matrix A

1

2

1

3

2

1

3

2

1

2

3

2

4

1

7

4

=

Orthonormal Matrix Q

0.2132

0.6617

-0.6199

0.3638

0.4264

0.0389

-0.31

-0.8489

0.2132

0.6617

0.7085

-0.1213

0.8528

-0.3503

0.1328

0.3638

*

Matrix R

4.6904

2.132

8.1016

5.33

0

2.3355

0.3114

1.9852

0

0

1.5055

-0.5314

0

0

0

0.6063

Use the form below to factor some other square matrix A, up to 5 by 5.

If your matrix is smaller, specify its dimension and fill in / change the numbers in the upper left hand corner, and then click "Factor".

Gram-Schmidt Decomposition.

Factor a Matrix with # Rows m >= # Columns n

Matrices with more rows than columns are found in least-squares problems. Three methods to decompose m by n matrices with m >= n are Householder, Jacobi, and singular value decomposition. Below I describe Householder and Jacobi transformations, leaving singular value decompositions for the case where m < n, ie where A has more columns than rows.

Householder Decomposition

Householder decomposition is used on the Statistics web pages to solve least-squares (multiple regression) problems.

Given an m by n matrix A, with rank(A) = n and with m >= n, factor A into the product of an orthonormal matrix Q and an upper triangular matrix R, such that:

A = Q * R.

The Householder decomposition is based on the fact that for any two different vectors, v and w, with ||v|| = ||w||, (ie different vectors of equal length), a reflection matrix H exists such that:

H * v = w.

To obtain the matrix H, define the vector u by:

u = (v - w) / (||v - w||).

The matrix H defined by:

H = I - 2 * u * u^{T}

is the required reflection matrix.

Observe that in the following 3-dimensional diagram,

q = v - w,

p = (1/2)*(v + w)

v = (1/2)* [(v - w) + (v + w)]

u = (v - w) / ||v - w|| = q / ||q||

u^{T} * (2 * p) = 0, since they are perpendicular

To prove that H * v = w, observe that

H * v = (I - 2 * u * u^{T}) * v = v - 2 * u * (u^{T} * v), and

2 * u* (u^{T} * v) = u * (u^{T} * [(v - w) + (2*p)]) = (v - w) * [ (v - w)^{T}/ ||(v - w)||] * [ (v - w)/ ||(v - w)||] = v - w.

Consequently,
H * v = v - (v - w) = w.

Note that H = I - 2 * P, where P = u * u^{T} is the projection matrix onto the yellow line in the diagram above, along the vector (v + w).

Factor the following matrix as A = Q * R.

A =

16

2

3

5

11

10

9

7

6

4

14

15

Given the matrix A, of dimension 4 by 3, the idea is to turn A into an upper triangular matrix using a sequence of 3 Householder transformations,
H_{1}, H_{2}, H_{3}.

Starting with A_{1} = A, let v_{1} = the first column of A_{1}, and w_{1}^{T} = (norm(v_{1}), 0,...0), ie a column vector whose first component is the norm of v_{1} with the remaining components equal to 0. The Householder transformation H_{1} = I - 2 * u_{1} * u_{1}^{T} with u_{1} = v_{1} - w_{1} / ||v_{1} - w_{1}|| will turn the first column of A_{1} into w_{1} as with H_{1} * A1 = A_{2}. At each stage k, v_{k} = the kth column of A_{k} on and below the diagonal with all other components equal to 0, and w_{k}'s kth component equals the norm of v_{k} with all other components equal to 0. Letting H_{k} * A_{k} = A_{k+1}, the components of the kth column of A_{k+1} below the diagonal are each 0. These calculations are listed below for each stage for the matrix A.

Use the form below to factor some other m by n matrix A, up to 5 by 4.

If your matrix is smaller, specify its dimension and fill in / change the numbers in the upper left hand corner, and then click "Factor".

Householder Decomposition.

Jacobi Decomposition

Jacobi (Givens Rotations) transformations are used to factor an m by n matrix A, with rank(A) = n and with m >= n, into the product of an orthonormal matrix Q and an upper triangular matrix R, such that:

A = Q * R.

Householder transformation are preferred to Jacobi transformations, because Householder transformations require fewer (almost one half) calculations.

A Jacobi transformation, J, is based on a rotation matrix, where the counter-clockwise angle of rotation, Ø, is selected to rotate a vector v^{T} = (v_{1}, v_{2}) into a vector w^{T} = (w_{1}, 0).

These values for the matrix J are computed directly from the vector v, without computing the angle Ø.

Starting in the upper left hand corner of a matrix A, such as:

A =

16

2

3

5

11

10

9

7

6

4

14

15

zero out all entries below the diagonal of A in the first column with one Jacobi transformation per entry. Similarly, zero out the entries below the diagonal in the other columns of A.

Use the form below to factor some other m by n matrix A, up to 5 by 4.

If your matrix is smaller, specify its dimension and fill in / change the numbers in the upper left hand corner, and then click "Factor".

Jacobi Decomposition.

Factor a Matrix with # Rows m <= # Columns n

Matrices with fewer rows than columns are found in under-determined problems, with fewer equations than unknown variables. Singular value decomposition is a convenient way to factor the matrix A where m < n, ie where A has more columns than rows.

However, any matrix can be factored into its singular value components. Furthermore, the matrices from the singular value decomposition can be combined to form the pseudoinverse (Moore-Penrose generalized inverse) of an arbitrary m by n matrix. If A is a square matrix of full rank, its pseudoinverse equals the matrix inverse. If A has independent columns, the pseudoinverse equals its left-inverse. If A has independent rows, its pseudoinverse equals its right-inverse.

Singular Value Decomposition

The Singular Value Decomposition (SVD) of a matrix A, m by n, is based on the following theorem of linear algebra:

Any m by n matrix A can be decomposed (factored) into:

A = U * S* V^{T}, where

- the m by m orthogonal matrix U consists of the eigenvectors of A * A^{T}.

- the n by n orthogonal matrix V consists of the eigenvectors of A^{T} * A.

- the m by n diagonal matrix S consists of the square roots of the eigenvalues of A * A^{T} (which equal the eigenvalues of A^{T} * A), arranged in descending order.

- the eigenvectors of A * A^{T} and A^{T} * A are arranged in the columns of U and V, respectively, in order of their shared eigenvalues^{1/2} along the diagonal of S.

- the diagonal entries of S are called the singular values of A, all of which are nonnegative.

- the number of positive singular values equals the rank of the matrix A.

The Pseudoinverse A^{+}

A^{+} = V * S^{+} * U^{T}

- the n by m diagonal matrix S^{+}consists of the reciprocals of the positive singular values on the diagonal of S, computed in the same order.

Example of a Singular Value Decomposition.

A = U * S * V^{T}

A

=

U

*

S

*

V^{T}

1

1

1

1

1

2

3

4

1

3

6

10

=

-0.1274

0.7453

-0.6544

-0.4061

0.5628

0.72

-0.9049

-0.3574

-0.231

*

13.341

0

0

0

0

1.3997

0

0

0

0

0.2395

0

*

-0.1078

-0.2739

-0.5078

-0.8096

0.6792

0.5705

0.2065

-0.413

-0.6907

0.3866

0.4995

-0.3521

-0.2236

0.6708

-0.6708

0.2236

The Pseudoinverse of A.

A^{+} = V * S^{+} * U^{T}

A^{+}

=

V

*

S^{+}

*

U^{T}

2.25

-1.8

0.5

-0.75

1.4

-0.5

-1.25

1.6

-0.5

0.75

-1.2

0.5

=

-0.1078

0.6792

-0.6907

-0.2236

-0.2739

0.5705

0.3866

0.6708

-0.5078

0.2065

0.4995

-0.6708

-0.8096

-0.413

-0.3521

0.2236

*

0.075

0

0

0

0.7145

0

0

0

4.1754

0

0

0

*

-0.1274

-0.4061

-0.9049

0.7453

0.5628

-0.3574

-0.6544

0.72

-0.231

Use the form below to factor some other m by n matrix A, up to 4 by 5.

If your matrix is smaller, specify its dimension and fill in / change the numbers in the upper left hand corner, and then click "Factor".

Singular Value Decomposition.

Singular value decomposition works on a matrix of any dimensions, regardless of its rank.

References.

Ayres, Frank Jr. Matrices. New York: Schaum McGraw-Hill, 1962.

Ayres, Frank Jr. Modern Algebra. New York: Schaum McGraw-Hill 1965.

Bretscher, Otto. Linear Algebra with Applications. Upper Saddle River: Prentice Hall, 1997.

Burden, Richard L. and J. Douglas Faires. Numerical Analysis. 6th ed. Pacific Grove: Brooks/Cole, 1997.

Campbell, Hugh D. Matrices with Applications. New York: Appleton 1968.

Cohn, P. M. Linear Equations. London: Routledge, 1964.

Demmel, James W. Applied Numerical Linear Algebra. Philadelphia: Siam, 1997.

Dowling, Edward T. Mathematics for Economists. New York: Schaum McGraw-Hill, 1980.

Goult, R. J., R. F. Hoskins, J.A. Milner, and M. J. Pratt. Computational Methods in Linear Algebra. New York: John Wiley, 1974

Lipschutz, Seymour. Linear Algebra. New York, Schaum McGraw-Hill, 1968.

Mathews, John H. and Kurtis D. Fink. Numerical Methods Using MATLAB. 3rd ed. Upper Saddle River: Prentice Hall, 1999.

Press, William H., Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling. Numerical Recipes: The Art of Scientific Computing. Cambridge: Cambridge UP, 1989.

Spiegel, Murray R. Vector Analysis and an Introduction to Tensor Analysis. New York: Schaum, 1959.

Strang, Gilbert. Linear Algebra and Its Applications. 3d ed. San Diego: Harcourt, 1976.

Varah, James. Numerical Linear Algebra: Computer Science 402 Lecture Notes. Vancouver: University of B.C., 2000.

Watkins, David S. Fundamentals of Matrix Computations. New York: John Wiley, 1991.