     Egwald Mathematics: Linear Algebra

Matrices

by

Elmer G. Wiens

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: Definition of a Matrix.

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

 A = a1,1 a1,2     . . . . .    a1,n-1 a1,n a2,1 a2,2     . . . . .    a2,n-1 a2,n               . . . . . . . . . . .           am,1 am,2     . . .     am,n-1 am,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 matrix O = [oi, j], with oi, j = 0 for all i and j.

The Identity Matrix.

A square matrix is an identity matrix I = [ii, j], if  ii, j = 1 where i = j, and ii, 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 * a1,1 s * a1,2     . . . . .    s * a1,n s * a2,1 s * a2,2     . . . . .    s * a2,n               . . . . . . . . . . .                s * am,1 s * am,2     . . .     s * am,n

or s * A = [s * ai, j], with each entry multiplied by the number s.

Given matrices A = [ai, j] and B = [bi, j ] of dimension m by n, their sum A + B is another matrix C = [ci, j ] of dimension m by n, where:

ci, j = ai, j + bi, 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 = [ai, j] of dimension m by r and B = [bi, j ] of dimension r by n, their product A * B is another matrix C = [ci, j ] of dimension m by n, where:

ci, j = ai, 1 * b1, j + ai, 2 * b2, j + . . . . . ar, 1 * b1, r , for i = 1, . . . m; j = 1, . . . n

For example,

 1   3   -2   3 * 8   -1   4;  2 = (8 + 12)    (-1 + 6) (-16 + 12)    (2 + 6) = 20   5  -4   8

Two matrices are conformable for multiplication if their dimensions permit multiplication.

Two matrices are commutative under multiplication if A * B = B * A.

Product of a Matrix and a Vector.

The product of a column vector x = [xi] and a matrix A = [ai, j] is a column vector z, where :

z = A * x.

For example,

 8 1 6 3 5 7 4 9 2
*
 1 1 1
=
 8 + 1 + 6 3 + 5 + 7 4 + 9 + 7
=
 15 15 15

The product of a row vector y = [yi] and a matrix A = [ai, j] is a row vector w, where :

w = y * A.

For example,

 1 1 1
*
 8 1 6 3 5 7 4 9 2
=
 (8 + 3 + 4) (1 + 5 + 9) (6 + 7 + 2)
=
 15 15 15

The Transpose of a Matrix.

The transpose of an m by n matrix A = [ai, j] is an n by m matrix AT = [aj, i] formed by interchanging the rows and columns of A.

 A = 10  5   7  9  22  7   5  6  3   15  4  4 AT = 10  22  3   5   7   15  7   5   4   9   6   4

The transpose of a column vector x is a row vector xT, and vice versa.

From this point, vectors denoted by underlined small case letters, x, are column vectors, while row vectors have the superscript T, xT.

Given matrices A = [ai, j] of dimension m by r and B = [bi, j ] of dimension r by n, the transpose of their product obeys:

(A * B)T = BT * AT.

Vector Operations.

Vectors and their operations are described on the linear algebra: vectors web page.

Inner Product.

The inner product of a column vector x = [xi] (n by 1) and a row vector yT = [yi] (1 by n) is their product as matrices:

yT * x = y1 * x1 + y2 * x2 + . . . . + yn * xn

For example,

 2 4 -3
*
 3 -6 5
= 2*3 + 4*(-6) + (-3)*5 = -33

Outer Product.

The outer product of a column vector x = [xi] (n by 1) and a row vector yT = [yi] (1 by n) is a matrix Z = [zi, j] = x * yT, where zi, j = xi * yj.

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 = [xi] is the square root of  xT * x, ie ||x|| = (xT * x)1/2.

