Learn more about MATLAB

Systems of Linear Equations

On this page…
Computational Considerations The mldivide Algorithm General Solution Square Systems Overdetermined Systems Using Multithreaded Computation with Systems of Linear Equations Iterative Methods for Solving Systems of Linear Equations

Computational Considerations

One of the most important problems in technical computing is the solution of systems of simultaneous linear equations.

In matrix notation, the general problem takes the following form: Given two matrices A and B, does there exist a unique matrix X so that AX = B or XA = B?

It is instructive to consider a 1-by-1 example. For example, does the equation

have a unique solution?

The answer, of course, is yes. The equation has the unique solution x = 3. The solution is easily obtained by division:

The solution is not ordinarily obtained by computing the inverse of 7, that is 7^-1 = 0.142857..., and then multiplying 7^-1 by 21. This would be more work and, if 7^-1 is represented to a finite number of digits, less accurate. Similar considerations apply to sets of linear equations with more than one unknown; the MATLAB software solves such equations without computing the inverse of the matrix.

Although it is not standard mathematical notation, MATLAB uses the division terminology familiar in the scalar case to describe the solution of a general system of simultaneous equations. The two division symbols, slash, /, and backslash, \, correspond to the two MATLAB functions mldivide and mrdivide. mldivide and mrdivide are used for the two situations where the unknown matrix appears on the left or right of the coefficient matrix:

X = B/A	Denotes the solution to the matrix equation XA = B.
X = A\B	Denotes the solution to the matrix equation AX = B.

Think of "dividing" both sides of the equation AX = B or XA = B by A. The coefficient matrix A is always in the "denominator."

The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows. The solution X then has the same number of columns as B and its row dimension is equal to the column dimension of A. For X = B/A, the roles of rows and columns are interchanged.

In practice, linear equations of the form AX = B occur more frequently than those of the form XA = B. Consequently, the backslash is used far more frequently than the slash. The remainder of this section concentrates on the backslash operator; the corresponding properties of the slash operator can be inferred from the identity:

(B/A)' = (A'\B')

The coefficient matrix A need not be square. If A is m-by-n, there are three cases:

m = n	Square system. Seek an exact solution.
m > n	Overdetermined system. Find a least squares solution.
m < n	Underdetermined system. Find a basic solution with at most m nonzero components.

Back to Top

The mldivide Algorithm

The mldivide operator employs different algorithms to handle different kinds of coefficient matrices. The various cases are diagnosed automatically by examining the coefficient matrix.

Permutations of Triangular Matrices

mldivide checks for triangularity by testing for zero elements. If a matrix A is triangular, MATLAB software uses a substitution to compute the solution vector x. If A is a permutation of a triangular matrix, MATLAB software uses a permuted substitution algorithm.

Square Matrices

If A is symmetric and has real, positive diagonal elements, MATLAB attempts a Cholesky factorization. If the Cholesky factorization fails, MATLAB performs a symmetric, indefinite factorization. If A is upper Hessenberg, MATLAB uses Gaussian elimination to reduce the system to a triangular matrix. If A is square but is neither permuted triangular, symmetric and positive definite, or Hessenberg, MATLAB performs a general triangular factorization using LU factorization with partial pivoting (see lu).

Rectangular Matrices

If A is rectangular, mldivide returns a least-squares solution. MATLAB solves overdetermined systems with QR factorization (see qr). For an underdetermined system, MATLAB returns the solution with the maximum number of zero elements.

The mldivide function reference page contains a more detailed description of the algorithm.

Back to Top

General Solution

The general solution to a system of linear equations AX = b describes all possible solutions. You can find the general solution by:

1. Solving the corresponding homogeneous system AX = 0. Do this using the null command, by typing null(A). This returns a basis for the solution space to AX = 0. Any solution is a linear combination of basis vectors.

2. Finding a particular solution to the nonhomogeneous system AX = b.

You can then write any solution to AX = b as the sum of the particular solution to AX = b, from step 2, plus a linear combination of the basis vectors from step 1.

The rest of this section describes how to use MATLAB to find a particular solution to AX = b, as in step 2.

Back to Top

Square Systems

The most common situation involves a square coefficient matrix A and a single right-hand side column vector b.

Nonsingular Coefficient Matrix

If the matrix A is nonsingular, the solution, x = A\b, is then the same size as b. For example:

A = pascal(3);
u = [3; 1; 4];
x = A\u

x =
      10
     -12
       5

It can be confirmed that A*x is exactly equal to u.

If A and B are square and the same size, X = A\B is also that size:

B = magic(3);
X = A\B

X =
      19    -3    -1
     -17     4    13
       6     0    -6

It can be confirmed that A*X is exactly equal to B.

Both of these examples have exact, integer solutions. This is because the coefficient matrix was chosen to be pascal(3), which has a determinant equal to 1.

Singular Coefficient Matrix

A square matrix A is singular if it does not have linearly independent columns. If A is singular, the solution to AX = B either does not exist, or is not unique. The backslash operator, A\B, issues a warning if A is nearly singular and raises an error condition if it detects exact singularity.

If A is singular and AX = b has a solution, you can find a particular solution that is not unique, by typing

P = pinv(A)*b

P is a pseudoinverse of A. If AX = b does not have an exact solution, pinv(A) returns a least-squares solution.

For example:

A = [ 1     3     7
     -1     4     4
      1    10    18 ]

is singular, as you can verify by typing

det(A)

ans =
     0

Note   For information about using pinv to solve systems with rectangular coefficient matrices, see Pseudoinverses.

