Matrix Computations

This tutorial reviews the functions that Mathematica provides for carrying out matrix computations. Further information on these functions can be found in standard mathematical texts by such authors as Golub and van Loan or Meyer. The operations described in this tutorial are unique to matrices; an exception is the computation of norms, which also extends to scalars and vectors.

Basic Operations

This section gives a review of some basic concepts and operations that will be used throughout the tutorial to discuss matrix operations.

Norms

The norm of a mathematical object is a measurement of the length, size, or extent of the object. In Mathematica norms are available for scalars, vectors, and matrices.

Computing norms in Mathematica.

For numbers, the norm is the absolute value:

Norm[-4]

Norm[ 4.  +6. I]

Vector Norms

For vector spaces, norms allow a measure of distance. This allows the definition of many familiar concepts such as neighborhood, closeness of approximation, and goodness of fit. A vector norm is a function f:(ℝ^n)⟶ ℝ that satisfies the following relations.

({{f(x)>0, x∈(ℝ^n),           (f(x)==0 iff x==0)}, {f(x+y)≤f(x)+f(y), x,y∈(ℝ^n),}, {f(α x)==∣α∣f(x), α∈ℝ,x∈(ℝ^n)}})

Typically this function uses the notation ∥x(Subscript[∥, p]). The subscript is used to distinguish different norms, of which the p-norms are particularly important. For x=={(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n])}, the p-norm is defined as follows.

∥x(Subscript[∥, p])==(((∣(Subscript[x, 1])(∣^p)+…+∣(Subscript[x, n])(∣^p)))^(1/p))          p≥1

Some common p-norms are the 1-, 2-, and ∞-norms.

({{∥x(Subscript[∥, 1])==∣(Subscript[x, 1])∣+…+∣(Subscript[x, n])∣}, {∥x(Subscript[∥, 2])==(((∣(Subscript[x, 1])(∣^2)+…+∣(Subscript[x, n])(∣^2)))^(1/2))}, {∥x(Subscript[∥, ∞])==(max_(1≤i≤n))∣(Subscript[x, i])∣}})

In Mathematica, vector p-norms can be computed with the function (RefLink[Norm, paclet:ref/Norm]). The 1-, 2-, and ∞-norms are demonstrated in the following examples:

vec={1., 2., 5.6};

Norm[vec,1]

Norm[vec,2]

Norm[vec,∞]

The 2-norm is particularly useful and this is the default:

Norm[vec]

Norms are implemented for vectors with exact numerical entries:

Norm[ {1,2,3}]

Norms are also implemented for vectors with symbolic entries:

Norm[ {x,y,z}]

In addition, if a symbolic p is used, the result is symbolic:

Norm[vec,p]

Matrix Norms

Matrix norms are used to give a measure of distance in matrix spaces. This is necessary if you want to quantify what it means for one matrix to be near to another, for example, to say that a matrix is nearly singular.

Matrix norms also use the double bar notation of vector norms. One of the most common matrix norms is the Frobenius norm (also called the Euclidean norm).

∥A(Subscript[∥, E])==(Sqrt[(∑_(i==1)^m)(∑_(j==1)^n)∣(Subscript[a, ij])(∣^2)])

Other common norms are the p-norms. These are defined in terms of vector p-norms as follows.

∥A(Subscript[∥, p])==(sup_(x≠0))((∥A x(Subscript[∥, p]))/(∥x(Subscript[∥, p])))

Thus the matrix p-norms show the maximum expansion that a matrix can apply to any vector.

In Mathematica the Frobenius norm can be computed as follows:

mat={{1.,2.},{4.5,2.1}};

Norm[ mat,"Frobenius"]

The matrix 2-norm can be computed as follows:

Norm[ {{1.,2.},{4.5,2.1}}]

It is also possible to give an argument to specify the 1, 2, or ∞ matrix p-norms:

Norm[ {{1.,2.},{4.5,2.1}},1]

Norm[ {{1.,2.},{4.5,2.1}},2]

Norm[ {{1.,2.},{4.5,2.1}},∞]

This computes the expansion that a collection of different vectors gets from an input matrix. You can see that the maximum value is still less than the 2-norm:

Max[Table[ vec={Random[],Random[]};Norm[mat.vec]/Norm[vec],{100}]]

All the p-norms are supported for exact numerical matrices:

Norm[ {{1,2},{4,3}},1]

Norm[  {{1,2},{4,3}},2]

Norm[  {{1,2},{4,3}},∞]

However, only the 1 and ∞ matrix p-norms are supported for symbolic matrices:

Norm[ {{a,b},{c,d}},1]

Norm[  {{a,b},{c,d}},∞]

NullSpace

One of the fundamental subspaces associated with each matrix is the null space. Vectors in the null space of a matrix are mapped to zero by the action of the matrix.

For some matrices (for instance, nonsingular square matrices), the null space is empty:

mat={{1., 2.},{3., 4.}};
NullSpace[mat]

This matrix has a null space with one vector:

mat={{1., 2.},{2., 4.}};
null=NullSpace[mat]

mat.null[[1]]

This matrix has a null space with two basis vectors:

mat = {{1.,1.,1.}, {2.,2.,2.}, {3.,3.,3.}};
NullSpace[mat]

Rank

The rank of a matrix corresponds to the number of linearly independent rows or columns in the matrix.

If an m×n matrix has no linear dependencies between its rows, its rank is equal to min(m,n) and it is said to be full rank:

mat={{1., 2.},{3., 4.}};
MatrixRank[mat]

If a matrix has any dependencies in its rows then its rank is less than min(m,n). In this case the matrix is said to be rank deficient:

mat={{1., 2.},{2., 4.}};
MatrixRank[mat]

Note that the rank of a matrix is equal to the rank of its transpose. This means that the number of linearly independent rows is equal to the number of linearly independent columns:

MatrixRank[ Transpose[mat]]

For an m×n matrix A the following relations hold: (RefLink[Length, paclet:ref/Length])[(RefLink[NullSpace, paclet:ref/NullSpace])[A]]+ (RefLink[MatrixRank, paclet:ref/MatrixRank])[A]==n, and (RefLink[Length, paclet:ref/Length])[(RefLink[NullSpace, paclet:ref/NullSpace])[(A^T)]]+ (RefLink[MatrixRank, paclet:ref/MatrixRank])[(A^T)]==m. From this it follows that the null space is empty if and only if the rank is equal to n and that the null space of the transpose of A is empty if and only if the rank of A is equal to m.

If the rank is equal to the number of columns, it is said to have full column rank. If the rank is equal to the number of rows, it is said to have full row rank. One way to understand the rank of a matrix is to consider the row echelon form.

Reduced Row Echelon Form

A matrix can be reduced to a row echelon form by a combination of row operations that start with a pivot position at the top-left element and subtract multiples of the pivot row from following rows so that all entries in the column below the pivot are zero. The next pivot is chosen by going to the next row and column. If this pivot is zero and any nonzero entries are in the column beneath, the rows are exchanged and the process is repeated. The process finishes when the last row or column is reached.

A matrix is in row echelon form if any row that consists entirely of zeros is followed only by other zero rows and if the first nonzero entry in row i is in the column j, then elements in columns from 1 to j in all rows below i are zero. The row echelon form is not unique but its form (in the sense of the positions of the pivots) is unique. It also gives a way to determine the rank of a matrix as the number of nonzero rows.

A matrix is in reduced row echelon form if it is in row echelon form, each pivot is one, and all the entries above (and below) each pivot are zero. It can be formed by a procedure similar to the procedure for the row echelon form, also taking steps to reduce the pivot to one (by division) and reduce entries in the column above each pivot to zero (by subtracting multiples of the current pivot row). The reduced row echelon form of a matrix is unique.

The reduced row echelon form (and row echelon form) give a way to determine the rank of a matrix as the number of nonzero rows. In Mathematica the reduced row echelon form of a matrix can be computed by the function (RefLink[RowReduce, paclet:ref/RowReduce]).

The reduced row echelon form of this matrix only has one nonzero row. This means that the rank is 1:

mat={{1., 2.},{2., 4.}};
RowReduce[mat]//MatrixForm

This is a 3×2 random matrix whose columns are linearly independent:

mat=Table[Random[],{3},{2}];
mat//MatrixForm

The reduced row echelon form has two nonzero rows; thus, its rank should be 2:

RowReduce[mat]//MatrixForm

(RefLink[MatrixRank, paclet:ref/MatrixRank]) also computes the rank as 2:

MatrixRank[mat]

The reduced row echelon form of the transpose of the matrix also has two nonzero rows. This is consistent with the fact that the rank of a matrix is equal to the rank of its transpose, even when the matrix is rectangular:

RowReduce[Transpose[mat]]//MatrixForm

Inverse

The inverse of a square matrix A is defined by

A.(A^(-1))==I

where I is the identity matrix. The inverse can be computed in Mathematica with the function (RefLink[Inverse, paclet:ref/Inverse]).

Here is a sample 2×2 matrix:

mat={{1., 2.},{3.,5.}};
invMat=Inverse[mat];
invMat//MatrixForm

This demonstrates that the matrix satisfies the definition:

mat.invMat//Chop//MatrixForm

Not all matrices have an inverse; if a matrix does not have an inverse it is said to be singular. If a matrix is rank-deficient then it is singular:

mat={{1., 2.},{2., 4.}};
MatrixRank[mat]

Inverse[ mat]

The matrix inverse can in principle be used to solve the matrix equation

A.X==B