For example, if vT = (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 = AT.

Given matrices A and ATdefined above, then the matrices:

 A * AT = 255   344   169  344   594   215  169   215   266 AT * 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 = -AT.

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 = [ai, j] is tr(A) = a1, 1 + a2, 2 . . . . . + an, 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 Rn to R is called a quadratic form if it is defined by:

f(x) = xT * A * x,   for all vectors x in Rn.

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 Rn. ie

f(x) = xT * A * x > 0,   for all nonzero vectors x in Rn.

Example.

The identity matrix I is positive definite because xT *I * x = x12 + x22 + . . . + xn2 > 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 Rn, the matrix norm induced by ||x|| is defined by:

||A||  =  max{ ||A * x||  |   all x in Rn with ||x|| = 1 }

Properties of Matrix Norms.

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

 ||A|| = 0 if and only if A = O, the zero matrix. ||A|| > 0 for any non-zero matrix A. ||s * A|| = |s| * ||A||. ||A + B|| <= ||A|| + ||B||. ||A * B|| <= ||A|| * ||B||. ||A * x|| <= ||A|| * ||x||, for induced matrix norms.

The Frobenius Norm.

The Frobenius Norm, ||.||F looks much like the Euclidean vector norm, but it is not the induced norm of the Euclidean norm.

For a given square matrix A of dimension n:

||A||F = ( |a1,1| + |a1,2| + . . . . + |a1,n|) + ( |a2,1| + |a2,2| + . . . . + |a2,n|) + .   .   .   . + ( |an,1| + |an,2| + . . . . + |an,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:

Q0.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. Rotation Matrix Q0.52
v1: v2:
-2 * π <= Ø <= 2 * π

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-Ø = Q0 QØ * QØ = Q2*Ø 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:

Q0.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. Projection Matrix P0.524
v1: v2:
-2 * π <= Ø <= 2 * π

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:

H0.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. Reflection Matrix H0.79
v1: v2:
-2 * π <= Ø <= 2 * π

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 ai of A,  aiT * aj = 0 for any other column aj of A.   If each column ai of A has a norm of 1, then A is orthonormal. If A is orthonormal, then:

AT * 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 AT = A-1  if A is orthonormal A-1 exists if and only if rank(A) = n. (AT)-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 Tableau0 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.

Tableau0
AI
 R01 R02 R03
 1 1 0 2 2 1 0 1 1
 1 0 0 0 1 0 0 0 1

Tableau1
 R11 = R01 / 1 R12 = R02 - (2) * R11 R13 = R03 - (0) * R11
 1 1 0 0 0 1 0 1 1
 1 0 0 -2 1 0 0 0 1

Tableau2
Switch R3 and R2
 R21 = R11 - (1) * R22 R22 = R12 / 1 R23 = R13 - (0) * R22
 1 0 -1 0 1 1 0 0 1
 1 0 -1 0 0 1 -2 1 0

Tableau3
IA-1
 R31 = R21 - (-1) * R33 R32 = R22 - (1) * R33 R33 = R23 / 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 A

Matrix Determinant.

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

1. If A = [a], a 1 by 1 matrix, det(A) = a.

2. If A =
 a1, 1   a1, 2 a2, 1   a2, 2
det(A) = a1, 1 * a2, 2 - a2, 1 * a1, 2.
3. If A is an n by n matrix, the determinant Mi, j of the matrix formed by deleting A's ith row and jth column is called the minor determinant for ai, j.
4. The cofactor Ai, j = (-1)(i+j) * Mi, j.
5. The determinant det(A) = ai, 1 * Ai, 1 + ai, 2 * Ai, 2 + . . . . . + ai, n * Ai, n,   for any i = 1, . . . n. (Expanding by any row of A).
5. The determinant det(A) = a1, j * A1, j + a2, j * A2, j + . . . . . + an, j * An, 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(AT). 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:

A =
 5 -1 3 1

Eigenvalues.

The characteristic equation is the polynomial:

f(µ) = det
 ( (µ - 5) ; 1 ) -3 (µ - 1)
= (µ - 5) * (µ - 1) - (-3) * (1)   =   µ2 - 6 * µ + 8   =   (µ - 4) * (µ - 2)

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) vk associated with eigenvalue µk .

For the eigenvalue µ1 = 4:

1 * I - A ) * v   =   (4 * I - A ) * v   =
 -1 1 -3 3
*
 v1 v2
=
 - v1 + v2   =   0 - 3*v1 + 3*v2   =   0
 e1 e2

Equations e1 and e2 are the same, after dividing equation e2 by 3. Solving for v2 in terms of v1 yields:

v2 = v1.

This means that the vector v1T = (1, 1), and any scalar multiple s * v1, are eigenvectors associated with the eigenvalue µ1 = 4.

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

For the eigenvalue µ2 = 2:

2 * I - A ) * v   =   (2 * I - A ) * v   =
 -3 1 -3 1
*
 v1 v2
=
 - 3*v1 + v2   =   0 - 3*v1 + v2   =   0
 e1 e2

Equations e1 and e2 are the same. Solving for v2 in terms of v1 yields:

v2 = 3 * v1.

This means that the vector v2T = (1, 3), and any scalar multiple s * v2, are eigenvectors associated with the eigenvalue µ2 = 2.

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

The next diagram shows    v1,   A * v1,   v2,   A * v2,   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. A * u
u1: u2:

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 v1, and whose second column is the eigenvector v2. 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.

For example, the matrix:

A =
 3 1 -1 1

has the characteristic equation:

f(µ)   =   µ2 - trace(A) * µ + det(A) =   µ2 - 4 * µ + 4,