Exact Solutions.   For b =[5;2;12], the equation AX = b has an exact solution, given by

pinv(A)*b

ans =
    0.3850
   -0.1103
    0.7066

Verify that pinv(A)*b is an exact solution by typing

A*pinv(A)*b

ans =
    5.0000
    2.0000
   12.0000

Least-Squares Solutions.   On the other hand, if b = [3;6;0], AX = b does not have an exact solution. In this case, pinv(A)*b returns a least squares solution. If you type

A*pinv(A)*b

ans =
   -1.0000
    4.0000
    2.0000

you do not get back the original vector b.

You can determine whether AX = b has an exact solution by finding the row reduced echelon form of the augmented matrix [A b]. To do so for this example, enter

rref([A b])
ans =
    1.0000         0    2.2857         0
         0    1.0000    1.5714         0
         0         0         0    1.0000

Since the bottom row contains all zeros except for the last entry, the equation does not have a solution. In this case, pinv(A) returns a least-squares solution.

Back to Top

Overdetermined Systems

Overdetermined systems of simultaneous linear equations are often encountered in various kinds of curve fitting to experimental data. Here is a hypothetical example. A quantity y is measured at several different values of time, t, to produce the following observations:

t	y
0.0	0.82
0.3	0.72
0.8	0.63
1.1	0.60
1.6	0.55
2.3	0.50

Enter the data into MATLAB with the statements

t = [0 .3 .8 1.1 1.6 2.3]';
y = [.82 .72 .63 .60 .55 .50]';

Try modeling the data with a decaying exponential function:

The preceding equation says that the vector y should be approximated by a linear combination of two other vectors, one the constant vector containing all ones and the other the vector with components e^–t. The unknown coefficients, c₁ and c₂, can be computed by doing a least-squares fit, which minimizes the sum of the squares of the deviations of the data from the model. There are six equations in two unknowns, represented by the 6-by-2 matrix:

E = [ones(size(t)) exp(-t)]

E =
      1.0000    1.0000
      1.0000    0.7408
      1.0000    0.4493
      1.0000    0.3329
      1.0000    0.2019
      1.0000    0.1003

Use the backslash operator to get the least-squares solution:

c = E\y

c =
      0.4760
      0.3413

In other words, the least-squares fit to the data is

The following statements evaluate the model at regularly spaced increments in t, and then plot the result, together with the original data:

T = (0:0.1:2.5)';
Y = [ones(size(T)) exp(-T)]*c;
plot(T,Y,'-',t,y,'o')

E*c is not exactly equal to y, but the difference might well be less than measurement errors in the original data.

A rectangular matrix A is rank deficient if it does not have linearly independent columns. If A is rank deficient, the least-squares solution to AX = B is not unique. The backslash operator, A\B, issues a warning if A is rank deficient and produces a least-use the format command to display the solution in rational format. The particular solution is obtained with

format rat
p = R\b
p =
       0
      -3/7
       0
      29/7

One of the nonzero components is p(2) because R(:,2) is the column of R with largest norm. The other nonzero component is p(4) because R(:,4) dominates after R(:,2) is eliminated.

The complete solution to the underdetermined system can be characterized by adding an arbitrary vector from the null space, which can be found using the null function with an option requesting a rational basis:

Z = null(R,'r')
Z =
      -1/2           -7/6     
      -1/2            1/2     
       1             0      
       0             1 

It can be confirmed that R*Z is zero and that any vector x where

x = p + Z*q

for an arbitrary vector q satisfies R*x = b.

Back to Top

Using Multithreaded Computation with Systems of Linear Equations

MATLAB software supports multithreaded computation for a number of linear algebra and element-wise numerical functions. These functions automatically execute on multiple threads. For a function or expression to execute faster on multiple CPUs, a number of conditions must be true:

1. The function performs operations that easily partition into sections that execute concurrently. These sections must be able to execute with little communication between processes. They should require few sequential operations.

2. The data size is large enough so that any advantages of concurrent execution outweigh the time required to partition the data and manage separate execution threads. For example, most functions speed up only when the array contains than several thousand elements or more.

3. The operation is not memory-bound; processing time is not dominated by memory access time. As a general rule, complex functions speed up more than simple functions.

inv, lscov, linsolve, and mldivide show significant increase in speed on large double-precision arrays (on order of 10,000 elements or more) when multithreading is enabled. To learn more about mulithreading, see MATLAB Multiprocessing in the MATLAB Programming Fundamentals documentation.

Back to Top

Iterative Methods for Solving Systems of Linear Equations

If the coefficient matrix A is large and sparse, factorization methods are generally not efficient. Iterative methods generate a series of approximate solutions. MATLAB provides several iterative methods to handle large, sparse input matrices.

pcg

Preconditioned conjugate gradients method. This method is appropriate for Hermitian positive definite coefficient matrix A.

bicg

BiConjugate Gradients Method

bicgstab

BiConjugate Gradients Stabilized Method

bicgstabl

BiCGStab(l) Method

cgs

Conjugate Gradients Squared Method

gmres

Generalized Minimum Residual Method

lsqr

LSQR Method

minres

Minimum Residual Method. This method is appropriate for Hermitian coefficient matrix A.

qmr

Quasi-Minimal Residual Method

symmlq

Symmetric LQ Method

tfqmr

Transpose-Free QMR Method

Back to Top

Matrices in the MATLAB Environment

Inverses and Determinants

© 1984-2009- The MathWorks, Inc.    -   Site Help   -   Patents   -   Trademarks   -   Privacy Policy   -   Preventing Piracy   -   RSS