by multiplying it by the inverse of A.

X==(A^(-1)).B

However, it is always better to solve the matrix equation directly. Such techniques are discussed in the section "Solving Linear Systems".

PseudoInverse

When the matrix is singular or is not square it is still possible to find an approximate inverse that minimizes ∥A.(A^(-1))-I(Subscript[∥, p]). When the 2-norm is used this will find a least squares solution known as the pseudoinverse.

This finds the pseudoinverse for the singular matrix defined earlier:

mat={{1., 2.},{2., 4.}};
matInv=PseudoInverse[mat];
matInv//MatrixForm

The result is not the identity matrix, but it is "close" to the identity matrix:

mat.matInv//MatrixForm

This computes the pseudoinverse of a rectangular matrix:

mat={{1., 2.,3.},{2.,5, 4.}};
matInv=PseudoInverse[mat];
matInv//MatrixForm

The result is quite close to the identity matrix:

mat.matInv//MatrixForm

The solution found by the pseudoinverse is a least squares solution. These are discussed in more detail under "Least Squares Solutions".

Determinant

The determinant of an n×n matrix is defined as follows.

({{det(A)==(Subscript[a, 11])      if n==1}, {det(A)==(∑_(i==1)^n)(((-1))^(1+i))(Subscript[a, i1])det((Subscript[A, i1]))}})

It can be computed in Mathematica with the function (RefLink[Det, paclet:ref/Det]).

Here is a sample 2×2 matrix:

mat={{1.,2.},{3.,5.}};
Det[mat]

One useful property of the determinant is that (det(A) == 0) if and only if A is singular:

mat={{1.,2.},{2.,4.}};
Det[mat]

As already pointed out, if the matrix is singular, it does not have full rank:

MatrixRank[mat]

Minors

A minor of a matrix is the determinant of any k×k submatrix. This example uses the function (RefLink[Minors, paclet:ref/Minors]) to compute all the 2×2 minors:

mat={{1,2,3},{4,5,6},{7,8,9}};
Minors[mat]

The size of the submatrices can be controlled by a second option. Here the 1×1 minors are computed; these are just the entries of the matrix:

Minors[ mat,1]

Note that (RefLink[Minors, paclet:ref/Minors])[m] is equivalent to (RefLink[Minors, paclet:ref/Minors])[m,n-1]. In general, (RefLink[Minors, paclet:ref/Minors])[m,k] computes determinants of all the possible k×k submatrices that can be generated from the matrix m (this is done by picking different k rows and k columns). For the 2×2 minors of a 4×4 matrix you obtain a 6×6 matrix because there are only 6 different ways to pick 2 rows from 4 rows, and to pick 2 columns from 4 columns:

mat = N[{{1,4,5,-1}, {2,3,0,4}, {3,2,-1,-2}, {4,3,2,1}}];Minors[mat,2]

One of the properties of the rank of a matrix is that it is equal to the size of the largest nonsingular submatrix. This can be demonstrated with (RefLink[Minors, paclet:ref/Minors]). In this example, the following matrix has a rank of 2:

mat={{1,2,3},{4,5,6},{3,6,9}};
MatrixRank[mat]

The determinant of the matrix is zero:

Det[mat]

When the 2×2 minors are computed you can see that not all are zero, confirming that the rank of the matrix is 2:

Minors[mat]

Solving Linear Systems

One of the most important uses of matrices is to represent and solve linear systems. This section discusses how to solve linear systems with Mathematica. It makes strong use of (RefLink[LinearSolve, paclet:ref/LinearSolve]), the main function provided for this purpose.

Solving a linear system involves solving a matrix equation A.X==B. Because A is an m×n matrix this is a set of m linear equations in n unknowns.

(({{(Subscript[a, 11]), (Subscript[a, 12]), …, (Subscript[a, 1n])}, {(Subscript[a, 21]), (Subscript[a, 22]), …, (Subscript[a, 2n])}, {⋮, ⋮, ⋱, ⋮}, {(Subscript[a, m1]), (Subscript[a, m2]), …, (Subscript[a, mn])}}))(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {⋮}, {(Subscript[x, n])}}))==(({{(Subscript[b, 1])}, {(Subscript[b, 2])}, {⋮}, {(Subscript[b, m])}}))

When m==n the system is said to be square. If m>n there are more equations than unknowns and the system is overdetermined. If m<n there are fewer equations than unknowns and the system is underdetermined.

Classes of linear systems represented by rectangular matrices.

Note that even though you could solve the matrix equation by computing the inverse via X==(A^(-1)).B, this is not a recommended method. You should use a function that is designed to solve linear systems directly and in Mathematica this is provided by (RefLink[LinearSolve, paclet:ref/LinearSolve]).

Solving linear systems with (RefLink[LinearSolve, paclet:ref/LinearSolve]).

This uses (RefLink[LinearSolve, paclet:ref/LinearSolve]) to solve a matrix equation:

A={{1., 2.},{4.,6.}};
B={5.4,2.};
X=LinearSolve[A,B]

The solution can be tested to demonstrate that it solves the problem:

A.X

You can always generate the linear equations from the input matrices:

Thread[ A.{x1,x2}==B]

(RefLink[Solve, paclet:ref/Solve]) can then be used to find the solution:

Solve[%]

(RefLink[LinearSolve, paclet:ref/LinearSolve]) also works when the right-hand side of the matrix equation is also a matrix:

A={{1., 2.},{4.,6.}};
B={{5.4,2.},{5.6, 1.1}};
X=LinearSolve[A,B]

A.X

If the system is overdetermined, (RefLink[LinearSolve, paclet:ref/LinearSolve]) will find a solution if it exists:

A={{1.,2.},{5.,6.},{4.5,6.}};
B={5.,6.,8.25};

LinearSolve[A,B]

If no solution can be found, the system is inconsistent. In this case it might still be useful to find a best-fit solution. This is often done with a least squares technique, which is explored under "Least Squares Solutions".

If the system is underdetermined, there is either no solution or infinitely many solutions and (RefLink[LinearSolve, paclet:ref/LinearSolve]) will return one:

A={{1.,2.,3.},{5.,6.,7.}};
B={5.,6.};

LinearSolve[A,B]

A system is consistent if and only if the rank of the matrix A equals the rank of the matrix formed by adding B as an extra column to A (known as the augmented matrix). The following demonstrates that the previous system is consistent:

MatrixRank[A]==MatrixRank[ {{1.,2.,3.,5.},{5.,6.,7.,6.}}]

(RefLink[LinearSolve, paclet:ref/LinearSolve]) can work with all the different types of matrices that Mathematica supports, including symbolic and exact, machine number, and arbitrary precision. It also works with both dense and sparse matrices. When a sparse matrix is used in (RefLink[LinearSolve, paclet:ref/LinearSolve]), special techniques suitable for sparse matrices are used and the result is the solution vector:

A=SparseArray[{i_,j_}/;i≥j:>Random[],{3,3}];
B={1., 2., 1.1};
X=LinearSolve[A,B]

Options for (RefLink[LinearSolve, paclet:ref/LinearSolve]).

(RefLink[LinearSolve, paclet:ref/LinearSolve]) has three options that allow users more control. The (RefLink[Modulus, paclet:ref/Modulus]) and (RefLink[ZeroTest, paclet:ref/ZeroTest]) options are used by symbolic techniques and are discussed in the section "Symbolic and Exact Matrices". The (RefLink[Method, paclet:ref/Method]) option allows you to choose methods that are known to be suitable for a particular problem; the default setting of (RefLink[Automatic, paclet:ref/Automatic]) means that (RefLink[LinearSolve, paclet:ref/LinearSolve]) chooses a method based on the input.

Singular Matrices

A matrix is singular if its inverse does not exist. One way to test for singularity is to compute a determinant; this will be zero if the matrix is singular. For example, the following matrix is singular:

A={{1.,2.},{1.,2.}};
Det[A]

For many right-hand sides there is no solution:

B={4., 1.};
LinearSolve[A,B]

However, for certain values there will be a solution:

B={4., 4.};
LinearSolve[A,B]

For the first example, the rank of A does not equal that of the augmented matrix. This confirms that the system is not consistent and cannot be solved by (RefLink[LinearSolve, paclet:ref/LinearSolve]):

MatrixRank[A]==MatrixRank[{{1.,2.,4.},{1.,2.,1.}}]

For the second example, the rank of A does equal that of the augmented matrix. This confirms that the system is consistent and can be solved by (RefLink[LinearSolve, paclet:ref/LinearSolve]):

MatrixRank[A]==MatrixRank[{{1.,2.,4.},{1.,2.,4.}}]

In those cases where a solution cannot be found, it is still possible to find a solution that makes a best fit to the problem. One important class of best-fit solutions involves least squares solutions and is discussed in "Least Squares Solutions".

Homogeneous Equations

The homogeneous matrix equation involves a zero right-hand side.

A.x==0

This equation has a nonzero solution if the matrix is singular. To test if a matrix is singular, you can compute the determinant.

A={{1.,2.},{1.,2.}};
Det[A]

To find the solution of the homogeneous equation you can use the function (RefLink[NullSpace, paclet:ref/NullSpace]). This returns a set of orthogonal vectors, each of which solves the homogeneous equation. In the next example, there is only one vector:

null=NullSpace[A]

This demonstrates that the solution in fact solves the homogeneous equation:

A.null[[1]]

The function (RefLink[Chop, paclet:ref/Chop]) can be used to replace approximate numbers close to 0:

%//Chop

The solution to the homogeneous equation can be used to form an infinite number of solutions to the inhomogeneous equation. This solves an inhomogeneous equation:

B={4.,4.};
sol=LinearSolve[A,B]

The solution does in fact solve the equation:

A.sol

If you add to sol an arbitrary factor times the homogeneous solution, this new vector also solves the matrix equation:

A.(sol+Random[]  null[[1]])

Estimating and Calculating Accuracy

An important way to quantify the accuracy of the solution of a linear system is to compute the condition number κ(A). This is defined as follows for a suitable choice of norm.

κ(A)==∥A∥∥(A^(-1))∥

It can be shown that for the matrix equation Ax==B the relative error in x is κ(A) times the relative error in A and B. Thus, the condition number quantifies the accuracy of the solution. If the condition number of a matrix is large, the matrix is said to be ill-conditioned. You cannot expect a good solution from an ill-conditioned system. For some extremely ill-conditioned systems it is not possible to obtain any solution.

In Mathematica an approximation of the condition number can be computed with the function LinearAlgebra`MatrixConditionNumber.

This computes the condition number of a matrix:

mat={{1.,2.},{3., 4.}};
LinearAlgebra`MatrixConditionNumber[mat]

This matrix is singular and the condition number is ∞:

mat={{1.,2.},{1., 2.}};
LinearAlgebra`MatrixConditionNumber[mat]

This matrix has a large condition number and is said to be ill-conditioned:

mat=N[{{1,2,3},{4,5,6},
{3,6,9+(10^11) $MachineEpsilon}}];
LinearAlgebra`MatrixConditionNumber[mat]

If a matrix is multiplied with itself, the condition number increases:

mat2=mat.mat;
LinearAlgebra`MatrixConditionNumber[mat2]

When you solve a matrix equation involving an ill-conditioned matrix, the result may not be as accurate:

vec={3,4,5};
sol=LinearSolve[mat,vec];
Norm[mat.sol-vec]

The solution is less accurate for the matrix mat2:

sol=LinearSolve[mat2,vec];
Norm[mat2.sol-vec]

The solution to these problems is to avoid constructing matrices that may be ill-conditioned, for example, by avoiding multiplying matrices by themselves. Another method of solving problems more accurately is to use Mathematica's arbitrary-precision computation; this is discussed under "Arbitrary-Precision Matrices":

matAcc = SetPrecision[mat2,16];
sol=LinearSolve[matAcc,vec];
Norm[matAcc.sol-vec]

Symbolic and Exact Matrices

(RefLink[LinearSolve, paclet:ref/LinearSolve]) works for all the different types of matrices that can be represented in Mathematica. These are described in more detail in "Matrix Types".

Here is an example of working with purely symbolic matrices:

A={{x1,x2},{x3,x4}};
MatrixForm[A]

b={y1,y2};
X=LinearSolve[A,b]

The original matrix can be multiplied into the solution:

A.X

The result of multiplication needs some postprocessing with (RefLink[Simplify, paclet:ref/Simplify]) to reach the expected value. This demonstrates the need for intermediate processing in symbolic linear algebra computation. Without care these intermediate results can become extremely large and make the computation run slowly:

%//Simplify

In order to simplify the intermediate expressions, the (RefLink[ZeroTest, paclet:ref/ZeroTest]) option might be useful.

(RefLink[LinearSolve, paclet:ref/LinearSolve]) can also be used to solve exact problems:

A={{1,2},{3,4}};
MatrixForm[A]

b={6,8};
X=LinearSolve[A,b]

A.X

There are a number of methods that are specific to symbolic and exact computation: CofactorExpansion, DivisionFreeRowReduction, and OneStepRowReduction. These are discussed in the section "Symbolic Methods".

Row Reduction

Row reduction involves adding multiples of rows together so as to produce zero elements when possible. The final matrix is in reduced row echelon form as described previously. If the input matrix is a nonsingular square matrix, the result will be an identity matrix:

RowReduce[ {{1,2},{1,1}}]

This can be used as a way to solve a system of equations:

RowReduce[ {{1,2,x},{1,1,y}}]

This can be compared with the following use of (RefLink[LinearSolve, paclet:ref/LinearSolve]):

LinearSolve[ {{1,2},{1,1}},{x,y}]

Saving the Factorization

Many applications of linear systems involve the same matrix A but different right-hand sides B. Because a significant part of the effort involved in solving the system involves processing A, it is common to save the factorization and use it to solve repeated problems. In Mathematica, this is done by using a one-argument form of (RefLink[LinearSolve, paclet:ref/LinearSolve]); this returns a functional that you can apply to different vectors to obtain each solution:

A={{1., 2.},{4.,6.}};
fun=LinearSolve[A]

When you can apply the (RefLink[LinearSolveFunction, paclet:ref/LinearSolveFunction]) to a particular right-hand side the solution is obtained:

B={1., 20.};
X=fun[B]

This solves the matrix equation:

A.X

A different right-hand side yields a different solution:

B={11., 2.};
X=fun[B]

This new solution solves this matrix equation:

A.X

The one-argument form of (RefLink[LinearSolve, paclet:ref/LinearSolve]) works in a completely equivalent way to the two-argument form. It works with the same range of input matrices, for example, returning the expected results for symbolic, exact, or sparse matrix input. It also accepts the same options.

One issue with the one-argument form is that for certain input matrices the factorization cannot be saved. For example, if the system is overdetermined, the exact solution is not certain to exist. In this case a message is generated that warns you that the factorization will be repeated each time the functional is applied to a particular right-hand side:

A={{1.,2.},{5.,6.},{4.5,6.}};
fun=LinearSolve[A]

For this right-hand side there is a solution, and this is returned:

B={5.,6.,8.25};
fun[B]

However, for this right-hand side there is no solution, and an error is encountered:

B={5.,6.,1.};
fun[B]

Methods

(RefLink[LinearSolve, paclet:ref/LinearSolve]) provides a number of different techniques to solve matrix equations that are specialized to particular problems. You can select between these using the option (RefLink[Method, paclet:ref/Method]). In this way a uniform interface is provided to all the functionality that Mathematica provides for solving matrix equations.

The default setting of the option (RefLink[Method, paclet:ref/Method]) is (RefLink[Automatic, paclet:ref/Automatic]). As is typical for Mathematica, this means that the system will make an automatic choice of method to use. For (RefLink[LinearSolve, paclet:ref/LinearSolve]) if the input matrix is numerical and dense, then a method using LAPACK routines is used; if it is numerical and sparse, a multifrontal direct solver is used. If the matrix is symbolic then specialized symbolic routines are used.

The details of the different methods are now described.

LAPACK

LAPACK is the default method for solving dense numerical matrices. When the matrix is square and nonsingular the routines dgesv, dlange, and dgecon are used for real matrices and zgesv, zlange, and zgecon for complex matrices. When the matrix is nonsquare or singular dgelss is used for real matrices and zgelss for complex matrices. More information on LAPACK is available in the references.

If the input matrix uses arbitrary-precision numbers, then LAPACK algorithms extended for arbitrary-precision computation are used.

Multifrontal

The multifrontal method is a direct solver used by default if the input matrix is sparse; it can also be selected by specifying a method option.

This reads in a sample sparse matrix stored in the Matrix Market format. More information on the importing of matrices is under "Import and Export of Sparse Matrices":

sp=Import["LinearAlgebraExamples/Data/dwg961b.mtx","MatrixMarket"]

vec=Table[Random[],{961}];

This solves the matrix equation using the sample matrix:

LinearSolve[ sp,vec,Method->Multifrontal];

If the input matrix to the multifrontal method is dense, it is converted to a sparse matrix.

The implementation of the multifrontal method uses the "UMFPACK" library.

Krylov

The Krylov method is an iterative solver that is suitable for large sparse linear systems, such as those arising from numerical solving of PDEs. In Mathematica two Krylov methods are implemented: conjugate gradient (for symmetric positive definite matrices) and BiCGSTAB (for nonsymmetric systems). It is possible to set the method that is used and a number of other parameters by using appropriate suboptions.

Suboptions of the Krylov method.

The default method for Krylov, BiCGSTAB, is more expensive but more generally applicable. The conjugate gradient method is suitable for symmetric positive definite systems, always converging to a solution (though the convergence may be slow). If the matrix is not symmetric positive definite, the conjugate gradient may not converge to a solution.

In this example, the matrix is not symmetric positive definite and the conjugate gradient method does not converge:

mat={{1,100.},{-3.,.1}};
vec={1.2,3};
LinearSolve[mat,vec,Method->{Krylov,Method->ConjugateGradient}]

The default method, BiCGSTAB, converges and returns the solution:

sol=LinearSolve[mat,vec,Method->Krylov]

This tests the solution that was found:

Norm[mat.sol-vec]

Here is an example that shows the benefit of the Krylov method and the use of a preconditioner. First, a function that builds a structured banded sparse matrix is defined:

BandedMatrix[n_]:= Module[
{A,spar,diag, udiag, udiag2, ldiag, ldiag2},
diag=Table[{i,i}->14.,{i,1,n}];
udiag=Table[{i,i+1}->3.,{i,1,n-1}];
udiag2=Table[{i,i+8}->-1.,{i,1,n-8}];
ldiag=Table[{i+1,i}->5.,{i,1,n-1}];
ldiag2=Table[{i+500,i}->-4.,{i,1,n-500}];
spar=Flatten[{diag,udiag,udiag2,ldiag,ldiag2}];
A=SparseArray[spar,{n,n}]
];

This demonstrates the structure of the matrices generated by this function:

MatrixPlot[BandedMatrix[700]]

This builds a larger system, with a 10000×10000 matrix:

num = 10000;
mat=BandedMatrix[num];
vec = mat.Table[i,{i,num}];

This will use the default sparse solver, a multifrontal direct solver:

sol=LinearSolve[mat,vec];//Timing

The Krylov method is faster:

LinearSolve[mat,vec,Method->Krylov];//Timing

The use of the ILU preconditioner is even faster. Currently, the preconditioner can only work if the matrix has nonzero diagonal elements:

LinearSolve[mat,vec,Method->{Krylov,Preconditioner->ILU0}];//Timing

At present, only the ILU preconditioner is built-in. You can still define your own preconditioner by defining a function that is applied to the input matrix. An example that involves solving a diagonal matrix is shown next.

First the problem is set up and solved without a preconditioner:

num=10000;
mat=SparseArray[{{i_,i_}->N[i]},{num,num}];
rhs=mat.Table[1.,{i,1,num}];
LinearSolve[mat,rhs,Method->Krylov];//Timing

Now, a preconditioner function fun is used; there is a significant improvement to the speed of the computation:

fun[x_]:=x/Table[N[i],{i,num}];
LinearSolve[mat,rhs,Method->{Krylov,Preconditioner->fun}];//Timing

Generally, a problem will not be structured so that it can benefit so much from such a simple preconditioner. However, this example is useful since it shows how to create and use your own preconditioner.

If the input matrix to the Krylov method is dense, the result is still found because the method is based on matrix/vector multiplication.

The Krylov method can be used to solve systems that are too large for a direct solver. However, it is not a general solver, being particularly suitable for those problems that have some form of diagonal dominance.

Cholesky

The Cholesky method is suitable for solving symmetric positive definite systems:

sp=Import["LinearAlgebraExamples/Data/nos6.mtx","MatrixMarket"]

B=Table[Random[],{675}];

LinearSolve[sp,B];//Timing

The Cholesky method is faster for this matrix:

LinearSolve[sp, B,Method->Cholesky];//Timing

The Cholesky method is also more stable. This can be demonstrated with the following matrix:

SeedRandom[4];
mat=Table[Random[],{1000},{1000}];
mat=mat.Transpose[mat];

For this matrix the default LU decomposition method reports a potential error:

LinearSolve[mat];//Timing

The Cholesky method succeeds:

LinearSolve[mat,Method->Cholesky];//Timing

For dense matrices the Cholesky method uses LAPACK functions such as dpotrf and dpotrs for real matrices and zpotrf and zpotrs for complex matrices. For sparse matrices the Cholesky method uses the "TAUCS" library.

Symbolic Methods

There are a number of methods that are specific to symbolic and exact computation: CofactorExpansion, DivisionFreeRowReduction, and OneStepRowReduction. These can be demonstrated for the following symbolic matrix:

mat={{a13*a2,a5,6^(1/2)*a1,2*a14},{18*a12,-3*a16,3*a9,a1+a3+a4},
{7*a15,15/14-(12*a11)/23+(9*a2)/2+(24*a5)/7+a9/5,
-5/23+(17*a152)/29+(9*a3)/4+(30*a15*a5)/23+(3*a52)/29-
(5*a6)/24,(4*a8)/3},{-9/16-(2*a12)/17-(7*a122)/5+(9*a163)/13-(7*a13*a72)/11+(10*a73)/7,
a11+a4,(11*a11)/9,-9/26+(13*a14)/29+7*a142+a15/4+25*a153-
(4*a2)/3-a22/6+(4*a15*a22)/3+(29*a4)/23+(25*a22*a4)/21}};
vec={1,0,3,0};

The computation is repeated for all three methods:

ByteCount[LinearSolve[mat,vec, Method -> CofactorExpansion]]//Timing

ByteCount[LinearSolve[mat,vec, Method -> DivisionFreeRowReduction]]//Timing

ByteCount[LinearSolve[mat,vec, Method -> OneStepRowReduction]]//Timing

Least Squares Solutions

This section follows from the previous section. It is concerned with finding solutions to the matrix equation A.X==B, where A is an m×n matrix and m>n (i.e., when there are more equations than unknowns). In this case, there may well be no solution, and (RefLink[LinearSolve, paclet:ref/LinearSolve]) will issue an error message:

A={{1.,2.},{5.,6.},{4.5,6.}};
B={5.,6.,8.};
LinearSolve[A,B]

In these cases it is possible to find a best-fit solution that minimizes ∥Ax-B(Subscript[∥, p]). A particularly common choice of p is 2; this generates a least squares solution. These are often used because the function (1/2)∥Ax-B(Subscript[∥, 2]^2) is differentiable in x and because the 2-norm is preserved under orthogonal transformations.

If the rank of A is n (i.e., it has full column rank), it can be shown (Golub and van Loan) that there is a unique solution to the least squares problem and it solves the linear system.

(A^T).Ax==(A^T)B

These are called the normal equations. Although they can be solved directly to get a least squares solution, this is not recommended because if the matrix A is ill-conditioned then the product (A^T).A will be even more ill-conditioned. Ways of finding least squares solutions for full rank matrices are explored in "Examples: Full Rank Least Squares Solutions".

It should be noted that best-fit solutions that minimize other norms, such as 1 and ∞, are also possible. A number of ways to find solutions to overdetermined systems that are suitable for special types of matrices or to minimize other norms are given under "Minimization of 1 and Infinity Norms".

A general way to find a least squares solution to an overdetermined system is to use a singular value decomposition to form a matrix that is known as the pseudoinverse of a matrix. In Mathematica this is computed with (RefLink[PseudoInverse, paclet:ref/PseudoInverse]). This technique works even if the input matrix is rank deficient. The basis of the technique follows.

Solve A.x==B
find p==(SuperPlus[A]) (pseudoinverse)
then x==p.B

Here is an overdetermined system; you can see that the matrix A is rank deficient:

A={{1.,2.},{2.,4.},{3.,6.}};
B={5.,6.,8.};
MatrixRank[A]

The least squares solution can still be found using (RefLink[PseudoInverse, paclet:ref/PseudoInverse]):

X=PseudoInverse[A].B

This shows how close the solution actually is:

A.X

Data Fitting

Scientific measurements often result in collections of ordered pairs of data, {(Subscript[x, i]),(Subscript[y, i])}. In order to make predictions at times other than those that were measured, it is common to try to fit a curve through the data points. If the curve is a straight line, then the aim is to find unknown coefficients α and β for which

(Subscript[y, i])==α+β(Subscript[x, i])

for all the data pairs. This can be written in a matrix form as shown here and is immediately recognized as being equivalent to solving an overdetermined system of linear equations.

(({{1, (Subscript[x, 1])}, {1, (Subscript[x, 2])}, {⋮, ⋮}, {1, (Subscript[x, n])}}))(({{α}, {β}}))==(({{(Subscript[y, 1])}, {(Subscript[y, 2])}, {⋮}, {(Subscript[y, n])}}))

Fitting a more general curve, for example,

(Subscript[y, i])==(Subscript[c, 0])+(Subscript[c, 1])(Subscript[x, i])+(Subscript[c, 2])(Subscript[x, i]^2)+…

is equivalent to the matrix equation

(({{1, (Subscript[x, 1]), (Subscript[x, 1]^2), …}, {1, (Subscript[x, 2]), (Subscript[x, 2]^2), …}, {⋮, ⋮, ⋮, …}, {1, (Subscript[x, n]), (Subscript[x, n]^2), …}}))(({{(Subscript[c, 0])}, {(Subscript[c, 1])}, {⋮}}))==(({{(Subscript[y, 1])}, {(Subscript[y, 2])}, {⋮}, {(Subscript[y, n])}})).

In this case the left-hand matrix is a Vandermonde matrix. In fact any function that is linear in the unknown coefficients can be used.

An example of how to find best-fit parameters will now be given. This will start with the measurements shown next:

data={{0.2,4.59},{0.4,5.05},{0.6,6.2},{0.8,6.6},{1.,7.4}};

These can be plotted graphically:

plot1=ListPlot[ data,PlotStyle->PointSize[0.02]]

You could try to solve this with the (RefLink[PseudoInverse, paclet:ref/PseudoInverse]) technique introduced in "PseudoInverse". First, split the data into x and y components:

xData=data[[All,1]];
yData=data[[All,2]];

Now form the left-hand side matrix:

mat=Transpose[ {Table[1,{Length[xData]}],xData}];
mat//MatrixForm

The best-fit parameter can be found as follows:

PseudoInverse[mat].yData

Mathematica provides the function (RefLink[FindFit, paclet:ref/FindFit]) that can find best-fit solutions to data. The solution agrees with that found by the (RefLink[PseudoInverse, paclet:ref/PseudoInverse]):

sol=FindFit[data,α+β x,{α,β},x]

Finally, a plot is made that shows the data points and the fitted curve:

Plot[ (α+β x)/.sol,{x,0,1},Epilog->First[plot1]]

The function (RefLink[FindFit, paclet:ref/FindFit]) is quite general; it can also fit functions to data that is not linear in the parameters.

Eigensystem Computations

The solution of the eigenvalue problem is one of the major areas for matrix computations. It has many applications in physics, chemistry, and engineering. For an n×n matrix A the eigenvalues are the n roots of its characteristic polynomial, p(z)==det(zI-A). The set of roots, λ(A)=={(Subscript[λ, 1]),…,(Subscript[λ, n])}, are called the spectrum of the matrix. For each eigenvalue, λ∈λ(A), the vectors, x, that satisfy

A.x==λ x

are described as eigenvectors. The matrix of the eigenvectors, X, if it exists and is nonsingular, may be used as a similarity transformation to form a diagonal matrix with the eigenvalues on the diagonal elements. Many important applications of eigenvalues involve the diagonalization of matrices in this way.

diag(λ(A))==(X^(-1)).A.X

Mathematica has various functions for computing eigenvalues and eigenvectors.

Eigenvalues and eigenvectors.

Here is a sample 2×2 matrix:

mat={{1., 2.},{3., 4.}};

This computes the eigenvalues:

Eigenvalues[mat]

This computes the eigenvectors:

Eigenvectors[mat]

You can compute both the eigenvalues and eigenvectors with (RefLink[Eigensystem, paclet:ref/Eigensystem]):

{vals,vecs}=Eigensystem[mat]

You can confirm that the first eigenvalue and its eigenvector satisfy the definition:

mat.vecs[[1]]-vals[[1]] vecs[[1]]

You should note that (RefLink[Eigenvectors, paclet:ref/Eigenvectors]) and (RefLink[Eigensystem, paclet:ref/Eigensystem]) return a list of eigenvectors. This means that if you want a matrix with columns that are the eigenvectors you must take a transpose. Here, the eigenvalues and their corresponding eigenvectors are shown to meet the definition of the eigenvectors:

mat.Transpose[vecs] - Transpose[vecs].DiagonalMatrix[vals]//Chop

You can also obtain the characteristic polynomial of the test matrix:

poly=CharacteristicPolynomial[mat,x]//Expand

This finds the roots of the characteristic polynomial; they are seen to match the eigenvalues of the matrix:

Solve[poly==0,x]

Eigenvalue computations work for sparse matrices as well as for dense matrices. In the following example, one of the eigenvalues is zero:

sp=SparseArray[ {i_,j_}/;Abs[i-j]<2->2.,{5,5}];
sp//MatrixForm

Eigenvalues[sp]

An eigenvalue being zero means that the null space of the matrix is not empty:

NullSpace[sp]

Certain applications of eigenvalues do not require all of the eigenvalues to be computed. Mathematica provides a mechanism for obtaining only some eigenvalues.

This is particularly useful if the matrix is very large. Here is a 10000×10000 sparse matrix:

sp=SparseArray[ {i_,j_}/;Abs[i-j]<4:>N[i+j],{10000,10000}]

This is the largest eigenvalue:

Eigenvalues[sp,1]

Options for (RefLink[Eigensystem, paclet:ref/Eigensystem]).

(RefLink[Eigensystem, paclet:ref/Eigensystem]) has four options that allow users more control. The (RefLink[Cubics, paclet:ref/Cubics]), (RefLink[Quartics, paclet:ref/Quartics]), and (RefLink[ZeroTest, paclet:ref/ZeroTest]) options are used by symbolic techniques and are discussed under "Symbolic and Exact Matrices". The (RefLink[Method, paclet:ref/Method]) option allows knowledgeable users to choose methods that are particularly appropriate for their problems; the default setting of (RefLink[Automatic, paclet:ref/Automatic]) means that (RefLink[Eigensystem, paclet:ref/Eigensystem]) makes a choice of a method that is generally suitable.

Eigensystem Properties

This section describes certain properties of eigensystem computations. It should be remembered that because the eigenvalues of an n×n matrix can be associated with the roots of an n(^th) degree polynomial, there will always be n eigenvalues.

Even if a matrix only contains real numbers, its eigenvalues may not be real:

mat={{1,-1},{1,1}};
Eigensystem[mat]

The eigenvalues of a matrix may be repeated; eigenvalues that are not repeated are called simple. In this example, the repeated eigenvalue has as many eigenvectors as the number of times it is repeated; this is a semisimple eigenvalue:

mat={{2,-1,-2},{2,5,2},{-2,-1,2}};
Eigensystem[mat]

The repeated eigenvalue may have fewer eigenvectors than the number of times it is repeated. In this case it is referred to as a defective eigenvalue and the matrix is referred to as a defective matrix:

mat={{-3,1,-3},{20,3,10},{2,-2,4}};
Eigensystem[mat]

A matrix may have a zero eigenvalue, but can still have a complete set of eigenvectors:

mat={{14,4,2},{-11,2,1},{14,4,2}};
Eigensystem[mat]

A matrix may have a zero eigenvalue, but not have a complete set of eigenvectors:

mat={{0,1},{0,0}};
Eigensystem[mat]

If a matrix is real and symmetric, all of the eigenvalues are real, and there is an orthonormal basis of eigenvectors:

mat={{1., 5.},{5., 1.}};
{lam,vec}=Eigensystem[mat]

The eigenvectors are an orthonormal basis:

vec.Transpose[vec]

A symmetric real matrix that has positive eigenvalues is known as positive definite:

mat={{10., 5.},{5., 16.}};
Eigensystem[mat]

A symmetric real matrix that has nonzero eigenvalues is known as positive semidefinite:

mat={{5,5},{5,5}};
Eigensystem[mat]

If a matrix is complex and Hermitian (i.e., equal to its conjugate transpose), all of the eigenvalues are real and there is an orthonormal basis of eigenvectors:

mat={{1.,5.I},{-5.I,10.}};
{lam,vec}=Eigensystem[mat]

The eigenvectors are an orthonormal basis:

vec.Transpose[Conjugate[vec]]

A Hermitian matrix that has positive eigenvalues is known as positive definite:

mat={{1.,5.I},{-5.I,30.}};
Eigensystem[mat]

Diagonalizing a Matrix

One of the uses of eigensystem computations is to provide a similarity transformation that diagonalizes a matrix:

(X^(-1)).A.X==diag(λ(A))

mat={{1., 2.},{3., 4.}};
vecs=Eigenvectors[mat];
Inverse[Transpose[vecs]].mat.Transpose[vecs]//Chop//MatrixForm

Similarity transformations preserve a number of useful properties of a matrix such as determinant, rank, and trace.

If the matrix of eigenvectors is singular then the matrix cannot be diagonalized. This happens if the matrix is defective (i.e., there is an eigenvalue that does not have a complete set of eigenvectors) as in this example:

mat={{-3,1,-3},{20,3,10},{2,-2,4}};
{lam,vec}=Eigensystem[mat];
vec=Transpose[vec];
vec//MatrixForm

The matrix of column eigenvectors is singular and so a similarity transformation that diagonalizes the original matrix does not exist:

Inverse[vec]

Symbolic and Exact Matrices

As is typical for Mathematica, if a computation is done for a symbolic or exact matrix it will use symbolic computer algebra techniques and return a symbolic or exact result.

This is a sample 3×3 matrix of integers:

mat={{1,2,3},{3,4,10},{-20,1,4}};

This computes the eigenvalues; the result uses (RefLink[Root, paclet:ref/Root]) objects:

Eigenvalues[mat]

A machine-precision approximation to the result can be computed by using (RefLink[N, paclet:ref/N]):

N[%]

The result can be returned in terms of radicals by setting the (RefLink[Cubics, paclet:ref/Cubics]) option to (RefLink[True, paclet:ref/True]):

Eigenvalues[mat,Cubics->True]

Generalized Eigenvalues

For n×n matrices A,B the generalized eigenvalues are the n roots of its characteristic polynomial, p(z)==det(A-zB). For each generalized eigenvalue, λ∈λ(A,B), the vectors, z, that satisfy

A.x==λ B.x

are described as generalized eigenvectors.

Generalized eigenvalues and eigenvectors.

Here are two 2×2 matrices:

matA={{1.746,0.940},{1.246,1.898}};
matB={{0.780,0.563},{0.913,0.659}};

This computes the generalized eigenvalues:

Eigenvalues[ {matA,matB}]

If the two matrices share a common nonzero null space, there will be an indeterminate eigenvalue:

matA={{1.5,0},{0,0}};
matB={{2.,0},{1.,0}};

These two matrices share a common one-dimensional null space:

NullSpace[matA]

NullSpace[matB]

One of the generalized eigenvalues is (RefLink[Indeterminate, paclet:ref/Indeterminate]):

Eigenvalues[ {matA,matB}]

Methods

Eigenvalue computations are solved with a number of different techniques that are specialized to particular problems. You can select between these using the option (RefLink[Method, paclet:ref/Method]). In this way a uniform interface is provided to all the functionality that Mathematica provides.

The default setting of the option (RefLink[Method, paclet:ref/Method]) is (RefLink[Automatic, paclet:ref/Automatic]). As is typical for Mathematica, this means that the system will make an automatic choice of method to use. For eigenvalue computation when the input is an n×n matrix of machine numbers and the number of eigenvalues requested is less than 20% of n an Arnoldi method is used. Otherwise, if the input matrix is numerical then a method using LAPACK routines is used. If the matrix is symbolic then specialized symbolic routines are used.

The details of the different methods are now described.

LAPACK

LAPACK is the default method for computing the entire set of eigenvalues and eigenvectors. When the matrix is unsymmetric, dgeev is used for real matrices and zgeev is used for complex matrices. For symmetric matrices, dsyevr is used for real matrices and zheevr is used for complex matrices. For generalized eigenvalues the routine dggev is used for real matrices and zggev for complex matrices. More information on LAPACK is available in the references section.

It should be noted that the computation of eigenvalues for a symmetric matrix is faster because Mathematica detects this case and switches to the appropriate method:

mat=Table[Random[],{700},{700}];
matSym=mat+Transpose[mat];

Eigenvalues[mat];//Timing

Eigenvalues[matSym];//Timing

If the input matrix uses arbitrary-precision numbers, LAPACK algorithms extended for arbitrary-precision computation are used.

Arnoldi

The Arnoldi method is an iterative method used to compute a finite number of eigenvalues. The implementation of the Arnoldi method uses the "ARPACK" library. The most general problem that can be solved with this technique is to compute a few selected eigenvalues for

A.x==λ B.x.

Because it is an iterative technique and only finds a few eigenvalues, it is often suitable for sparse matrices. One of the features of the technique is the ability to apply a shift σ and solve for a given eigenvector x and eigenvalue λ (where σ≠λ).

(((A-σB))^(-1)).B.x==x ν, where ν==1/(λ-σ)

This is effective for finding eigenvalues near to σ.

The following suboptions can be given to the Arnoldi method.

Suboptions of the Arnoldi method.

Here is a 3×3 sparse matrix:

sp=SparseArray[N[{{1,1,1},{1,1,1},{1,1,1}}]];

This computes the largest eigenvalue, specifying that the Arnoldi method is used:

Eigenvalues[sp,1,Method->Arnoldi]

This computes the eigenvalue near to -1:

Eigenvalues[sp,1,Method->{Arnoldi,Shift->3.00001}]

If a shift is given that matches an eigenvalue an error may result:

Eigenvalues[sp,1,Method->{Arnoldi,Shift->3.}]

This tries to find eigenvalues of a 1500×1500 tridiagonal matrix with 2 on the diagonal and -1 off the diagonal. The algorithm does not converge:

num=1500;
sp=SparseArray[{{i_,j_}/;Abs[i-j]≤1->2.-3.*Abs[i-j]},{num,num}];
Eigenvalues[sp,2]//Timing

When the technique does not converge it returns a result, but this may not be very accurate. One way to accelerate the convergence would be to use the results generated as a shift:

Eigenvalues[sp,2,Method->{"Arnoldi",Shift->4}]//Timing

An alternative way to get convergence is to increase the number of iterations:

Eigenvalues[sp,2,Method->{"Arnoldi",MaxIterations->100000}]//Timing

Alternatively, the convergence may be achieved by increasing the Arnoldi basis size in addition to the number of iterations. In this example, the increased basis size makes the algorithm converge in a smaller number of iterations. Even though each iteration is more expensive, the overall computation is faster:

Eigenvalues[sp,2,Method->{"Arnoldi",BasisSize->30,MaxIterations->100000}]//Timing

The Arnoldi algorithm with a basis size of 30 is still often faster and less memory intensive than the dense eigenvalue computation. (It should be noted that specific timings on individual machines can vary.)

mat=Normal[sp];
Eigenvalues[mat,2]//Timing

For many large sparse systems that occur in practical computations the Arnoldi algorithm is able to converge quite quickly.

Symbolic Methods

Symbolic eigenvalue computations work by interpolating the characteristic polynomial.

Matrix Decompositions

This section will discuss a number of standard techniques for working with matrices. These are often used as building blocks for solving matrix problems. The decompositions fall into the categories of factorizations, orthogonal transformations, and similarity transformations.

LU Decomposition

The LU decomposition of a matrix is frequently used as part of a Gaussian elimination process for solving a matrix equation. In Mathematica there is no need to use an LU decomposition to solve a matrix equation, because the function (RefLink[LinearSolve, paclet:ref/LinearSolve]) does this for you, as discussed in the section "Solving Linear Systems". Note that if you want to save the factorization of a matrix, so that it can be used to solve the same left-hand side, you can use the one-argument form of (RefLink[LinearSolve, paclet:ref/LinearSolve]), as demonstrated in the section "Saving the Factorization".

The LU decomposition with partial pivoting of a square nonsingular matrix A involves the computation of the unique matrices L, U, and P such that

P.A==L.U

where L is lower triangular with ones on the diagonal, U is upper triangular, and P is a permutation matrix. Once it is found, it can be used to solve the matrix equation by finding the solution to the following equivalent system.

L.U.x==P.B

This can be done by first solving for y in the following.

L.y==P.B

Then solve for x to get the desired solution.

U.x==y

Because both L and U are triangular matrices it is particularly easy to solve these systems.

Despite the fact that in Mathematica you do not need to use the LU decomposition for solving matrix equations, Mathematica provides the function (RefLink[LUDecomposition, paclet:ref/LUDecomposition]) for computing the LU decomposition.

(RefLink[LUDecomposition, paclet:ref/LUDecomposition]) returns a list of three elements. The first element is a combination of upper and lower triangular matrices, the second element is a vector specifying rows used for pivoting (a permutation vector which is equivalent to the permutation matrix), and the third element is an estimate of the condition number.

In this example, the three results of the (RefLink[LUDecomposition, paclet:ref/LUDecomposition]) computation are shown:

mat={{5.,2.,8.},{1.,2.,3.},{3,4.,1.}};
{mLU,vecP,cond}=LUDecomposition[mat]

One of the features of the LU decomposition is that the lower and upper matrices are stored in the same matrix. The individual L and U factors can be obtained as follows. Note how they use a vectorized operation for element-by-element multiplication with a sparse matrix. This makes the process more efficient; these techniques are discussed further in "Programming Efficiency".

This generates the L matrix:

mL=Transpose[mLU]; num=Length[mL];
sp=SparseArray[ {i_,j_}/;i≥j->1,{num,num}];
mL=sp mL;
Do[mL[[i,i]]=1,{i,num}];
mL//MatrixForm

This generates the U matrix:

mU=Transpose[mLU]; num=Length[mU];
sp=SparseArray[ {i_,j_}/;i≤j->1,{num,num}];
mU=mU sp;
mU//MatrixForm

The original matrix was nonsingular and so this is a factorization. Therefore multiplying L and U recreates the original matrix permuted with the row pivot vector:

mL.mU

The following subtracts the product of the L and U matrices from the permuted original matrix. The difference is very close to zero:

mat[[vec]]-%

The diagonal elements in the U matrix are known as the pivots. If any zero pivots emerge during a computation of the LU factorization this means that the matrix is singular. In such a case (RefLink[LUDecomposition, paclet:ref/LUDecomposition]) produces a warning:

mat = {{1.,2.,3.},{1.,2.,3.},{1.,2.,3.}};
LUDecomposition[mat]

Cholesky Decomposition

The Cholesky decomposition of a symmetric positive definite matrix A is a factorization into a unique upper triangular G such that

A==(G^T).G.

This factorization has a number of uses, one of which is that, because it is a triangular factorization, it can be used to solve systems of equations involving symmetric positive definite matrices. (The way that a triangular factorization can be used to solve a matrix equation is shown in the section "LU Decomposition".) If a matrix is known to be of this form it is preferred over the LU factorization because the Cholesky factorization is faster to compute. If you want to solve a matrix equation using the Cholesky factorization you can do this directly from (RefLink[LinearSolve, paclet:ref/LinearSolve]) using the Cholesky method; this is described in a previous section.

The Cholesky factorization can be computed in Mathematica with the function (RefLink[CholeskyDecomposition, paclet:ref/CholeskyDecomposition]).

mat={{10.,5.},{5.,20.}};
matG=CholeskyDecomposition[mat];
matG//MatrixForm

This reconstructs the original matrix:

Transpose[matG].matG

If the matrix is not positive definite then the Cholesky decomposition does not exist:

mat1={{10.,5.},{5.,1.}};
CholeskyDecomposition[mat1]

There are a number of equivalent definitions of positive definite, such as the eigenvalues being all positive:

mat={{10.,5.},{5.,20.}};
Eigenvalues[mat]

One way to test if a matrix is positive definite is to see if the Cholesky decomposition exists.

Note that if the matrix is complex, the factorization is defined by the conjugate transpose. The following computes the Cholesky decomposition of a complex matrix:

mat={{11.,3.+1.I},{3.-1.I,3.}};
matG=CholeskyDecomposition[mat]

This demonstrates that the factorization is correct:

Conjugate[Transpose[matG]].matG//MatrixForm

Cholesky and LU Factorizations

The Cholesky factorization is related to the LU factorization as

(G^T)==L.(Sqrt[D])

where D is the diagonal matrix of pivots.

This can be demonstrated as follows:

mat={{10.,5.},{5.,20.}};
{mLU,vec,cond}=LUDecomposition[mat];
matG=CholeskyDecomposition[mat];
matG//MatrixForm

Now you can compute the L matrix:

mL=Transpose[mLU]; num=Length[mL];
sp=SparseArray[ {i_,j_}/;i≥j->1,{num,num}];
mL=sp mL;
Do[mL[[i,i]]=1,{i,num}];
mL//MatrixForm

This step computes the matrix of pivots:

sp=SparseArray[ {i_,i_}->1,{num,num}];
mD=mLU sp;
mD//MatrixForm

Finally, this computes L.(Sqrt[D]); its transpose is equal to the Cholesky decomposition.

(mL.Sqrt[mD])//Transpose//MatrixForm

Orthogonalization

Orthogonalization generates an orthonormal basis from an arbitrary basis. An orthonormal basis is a basis, {(Subscript[ν, 1]),(Subscript[ν, 2]),…}, for which

(Subscript[ν, i]).(Subscript[ν, j])=={({{1i==j}, {0i≠j}}).

In Mathematica a set of vectors can be orthogonalized with the function (RefLink[Orthogonalize, paclet:ref/Orthogonalize]).

This creates a set of three vectors that form a basis for (ℝ^3):

w1={1.,1,1.};
w2={1,1.2,0.9};
w3={1.,1.,0.8};

A plot visualizes the vectors; they all tend to lie in the same direction:

Show[Graphics3D[ {Hue[0],Map[ Line[ {{0,0,0},#}]&,{w1,w2,w3}]}],Axes->True]

This computes an orthonormal basis:

{v1,v2,v3}=Orthogonalize[{w1,w2,w3}]

The orthonormal vectors are obviously much more spread out:

Show[Graphics3D[ {Hue[0],Map[ Line[ {{0,0,0},#}]&,{v1,v2,v3}]}],Axes->True]

The vectors v1, v2, and v3 are orthonormal; thus, the dot product of each vector with itself is 1:

Map[#.#&,{v1,v2,v3}]

In addition, the dot product of a vector with another vector is 0:

v1.v2

This uses (RefLink[Outer, paclet:ref/Outer]) to compare all vectors with all other vectors:

Outer[ Dot,{v1,v2,v3},{v1,v2,v3},1]//Chop

By default a Gram–Schmidt orthogonalization is computed, but a number of other orthogonalizations can be computed.

Orthogonalization methods.

This uses a Householder method:

Orthogonalize[{w1,w2,w3},Method->Householder]

QR Decomposition

The QR decomposition of a rectangular matrix A with linearly independent columns involves the computation of matrices Q and R such that

A==Q . R

where the columns of Q form an orthonormal basis for the columns of A and R is upper triangular. It can be computed in Mathematica with the function (RefLink[QRDecomposition, paclet:ref/QRDecomposition]).

This computes the QR decomposition of a sample matrix. The result is a list of the two matrices:

mat={{1.,3.},{4.,8.}};
{matQ, matR}=QRDecomposition[ mat]

In fact the first argument that is returned is a list of the orthonormal columns. To generate the Q matrix it is necessary to compute the transpose:

matQ=Transpose[matQ];
{matQ//MatrixForm, matR//MatrixForm}

This is a matrix factorization, so Q:R is equal to A:

matQ.matR

In addition, the columns of Q are orthonormal. Because the result of (RefLink[QRDecomposition, paclet:ref/QRDecomposition]) returns the list of columns, they can be extracted with one part index. The following two examples demonstrate the orthonormal properties of the columns:

matQ[[All,1]].matQ[[All,2]]

matQ[[All,2]].matQ[[All,2]]

The QR decomposition is also defined for rectangular matrices. In this case m<n:

mat={{1.,7.,3.},{4., 5., 1.}};
{matQ, matR}=QRDecomposition[mat];matQ=Transpose[matQ];
{matQ//MatrixForm, matR//MatrixForm}

For an input matrix where m≤n, the result of (RefLink[QRDecomposition, paclet:ref/QRDecomposition]) is two matrices: an m×m matrix Q and an m×n matrix R. When the matrix represents an overdetermined system of equations, that is, m>n, (RefLink[QRDecomposition, paclet:ref/QRDecomposition]) returns an m×n matrix Q and an n×n matrix R. (Note that to generate the matrix Q, it is necessary to transpose the result from (RefLink[QRDecomposition, paclet:ref/QRDecomposition]).) This is demonstrated in the following example:

mat={{1.,7.},{3.,4.},{ 5., 1.}};
{matQ, matR}=QRDecomposition[mat];matQ=Transpose[matQ];
{matQ//MatrixForm, matR//MatrixForm}

This is often referred to as the thin QR decomposition; see for example Golub and van Loan. If you wish to compute the full QR decomposition, it is possible to pad out R with extra rows of zeros and Q with extra orthonormal columns. This can be done with the following function:

FullQRDecomposition[ mat_]:=
Module[ {q,r,n,m},
{q,r}=QRDecomposition[mat];
{n,m}=Dimensions[mat];
{ Join[q,NullSpace[q]], 
Join[r,Table[0,{n-Length[r]},{m}]]}
]

The Q and R computed by this function are m×m and m×n, respectively:

mat={{1.,7.},{3.,4.},{ 5., 1.}};
{matQ, matR}=FullQRDecomposition[mat];matQ=Transpose[matQ];
{matQ//MatrixForm, matR//MatrixForm}

matQ.matR

Solving Systems of Equations

For a square nonsingular matrix A the QR decomposition can be used to solve the matrix equation A.x==B, as is also the case for the LU decomposition. However, when the matrix is rectangular, the QR decomposition is also useful for solving the matrix equation.

One particular application of the QR decomposition is to find least squares solutions to overdetermined systems, by solving the system of normal equations

(A^T).A.x==(A^T).B.

Because (Q^T).Q==I, this can be simplified as

(A^T).A==(((Q.R))^T).(Q.R)==(R^T).R.

Thus the normal equations simplify to

(R^T).R.x==(R^T).(Q^T).B

and because (R^T) is nonsingular this simplifies to

R.x==(Q^T).b.

Because R is triangular this is a particularly easy system to solve; sample code to implement this technique is given in "Examples: Least Squares QR".

Singular Value Decomposition

The singular value decomposition of a rectangular matrix A involves the computation of orthogonal matrices U and V and a diagonal matrix D such that

A==U.D.(V^T).

The diagonal elements of the matrix D are called the singular values of A. In Mathematica, the function (RefLink[SingularValueList, paclet:ref/SingularValueList]) computes the singular values and the function (RefLink[SingularValueDecomposition, paclet:ref/SingularValueDecomposition]) computes the full singular value decomposition.

mat={{1., 2.},{6., 1.}};
SingularValueList[mat]

{matU,matD,matV}=SingularValueDecomposition[mat];
{matU//MatrixForm,matD//MatrixForm,matV//MatrixForm}

This is a factorization so the original matrix can be restored:

matU.matD.Transpose[matV]//MatrixForm

Because the matrices U and V are orthogonal, they can be used as an orthogonal transformation of the original matrix to generate the diagonal matrix with the singular values on the diagonal:

Transpose[matU].mat.matV//MatrixForm

If the matrix is singular then some of the singular values will be zero:

mat={{1., 2.},{1., 2.}};
{matU,matD,matV}=SingularValueDecomposition[mat];
matD//MatrixForm

(RefLink[SingularValueList, paclet:ref/SingularValueList]) only returns the nonzero singular values:

SingularValueList[mat]

Note that if the matrix is complex, the definition of the singular value decomposition uses the conjugate transpose.

There are many applications of the singular value decomposition. The singular values of a matrix give information on the linear transformation represented by A. For example, the action of A on a unit sphere generates an ellipsoid with semi-axes given by the singular values. The singular values can be used to compute the rank of a matrix; the number of nonzero singular values is equal to the rank. The singular values can be used to compute the 2-norm of a matrix, and the columns of the V matrix that correspond to zero singular values are the null space of the matrix.

Certain applications of singular values do not require all of the singular values to be computed. Mathematica provides a mechanism for obtaining only some singular values.

This returns the smallest singular value of the input matrix; because it is zero, this demonstrates the matrix is not full rank:

mat={{1., 2.},{1., 2.}};
SingularValueList[mat,-1]

Generalized Singular Values

For an m×n matrix A and a p×n matrix B the generalized singular values are given by the pair of factorizations

A==U.(Subscript[Σ, 1]).(Q^T) and B==V.(Subscript[Σ, 2]).(Q^T)

where U is m×m, V is p×p, and Q is n×n; U and V are orthogonal, and Q is invertible.

These are some sample matrices:

matA={{1.,2.},{4.,3.}};
matB={{4.,1.},{3.,2.}};

This returns the generalized singular values of the matrices matA and matB:

SingularValueList[ {matA,matB}]

This computes the full generalized singular value decomposition of the matrices matA and matB:

{{matU,matV},{dU,dV},matQ}=SingularValueDecomposition[ {matA,matB}];

This demonstrates the identity for matA:

matU.dU.Transpose[matQ]//MatrixForm

This demonstrates the identity for matB:

matV.dV.Transpose[matQ]//MatrixForm

Finally, the singular values are computed by dividing the corresponding diagonal elements:

Tr[dU,List]/Tr[dV,List]

Options

The functions (RefLink[SingularValueList, paclet:ref/SingularValueList]) and (RefLink[SingularValueDecomposition, paclet:ref/SingularValueDecomposition]) both take a (RefLink[Tolerance, paclet:ref/Tolerance]) option.

The option controls the size at which singular values are treated as zero. By default, values that are tol times smaller than the largest singular value are dropped, where tol is 100×(RefLink[$MachineEpsilon, paclet:ref/$MachineEpsilon]) for machine-number matrices. For arbitrary-precision matrices, it is (10^(-p)), where p is the precision of the matrix.

The smallest singular value of this matrix is just larger than the default setting for the tolerance. It is not treated as being equivalent to a zero singular value:

mat={{1.,2.},{2.,4.-6 (10^(-13))}};
num=SingularValueList[mat]

Increasing the setting for the tolerance causes the small singular value to be treated as zero:

SingularValueList[mat,Tolerance->(10^(-12))]

Schur Decomposition

The Schur decomposition of a square matrix A involves finding a unitary matrix U that can be used for a similarity transformation of A to form a block upper triangular matrix T with 1×1 and 2×2 blocks on the diagonal (the 2×2 blocks correspond to complex conjugate pairs of eigenvalues for a real matrix A). A block upper triangular matrix of this form can be called upper quasi-triangular.

(U^T).A.U==T

The Schur decomposition always exists and so the similarity transformation of A to upper triangular always exists. This contrasts with the eigensystem similarity transformation, used to diagonalize a matrix, which does not always exist.

This computes the Schur decomposition of a sample matrix. The result is a list of the two matrices:

mat={{1., 4.},{2.,3.}};
{matU,matT}=SchurDecomposition[mat];
{matU//MatrixForm, matT//MatrixForm}

The matrix U can be used for a similarity transformation on the matrix to generate the upper triangular result T:

Transpose[matU].mat.matU//MatrixForm

For this particular matrix, a similarity transformation that generates an even simpler form can be found, because the matrix can be diagonalized:

matV=Transpose[Eigenvectors[mat]];
invV=Inverse[matV];
invV.mat.matV//MatrixForm

This matrix cannot be diagonalized because the matrix of eigenvectors is singular:

matNoDiag={{0.,0.,1.},{1.,0.,0.},{0.,0.,0.}};
matV=Transpose[Eigenvectors[matNoDiag]];
invV = Inverse[ matV]

However, the Schur decomposition can be found and the matrix transformed to an upper triangular form:

{matQ,matT}=SchurDecomposition[matNoDiag];
{matQ//MatrixForm,matT//MatrixForm}

In this example, the matrix has complex eigenvalues:

mat={{1.,-1.},{1.,1.}};
Eigenvalues[mat]

Now the resulting matrix T has a 2×2 block on the diagonal:

{matU,matT}=SchurDecomposition[mat];
{matU//MatrixForm,matT//MatrixForm}

The matrix U can be used for a similarity transformation on the matrix to generate the upper quasi-triangular result. Note that an upper triangular result (with 1×1 blocks on the diagonal) that may involve complex numbers can be obtained by using the option (RefLinkPlain[RealBlockForm, paclet:tutorial/LinearAlgebraMatrixComputations#35657259]). When the result is upper triangular (i.e., has 1×1 blocks on the diagonal) the eigenvalues of the matrix are always found on the diagonal:

Transpose[matU].mat.matU//MatrixForm

Note that if the matrix is complex the definition of the Schur decomposition uses the conjugate transpose and returns an upper triangular result. This is demonstrated for the following complex matrix:

mat={{2., 4.+4.I},{1.,2.}};
{matU,matT}=SchurDecomposition[mat];
{matU//MatrixForm, matT//MatrixForm}

This demonstrates that the result satisfies the definition of the Schur decomposition:

matU.matT.Conjugate[Transpose[matU]]//Chop//MatrixForm

The diagonal of T contains the eigenvalues of A:

Eigenvalues[mat]

Generalized Schur Decomposition

For n×n matrices A and B, the generalized Schur decomposition is defined as

A==Q.S.(Z^T) and B==Q.T.(Z^T)

where Q and Z are unitary, T is upper triangular, and S is upper quasi-triangular.

These are some sample matrices:

matA={{1.,2.},{4.,3.}};
matB={{4.,1.},{3.,2.}};

This returns the generalized Schur decomposition:

{matQ,matS,matZ,matT}=SchurDecomposition[ {matA,matB}]

This demonstrates the results are consistent with the definition of the decomposition:

matQ.matS.Transpose[matZ]//MatrixForm

matQ.matT.Transpose[matZ]//MatrixForm

For real input, a result involving complex numbers and an upper triangular result can be obtained with the option RealBlockForm.

Options

(RefLink[SchurDecomposition, paclet:ref/SchurDecomposition]) takes two options.

Options for (RefLink[SchurDecomposition, paclet:ref/SchurDecomposition]).

The option Pivoting can be used to allow pivoting and scaling to improve the quality of the result. When it is set to (RefLink[True, paclet:ref/True]), pivoting and scaling may be used and a matrix P that represents the changes to A is returned. With this new matrix the definition of the Schur decomposition can be seen as follows.

(U^T).(P^(-1)).A.P.U==T

The use of pivoting and scaling is now demonstrated for the following matrix. When Pivoting is set to (RefLink[True, paclet:ref/True]), pivoting and scaling are used if necessary, and an extra matrix (which here only represents scaling) is returned:

mat=N[{{10^16,100,1},{2 10^16,115,2},{3 10^16,120,3}}];
{matU,matT,matP}=SchurDecomposition[mat, Pivoting->True];
{matU//MatrixForm,matT//MatrixForm,matP//MatrixForm}

This demonstrates the transformation to an upper triangular form:

Transpose[matU].Inverse[matP].mat.matP.matU//MatrixForm

The diagonal elements of T are the eigenvalues of A:

Eigenvalues[mat]

When the Schur decomposition is computed without pivoting and scaling, the diagonal elements of T are not as close to the eigenvalues of A. This demonstrates the utility of the Pivoting option:

{matU,matT}=SchurDecomposition[mat];
matT//MatrixForm

The option RealBlockForm controls the generation of the upper quasi-triangular result. If this is set to (RefLink[True, paclet:ref/True]), a result that may have 2×2 blocks on the diagonal is generated. If it is set to (RefLink[False, paclet:ref/False]), the result is an upper triangular with 1×1 blocks on the diagonal (but which may involve complex numbers). This is demonstrated for the following real matrix, which has complex eigenvalues and a Schur decomposition with 2×2 blocks on the diagonal:

mat={{1.0000000000000004,-3.302775637731996},{0.3027756377319948,1.0000000000000004}};
Eigenvalues[mat]

Setting RealBlockForm to (RefLink[False, paclet:ref/False]) generates a matrix T that is upper triangular; both U and T are complex:

{matU,matT}=SchurDecomposition[mat, RealBlockForm->False];
{matU//MatrixForm,matT//MatrixForm}

This generates the upper triangular result T:

Conjugate[Transpose[matU]].mat.matU//MatrixForm//Chop

Jordan Decomposition

The Jordan decomposition of a square matrix A involves finding the nonsingular matrix X that can be used for a similarity transformation of A to generate a matrix J (known as the Jordan form) that has a particularly simple triangular structure.

(X^(-1)).A.X==J

The Jordan decomposition always exists, but it is hard to compute with floating-point arithmetic. However, computation with exact arithmetic avoids these problems.

This demonstrates the Jordan form for a sample matrix:

mat={{5,4,0,0,4,3},{2,3,1,0,5,1},{0,-1,2,0,2,0},{-8,-8,-1,2,-12,-7},{0,0,0,0,-1,0},{-8,-8,-1,0,-9,-5}};
{matX,matJ}=JordanDecomposition[mat];
matX.matJ.Inverse[matX]//MatrixForm

The Jordan form has the eigenvalues of the matrix along its diagonal. Any defective eigenvalues are grouped into blocks by 1s just above the diagonal. The Jordan form of the above matrix is shown here:

matJ//MatrixForm

The Jordan form shows that there are two eigenvalues: -1 and 2. The eigenvalue -1 is repeated twice and has a complete set of eigenvectors. The eigenvalue 2 is repeated four times. It appears once with its own eigenvector, and then three times with only one full eigenvector. This is demonstrated when the eigensystem for the matrix is computed:

Eigensystem[mat]

Functions of Matrices

The computation of functions of matrices is a general problem with applications in many areas such as control theory. The function f(A) of a square matrix A is not the same as the application of the function to each element in the matrix. Clearly element-wise application would not maintain properties consistent with the application of the function to a scale.

For example, each element of the following matrix is raised to the zero power:

mat={{1.,2.},{4.,5.}};
(mat^0)

However, a much better result would be the identity matrix. Mathematica has a function for raising a matrix to a power, and when a matrix is raised to the zero power the result is the identity matrix:

MatrixExp[mat,0]

There are a number of ways to define functions of matrices; one useful way is to consider a series expansion. For the exponential function this works as follows.

(e^A)==I+A+((A^2)/(2!))+((A^3)/(3!))+…

One way to compute this series involves diagonalizing A, so that A==P.D.(P^(-1)) and (A^k)==P.(D^k).(P^(-1)). Therefore, the exponential of A can be computed as follows.

(e^A)==(∑_(k==0)^∞)((A^k)/(k!))==P.((∑_(k==0)^∞)((D^k)/(k!))).(P^(-1))==P.diag((e^(Subscript[λ, 1])),…,(e^(Subscript[λ, n]))).(P^(-1))

This technique can be generalized to functions f(λ) of the eigenvalues of A. Note that while this is one way to define functions of matrices, it does not provide a good way to compute them.

Mathematica does not have a function for computing general functions of matrices, but it has some specific functions.

Here is a sample matrix:

mat={{1., 2.},{5.,1.}};

This raises the matrix to the fourth power:

MatrixPower[mat,4]

The result is equivalent to squaring the square of the matrix:

matSq=mat.mat;matSq.matSq

This computes the exponential of the matrix:

MatrixExp[mat]

It is equivalent to the computation that uses the eigensystem of the matrix. (It should be noted that this is not an efficient way to compute a function of a matrix; the example here is only for exposition.)

{eigs,vecs}=Eigensystem[mat];
vecs=Transpose[vecs];
vecsInv=Inverse[vecs];
vecs.DiagonalMatrix[Exp[eigs]].vecsInv

A technique for computing parametrized functions of matrices by solving differential equations is given in "Examples: Matrix Functions with NDSolve".