with eigenvalues µ1 = 2, and µ2 = 2. However, the only eigenvector is vT1 = (1, -1) and its multiples.

Recall that the eigenvector vT1 = (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
*
 v1 v2
=
 - 1*v1 - 1*v2   =   0 v1 +  v2   =   0
 e1 e2

Equations e1 and e2 are the same. Solving for v2 in terms of v1 yields:

v2 = -v1.

This means that the vector v1T = (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 v2 for the eigenvalues µ1,2 = 2 by solving the system of equations:

(A - 2 * I)*v2 = v1, or

(2 * I - A )*v2 = - v1.

Expanding:

(2 * I - A )* v   = -v1
 -1 -1 1 1
*
 -v1 -v2
=
 - 1*v1 - 1*v2   =   -1 v1 +  v2   =   1
 e1 e2

Adding equation e1 to equation e2, and multiplying e1 by -1, yields the echelon form:

 1*v1 + 1*v2   =   1 0 * v1 + 0 *v2   =   0
 e1 e2

The particular solution vTp = (1, 0) = v2 is the generalized eigenvalue of the matrix A associated with the pair of equal eigenvalues µ1,2 = 2,  while the homogeneous solution vh is a multiple of the eigenvector v1:

v=vp + s * vh
 v1 v2
=
 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 v1 satisfies:

(A - 2 * I)*v1 = o, or   A * v1 = 2 * v1.

As the generalized eigenvector constructed from v1, v2 satisfies:

(A - 2 * I)*v2 = v1, or   A * v2 = 2 * v2 + 1 * v1.

Jordan Decomposition.

Let V be the matrix whose first column is the eigenvector v1, and whose second column is the generaralized eigenvector v2. 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 v1 and v2, A * V = V * J. Because v1 and v2 are linearly independent, the matrix V has an inverse. The Jordan Decomposition is:

A   =   V * J * V(-1)   =
 1 1 -1 0
*
 2 1 0 2
*
 0 -1 1 1
=
 3 1 -1 1

Polynomials and their roots are described on the linear algebra: polynomials web page.

Systems of linear equations and their solutions are described on the linear algebra: linear equations web page.

Matrix Decompositions.

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 = PT * L * U.

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

Tableau0
AL
 R01 R02 R03 R04 R05
 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

Tableau1
UL
Switch R1 and R5
 R11 R12 = R02 - (0.7467) * R11 R13 = R03 - (0.7867) * R11 R14 = R04 - (0.8667) * R11 R15 = R05 - (0.7333) * R11
 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

Tableau2
UL
Switch R2 and R5
 R21 R22 R23 = R13 - (0.35) * R22 R24 = R14 - (0.125) * R22 R25 = R15 - (0.65) * R22
 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

Tableau3
UL
Switch R3 and R5
 R31 R32 R33 R34 = R24 - (0.1923) * R33 R35 = R25 - (0.5385) * R33
 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

Tableau4
UL
Switch R4 and R5
 R41 R42 R43 R44 R45 = R35 - (0.3571) * R44
 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

PA
 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

LU
 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.

det(A) = 1 * (75) * (-2.6667) * (-0.65) * (-0.5385) * (-0.3571) = 25

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.

Matrix A
 Square dimension of A: 5

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 =
a1a2a3a4
=
 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 { q1, q2, q3, q4}, from the set of 4 linearly independent vectors { a1, a2, a3, a4}.

Begin by writing q1 = a1. Consider a2. The projection of a2 on q1 is the vector:

p1 = ((q1 * a2) / (q1 * q1)) * q1.

The vector

q2 = a2 - p1

is orthogonal to q1, and {q1, q2} span the same vector subspace as {a1, a2}.

Repeat this process with the remaining { a3, a4} vectors. At stage k for a general matrix A, subtract the projections of ak on the orthogonal vectors {q1, ... qk-1} from ak to construct the next orthogonal vector qk.

Do this at each stage k, by constructing new projection vectors, pj, obtained by projecting ak on each qj as given by:

pj = ((qj * ak) / (qj * qj)) * qj,   for j = 1, . . . , k-1.

The required vector qk orthogonal to {q1, ... qk-1} is:

qk = ak  -  p1  -  p2  -  .  .  .  -  pk-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.

q1 =
 1 2 1 4

q2 =
 2 1 2 1
- (0.4545 ) *
 1 2 1 4
=
 1.54545 0.09091 1.54545 -0.81818

q3 =
 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

q4 =
 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
11.5455-0.93330.2206
20.0909-0.4667-0.5147
11.54551.0667-0.0735
4-0.81820.20.2206
Orthonormal Matrix Q
0.21320.6617-0.61990.3638
0.42640.0389-0.31-0.8489
0.21320.66170.7085-0.1213
0.8528-0.35030.13280.3638

Gram-Schmidt Decomposition A = Q * R

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

QT * A = QT * Q * R = R.

Expanding this last equation yields:

R =
 (q1T * a1) (q1T * a2) (q1T * a3) (q1T * a4) (q2T * a1) (q2T * a2) (q2T * a3) (q2T * a4) (q3T * a1) (q3T * a2) (q3T * a3) (q3T * a4) (q4T * a1) (q4T * a2) (q4T * a3) (q4T * a4)
=
 (q1T * a1) (q1T * a2) (q1T * a3) (q1T * a4) (q2T * a2) (q2T * a3) (q2T * a4) (q3T * a3) (q3T * a4) (q4T * a4)

The entries below the diagonal of R are zero, because aj is orthogonal to qk for j < k.

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

Matrix A
1213
2132
1232
4174
=
Orthonormal Matrix Q
0.21320.6617-0.61990.3638
0.42640.0389-0.31-0.8489
0.21320.66170.7085-0.1213
0.8528-0.35030.13280.3638
*
Matrix R
4.69042.1328.10165.33
02.33550.31141.9852
001.5055-0.5314
0000.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.

Matrix A
 Square dimension of A: 4

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 * uT

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|| uT  *  (2 * p) = 0, since they are perpendicular To prove that H * v = w, observe that

 H * v = (I - 2 * u * uT) * v = v - 2 * u * (uT * v), and 2 * u* (uT * v) = u * (uT * [(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 * uT 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, H1, H2, H3.

Starting with A1 = A, let v1 = the first column of A1, and w1T = (norm(v1), 0,...0), ie a column vector whose first component is the norm of v1 with the remaining components equal to 0. The Householder transformation H1 = I - 2 * u1 * u1T with u1 = v1 - w1 / ||v1 - w1|| will turn the first column of A1 into w1 as with H1 * A1 = A2. At each stage k, vk = the kth column of Ak on and below the diagonal with all other components equal to 0, and wk's kth component equals the norm of vk with all other components equal to 0. Letting Hk * Ak = Ak+1, the components of the kth column of Ak+1 below the diagonal are each 0. These calculations are listed below for each stage for the matrix A.

v1w1u1H1 = I - 2*u1*u1T*A1=A2
 16 5 9 4
 -19.4422 0 0 0
 0.9547 0.1347 0.2424 0.1077
 -0.823 -0.2572 -0.4629 -0.2057 -0.2572 0.9637 -0.0653 -0.029 -0.4629 -0.0653 0.8825 -0.0522 -0.2057 -0.029 -0.0522 0.9768
*
 16 2 3 5 11 10 9 7 6 4 14 15
=
 -19.4422 -10.5955 -10.9041 0 9.2231 8.0385 0 3.8016 2.4693 0 12.5785 13.4308

v2w2u2H2 = I - 2*u2*u2T*A2=A3
 0 9.2231 3.8016 12.5785
 0 -16.0541 0 0
 0 0.8873 0.1334 0.4415
 1 0 0 0 0 -0.5745 -0.2368 -0.7835 0 -0.2368 0.9644 -0.1178 0 -0.7835 -0.1178 0.6101
*
 -19.4422 -10.5955 -10.9041 0 9.2231 8.0385 0 3.8016 2.4693 0 12.5785 13.4308
=
 -19.4422 -10.5955 -10.9041 0 -16.0541 -15.7259 0 0 -1.1048 0 0 1.6051

v3w3u3H3 = I - 2*u3*u3T*A3=A4
 0 0 -1.1048 1.6051
 0 0 1.9486 0
 0 0 -0.8851 0.4653
 1 0 0 0 0 1 0 0 0 0 -0.567 0.8237 0 0 0.8237 0.567
*
 -19.4422 -10.5955 -10.9041 0 -16.0541 -15.7259 0 0 -1.1048 0 0 1.6051
=
 -19.4422 -10.5955 -10.9041 0 -16.0541 -15.7259 0 0 1.9486 0 0 0

Writing QT = H3 * H2 * H1and R = A4, then:

A = Q * R

A=Q*R
 16 2 3 5 11 10 9 7 6 4 14 15
=
 -0.823 0.4186 0.3123 -0.2236 -0.2572 -0.5155 -0.4671 -0.6708 -0.4629 -0.1305 -0.5645 0.6708 -0.2057 -0.7363 0.6046 0.2236
*
 -19.4422 -10.5955 -10.9041 0 -16.0541 -15.7259 0 0 1.9486 0 0 0

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.

Matrix A
 Row dimension of A: m:   4 Column dimension of A: n:   3

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 vT = (v1, v2) into a vector wT = (w1, 0).

JT * v =
 cos(Ø) -sin(Ø) sin(Ø) cos(Ø)
 v1 v2
=
 w1 0

To get w2 = 0 requires;

cos(Ø) = v1 / ||v||,

sin(Ø) = -v2 / ||v||, and

w1 = norm(v) = ||v|| = sqrt(v1*v1 + v2*v2).

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.

v1w1J1T*A1=A2
 16 5 0 0
 16.7631 0 0 0
 0.9545 0.2983 0 0 -0.2983 0.9545 0 0 0 0 1 0 0 0 0 1
*
 16 2 3 5 11 10 9 7 6 4 14 15
=
 16.7631 5.19 5.8462 0 9.9027 8.65 9 7 6 4 14 15

v2w2J2T*A2=A3
 16.7631 0 9 0
 19.0263 0 0 0
 0.881 0 0.473 0 0 1 0 0 -0.473 0 0.881 0 0 0 0 1
*
 16.7631 5.19 5.8462 0 9.9027 8.65 9 7 6 4 14 15
=
 19.0263 7.8838 7.9889 0 9.9027 8.65 0 3.7123 2.5209 4 14 15

v3w3J3T*A3=A4
 19.0263 0 0 4
 19.4422 0 0 0
 0.9786 0 0 0.2057 0 1 0 0 0 0 1 0 -0.2057 0 0 0.9786
*
 19.0263 7.8838 7.9889 0 9.9027 8.65 0 3.7123 2.5209 4 14 15
=
 19.4422 10.5955 10.9041 0 9.9027 8.65 0 3.7123 2.5209 0 12.0785 13.0355

v4w4J4T*A4=A5
 0 9.9027 3.7123 0
 0 10.5757 0 0
 1 0 0 0 0 0.9364 0.351 0 0 -0.351 0.9364 0 0 0 0 1
*
 19.4422 10.5955 10.9041 0 9.9027 8.65 0 3.7123 2.5209 0 12.0785 13.0355
=
 19.4422 10.5955 10.9041 0 10.5757 8.9844 0 0 -0.6759 0 12.0785 13.0355

v5w5J5T*A5=A6
 0 10.5757 0 12.0785
 0 16.0541 0 0
 1 0 0 0 0 0.6588 0 0.7524 0 0 1 0 0 -0.7524 0 0.6588
*
 19.4422 10.5955 10.9041 0 10.5757 8.9844 0 0 -0.6759 0 12.0785 13.0355
=
 19.4422 10.5955 10.9041 0 16.0541 15.7259 0 0 -0.6759 0 0 1.8276

v6w6J6T*A6=A7
 0 0 -0.6759 1.8276
 0 0 1.9486 0
 1 0 0 0 0 1 0 0 0 0 -0.3469 0.9379 0 0 -0.9379 -0.3469
*
 19.4422 10.5955 10.9041 0 16.0541 15.7259 0 0 -0.6759 0 0 1.8276
=
 19.4422 10.5955 10.9041 0 16.0541 15.7259 0 0 1.9486 0 0 0

Writing QT = J6T * J5T * J4T * J3T * J2T * J1Tand R = A7, then:

Q * QT * A = A = Q * R

A=Q*R
 16 2 3 5 11 10 9 7 6 4 14 15
=
 0.823 -0.4186 0.3123 0.2236 0.2572 0.5155 -0.4671 0.6708 0.4629 0.1305 -0.5645 -0.6708 0.2057 0.7363 0.6046 -0.2236
*
 19.4422 10.5955 10.9041 0 16.0541 15.7259 0 0 1.9486 0 0 0

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.

Matrix A
 Row dimension of A: m:   4 Column dimension of A: n:   3

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* VT,   where - the m by m orthogonal matrix U consists of the eigenvectors of A * AT. - the n by n orthogonal matrix V consists of the eigenvectors of AT * A. - the m by n diagonal matrix S consists of the square roots of the eigenvalues of A * AT (which equal the eigenvalues of AT * A), arranged in descending order. - the eigenvectors of A * AT and AT * A are arranged in the columns of U and V, respectively, in order of their shared eigenvalues1/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+ * UT - 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 * VT

A = U*S*VT
 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+ * UT

A+ = V*S+*UT
 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.

Matrix A
 Row dimension of A: m:   3 Column dimension of A: n:   4

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. Copyright © Elmer G. Wiens:   Egwald Web Services All Rights Reserved.    Inquiries 