Gaussian elimination
Encyclopedia
In linear algebra
, Gaussian elimination is an algorithm
for solving systems of linear equations. It can also be used to find the rank
of a matrix
, to calculate the determinant
of a matrix, and to calculate the inverse of an invertible square matrix. The method is named after Carl Friedrich Gauss
, but it was not invented by him.
Elementary row operations
are used to reduce a matrix to what is called triangular form (in numerical analysis) or row echelon form
(in abstract algebra). Gauss–Jordan elimination
, an extension of this algorithm, reduces the matrix further to diagonal form, which is also known as reduced row echelon form. Gaussian elimination alone is sufficient for many applications, and requires fewer calculations than the Gauss–Jordan version.
. Its use is illustrated in eighteen problems, with two to five equations. The first reference to the book by this title is dated to 179 CE, but parts of it were written as early as approximately 150 BCE. It was commented on by Liu Hui
in the 3rd century.
The method in Europe stems from the notes of Isaac Newton. In 1670, he wrote that all the algebra books known to him lacked a lesson for solving simultaneous equations, which Newton then supplied. Cambridge University eventually published the notes as Arithmetica Universalis in 1707 long after Newton left academic life. The notes were widely imitated, which made (what is now called) Gaussian elimination a standard lesson in algebra textbooks by the end of the 18th century. Carl Friedrich Gauss
in 1810 devised a notation for symmetric elimination that was adopted in the 19th century by professional hand computers to solve the normal equations of least-squares problems. The algorithm that is taught in high school was named for Gauss only in the 1950s as a result of confusion over the history of the subject.
equation with no solution, indicating the system has no solution. This is accomplished through the use of elementary row operations
. The second step uses back substitution to find the solution of the system above.
Stated equivalently for matrices, the first part reduces a matrix to row echelon form
using elementary row operations
while the second reduces it to reduced row echelon form, or row canonical form.
Another point of view, which turns out to be very useful to analyze the algorithm, is that Gaussian elimination computes a matrix decomposition
. The three elementary row operations used in the Gaussian elimination (multiplying rows, switching rows, and adding multiples of rows to other rows) amount to multiplying the original matrix with invertible matrices from the left. The first part of the algorithm computes an LU decomposition
, while the second part writes the original matrix as the product of a uniquely determined invertible matrix and a uniquely determined reduced row-echelon matrix.
The algorithm is as follows: eliminate x from all equations below
, and then eliminate y from all equations below 
. This will put the system into triangular form. Then, using back-substitution, each unknown can be solved for.
In the example, x is eliminated from
by adding 
to 
. x is then eliminated from 
by adding 
to 
. Formally:
The result is:
Now y is eliminated from
by adding 
to 
:
The result is:
This result is a system of linear equations in triangular form, and so the first part of the algorithm is complete.
The last part, back-substitution, consists of solving for the knowns in reverse order. It can thus be seen that
Then,
can be substituted into 
, which can then be solved to obtain
Next, z and y can be substituted into
, which can be solved to obtain
The system is solved.
Some systems cannot be reduced to triangular form, yet still have at least one valid solution: for example, if y had not occurred in
and 
after the first step above, the algorithm would have been unable to reduce the system to triangular form. However, it would still have reduced the system to echelon form. In this case, the system does not have a unique solution, as it contains at least one free variable. The solution set can then be expressed parametrically (that is, in terms of the free variables, so that if values for the free variables are chosen, a solution will be generated).
In practice, one does not usually deal with the systems in terms of equations but instead makes use of the augmented matrix
(which is also suitable for computer manipulations). For example:
Therefore, the Gaussian Elimination algorithm applied to the augmented matrix
begins with:
which, at the end of the first part (Gaussian elimination, zeros only under the leading 1) of the algorithm, looks like this:
That is, it is in row echelon form
.
At the end of the algorithm, if the Gauss–Jordan elimination(zeros under and above the leading 1) is applied:
That is, it is in reduced row echelon form, or row canonical form.
is augmented to the right of A, forming a n by 2n matrix (the block matrix
B = [A, I]). Through application of elementary row operations and the Gaussian elimination algorithm, the left block of B can be reduced to the identity matrix I, which leaves A−1 in the right block of B. Note that this is the Gauss-Jordan elimination for finding inverses.
If the algorithm is unable to reduce A to triangular form, then A is not invertible.
matrix 
. If we get "stuck" in a given column, we move to the next column. In this way, for example, some 
matrices can be transformed to a matrix that has a reduced row echelon form like
(the *'s are arbitrary entries). This echelon matrix
contains a wealth of information about 
: the rank of 
is 5 since there are 5 non-zero rows in 
; the vector space
spanned by the columns of
has a basis consisting of the first, third, fourth, seventh and ninth column of 
(the columns of the ones in 
), and the *'s tell you how the other columns of 
can be written as linear combinations of the basis columns.
This algorithm can be used on a computer for systems with thousands of equations and unknowns. However, the cost becomes prohibitive for systems with millions of equations. These large systems are generally solved using iterative method
s. Specific methods exist for systems whose coefficients follow a regular pattern (see system of linear equations).
The Gaussian elimination can be performed over any field
.
Gaussian elimination is numerically stable
for diagonally dominant or positive-definite
matrices. For general matrices, Gaussian elimination is usually considered to be stable in practice if you use partial pivoting
as described below, even though there are examples for which it is unstable.
The formal algorithm to compute
from 
follows. We write 
for the entry in row 
, column 
in matrix 
with 1 being the first index. The transformation is performed "in place", meaning that the original matrix 
is lost and successively replaced by 
.
This algorithm differs slightly from the one discussed earlier, because before eliminating a variable, it first exchanges rows to move the entry with the largest absolute value
to the "pivot position". Such "partial pivoting" improves the numerical stability of the algorithm; some variants are also in use.
Upon completion of this procedure the augmented matrix will be in row-echelon form and may be solved by back-substitution.
With the increasing popularity of multi-core processors, programmers now exploit thread-level parallel Gaussian elimination algorithms to increase the speed of computing. The shared-memory programming model (as opposed to the message exchange model) pseudocode is listed below.
// Note this code sample has been mangled (missing attr init/scope? bad gauss indentation?)...
// What is original source? Can we get valid C or else simplified pseudo code instead of this hybrid?
void parallel(int num_threads,int matrix_dimension)
int i;
for(i=0;i
pthread_attr_destroy(&attr); // Free attribute and wait for the other threads
for(i=0;i
void *gauss(int thread_id)
int i,k,j;
for(k=0;k(k%num_thread)) //interleaved-row work distribution
thread_id)
for(j=k+1;j
M[i][k]=0;}
barrier(num_thread,&mybarrier);
return NULL;
void barrier(int num_thread, barrier_t * mybarrier)
pthread_mutex_lock(&(mybarrier->barrier_mutex));
mybarrier->cur_count++;
if(mybarrier->cur_count!=num_thread)
pthread_cond_wait(&(mybarrier->barrier_cond),&(mybarrier->barrier_mutex));
else
mybarrier->cur_count=0;
pthread_cond_broadcast(&(mybarrier->barrier_cond));
pthread_mutex_unlock(&(mybarrier->barrier_mutex));
Linear algebra
Linear algebra is a branch of mathematics that studies vector spaces, also called linear spaces, along with linear functions that input one vector and output another. Such functions are called linear maps and can be represented by matrices if a basis is given. Thus matrix theory is often...
, Gaussian elimination is an algorithm
Algorithm
In mathematics and computer science, an algorithm is an effective method expressed as a finite list of well-defined instructions for calculating a function. Algorithms are used for calculation, data processing, and automated reasoning...
for solving systems of linear equations. It can also be used to find the rank
Rank (linear algebra)
The column rank of a matrix A is the maximum number of linearly independent column vectors of A. The row rank of a matrix A is the maximum number of linearly independent row vectors of A...
of a matrix
Matrix (mathematics)
In mathematics, a matrix is a rectangular array of numbers, symbols, or expressions. The individual items in a matrix are called its elements or entries. An example of a matrix with six elements isMatrices of the same size can be added or subtracted element by element...
, to calculate the determinant
Determinant
In linear algebra, the determinant is a value associated with a square matrix. It can be computed from the entries of the matrix by a specific arithmetic expression, while other ways to determine its value exist as well...
of a matrix, and to calculate the inverse of an invertible square matrix. The method is named after Carl Friedrich Gauss
Carl Friedrich Gauss
Johann Carl Friedrich Gauss was a German mathematician and scientist who contributed significantly to many fields, including number theory, statistics, analysis, differential geometry, geodesy, geophysics, electrostatics, astronomy and optics.Sometimes referred to as the Princeps mathematicorum...
, but it was not invented by him.
Elementary row operations
Elementary row operations
In mathematics, an elementary matrix is a simple matrix which differs from the identity matrix by one single elementary row operation. The elementary matrices generate the general linear group of invertible matrices...
are used to reduce a matrix to what is called triangular form (in numerical analysis) or row echelon form
Row echelon form
In linear algebra a matrix is in row echelon form if* All nonzero rows are above any rows of all zeroes, and...
(in abstract algebra). Gauss–Jordan elimination
Gauss–Jordan elimination
In linear algebra, Gauss–Jordan elimination is an algorithm for getting matrices in reduced row echelon form using elementary row operations. It is a variation of Gaussian elimination. Gaussian elimination places zeros below each pivot in the matrix, starting with the top row and working downwards....
, an extension of this algorithm, reduces the matrix further to diagonal form, which is also known as reduced row echelon form. Gaussian elimination alone is sufficient for many applications, and requires fewer calculations than the Gauss–Jordan version.
History
The method of Gaussian elimination appears in Chapter Eight, Rectangular Arrays, of the important Chinese mathematical text Jiuzhang suanshu or The Nine Chapters on the Mathematical ArtThe Nine Chapters on the Mathematical Art
The Nine Chapters on the Mathematical Art is a Chinese mathematics book, composed by several generations of scholars from the 10th–2nd century BCE, its latest stage being from the 1st century CE...
. Its use is illustrated in eighteen problems, with two to five equations. The first reference to the book by this title is dated to 179 CE, but parts of it were written as early as approximately 150 BCE. It was commented on by Liu Hui
Liu Hui
Liu Hui was a mathematician of the state of Cao Wei during the Three Kingdoms period of Chinese history. In 263, he edited and published a book with solutions to mathematical problems presented in the famous Chinese book of mathematic known as The Nine Chapters on the Mathematical Art .He was a...
in the 3rd century.
The method in Europe stems from the notes of Isaac Newton. In 1670, he wrote that all the algebra books known to him lacked a lesson for solving simultaneous equations, which Newton then supplied. Cambridge University eventually published the notes as Arithmetica Universalis in 1707 long after Newton left academic life. The notes were widely imitated, which made (what is now called) Gaussian elimination a standard lesson in algebra textbooks by the end of the 18th century. Carl Friedrich Gauss
Carl Friedrich Gauss
Johann Carl Friedrich Gauss was a German mathematician and scientist who contributed significantly to many fields, including number theory, statistics, analysis, differential geometry, geodesy, geophysics, electrostatics, astronomy and optics.Sometimes referred to as the Princeps mathematicorum...
in 1810 devised a notation for symmetric elimination that was adopted in the 19th century by professional hand computers to solve the normal equations of least-squares problems. The algorithm that is taught in high school was named for Gauss only in the 1950s as a result of confusion over the history of the subject.
Algorithm overview
The process of Gaussian elimination has two parts. The first part (Forward Elimination) reduces a given system to either triangular or echelon form, or results in a degenerateDegeneracy (mathematics)
In mathematics, a degenerate case is a limiting case in which a class of object changes its nature so as to belong to another, usually simpler, class....
equation with no solution, indicating the system has no solution. This is accomplished through the use of elementary row operations
Elementary row operations
In mathematics, an elementary matrix is a simple matrix which differs from the identity matrix by one single elementary row operation. The elementary matrices generate the general linear group of invertible matrices...
. The second step uses back substitution to find the solution of the system above.
Stated equivalently for matrices, the first part reduces a matrix to row echelon form
Row echelon form
In linear algebra a matrix is in row echelon form if* All nonzero rows are above any rows of all zeroes, and...
using elementary row operations
Elementary row operations
In mathematics, an elementary matrix is a simple matrix which differs from the identity matrix by one single elementary row operation. The elementary matrices generate the general linear group of invertible matrices...
while the second reduces it to reduced row echelon form, or row canonical form.
Another point of view, which turns out to be very useful to analyze the algorithm, is that Gaussian elimination computes a matrix decomposition
Matrix decomposition
In the mathematical discipline of linear algebra, a matrix decomposition is a factorization of a matrix into some canonical form. There are many different matrix decompositions; each finds use among a particular class of problems.- Example :...
. The three elementary row operations used in the Gaussian elimination (multiplying rows, switching rows, and adding multiples of rows to other rows) amount to multiplying the original matrix with invertible matrices from the left. The first part of the algorithm computes an LU decomposition
LU decomposition
In linear algebra, LU decomposition is a matrix decomposition which writes a matrix as the product of a lower triangular matrix and an upper triangular matrix. The product sometimes includes a permutation matrix as well. This decomposition is used in numerical analysis to solve systems of linear...
, while the second part writes the original matrix as the product of a uniquely determined invertible matrix and a uniquely determined reduced row-echelon matrix.
Example
Suppose the goal is to find and describe the solution(s), if any, of the following system of linear equations:The algorithm is as follows: eliminate x from all equations below
, and then eliminate y from all equations below . This will put the system into triangular form. Then, using back-substitution, each unknown can be solved for.In the example, x is eliminated from
by adding to . x is then eliminated from by adding to . Formally:The result is:
Now y is eliminated from
by adding to :The result is:
This result is a system of linear equations in triangular form, and so the first part of the algorithm is complete.
The last part, back-substitution, consists of solving for the knowns in reverse order. It can thus be seen that
Then,
can be substituted into , which can then be solved to obtain
Next, z and y can be substituted into
, which can be solved to obtain
The system is solved.
Some systems cannot be reduced to triangular form, yet still have at least one valid solution: for example, if y had not occurred in
and after the first step above, the algorithm would have been unable to reduce the system to triangular form. However, it would still have reduced the system to echelon form. In this case, the system does not have a unique solution, as it contains at least one free variable. The solution set can then be expressed parametrically (that is, in terms of the free variables, so that if values for the free variables are chosen, a solution will be generated).In practice, one does not usually deal with the systems in terms of equations but instead makes use of the augmented matrix
Augmented matrix
In linear algebra, an augmented matrix is a matrix obtained by appending the columns of two given matrices, usually for the purpose of performing the same elementary row operations on each of the given matrices.Given the matrices A and B, where:A =...
(which is also suitable for computer manipulations). For example:
Therefore, the Gaussian Elimination algorithm applied to the augmented matrix
Augmented matrix
In linear algebra, an augmented matrix is a matrix obtained by appending the columns of two given matrices, usually for the purpose of performing the same elementary row operations on each of the given matrices.Given the matrices A and B, where:A =...
begins with:
which, at the end of the first part (Gaussian elimination, zeros only under the leading 1) of the algorithm, looks like this:
That is, it is in row echelon form
Row echelon form
In linear algebra a matrix is in row echelon form if* All nonzero rows are above any rows of all zeroes, and...
.
At the end of the algorithm, if the Gauss–Jordan elimination(zeros under and above the leading 1) is applied:
That is, it is in reduced row echelon form, or row canonical form.
Finding the inverse of a matrix
Suppose A is a n by n Invertible matrix. Then the inverse of A can be found as follows. The n by n identity matrixIdentity matrix
In linear algebra, the identity matrix or unit matrix of size n is the n×n square matrix with ones on the main diagonal and zeros elsewhere. It is denoted by In, or simply by I if the size is immaterial or can be trivially determined by the context...
is augmented to the right of A, forming a n by 2n matrix (the block matrix
Block matrix
In the mathematical discipline of matrix theory, a block matrix or a partitioned matrix is a matrix broken into sections called blocks. Looking at it another way, the matrix is written in terms of smaller matrices. We group the rows and columns into adjacent 'bunches'. A partition is the rectangle...
B = [A, I]). Through application of elementary row operations and the Gaussian elimination algorithm, the left block of B can be reduced to the identity matrix I, which leaves A−1 in the right block of B. Note that this is the Gauss-Jordan elimination for finding inverses.
If the algorithm is unable to reduce A to triangular form, then A is not invertible.
General algorithm to compute ranks and bases
The Gaussian elimination algorithm can be applied to any matrix . If we get "stuck" in a given column, we move to the next column. In this way, for example, some matrices can be transformed to a matrix that has a reduced row echelon form like(the *'s are arbitrary entries). This echelon matrix
contains a wealth of information about : the rank of is 5 since there are 5 non-zero rows in ; the vector spaceVector space
A vector space is a mathematical structure formed by a collection of vectors: objects that may be added together and multiplied by numbers, called scalars in this context. Scalars are often taken to be real numbers, but one may also consider vector spaces with scalar multiplication by complex...
spanned by the columns of
has a basis consisting of the first, third, fourth, seventh and ninth column of (the columns of the ones in ), and the *'s tell you how the other columns of can be written as linear combinations of the basis columns.Analysis
Gaussian elimination to solve a system of n equations for n unknowns requires n(n+1) / 2 divisions, (2n3 + 3n2 − 5n)/6 multiplications, and (2n3 + 3n2 − 5n)/6 subtractions, for a total of approximately 2n3 / 3 operations. Thus it has arithmetic complexity of O(n3). However, the intermediate entries can grow exponentially large, so it has exponential bit complexity.This algorithm can be used on a computer for systems with thousands of equations and unknowns. However, the cost becomes prohibitive for systems with millions of equations. These large systems are generally solved using iterative method
Iterative method
In computational mathematics, an iterative method is a mathematical procedure that generates a sequence of improving approximate solutions for a class of problems. A specific implementation of an iterative method, including the termination criteria, is an algorithm of the iterative method...
s. Specific methods exist for systems whose coefficients follow a regular pattern (see system of linear equations).
The Gaussian elimination can be performed over any field
Field (mathematics)
In abstract algebra, a field is a commutative ring whose nonzero elements form a group under multiplication. As such it is an algebraic structure with notions of addition, subtraction, multiplication, and division, satisfying certain axioms...
.
Gaussian elimination is numerically stable
Numerical stability
In the mathematical subfield of numerical analysis, numerical stability is a desirable property of numerical algorithms. The precise definition of stability depends on the context, but it is related to the accuracy of the algorithm....
for diagonally dominant or positive-definite
Positive-definite matrix
In linear algebra, a positive-definite matrix is a matrix that in many ways is analogous to a positive real number. The notion is closely related to a positive-definite symmetric bilinear form ....
matrices. For general matrices, Gaussian elimination is usually considered to be stable in practice if you use partial pivoting
Pivot element
The pivot or pivot element is the element of a matrix, an array, or some other kind of finite set, which is selected first by an algorithm , to do certain calculations...
as described below, even though there are examples for which it is unstable.
Higher order tensors
Gaussian elimination does not generalize in any simple way to higher order tensors (matrices are order 2 tensors); even computing the rank of a tensor of order greater than 2 is a difficult problem.Pseudocode
As explained above, Gaussian elimination writes a given m × n matrix A uniquely as a product of an invertible m × m matrix S and a row-echelon matrix T. Here, S is the product of the matrices corresponding to the row operations performed.The formal algorithm to compute
from follows. We write for the entry in row , column in matrix with 1 being the first index. The transformation is performed "in place", meaning that the original matrix is lost and successively replaced by .
for k = 1 ... m:
Find pivot for column k:
i_max := argmax (i = k ... m, abs(A[i, k]))
if A[i_max, k] = 0
error "Matrix is singular!"
swap rows(k, i_max)
Do for all rows below pivot:
for i = k + 1 ... m:
Do for all remaining elements in current row:
for j = k + 1 ... n:
A[i, j] := A[i, j] - A[k, j] * (A[i, k] / A[k, k])
Fill lower triangular matrix with zeros:
A[i, k] := 0
This algorithm differs slightly from the one discussed earlier, because before eliminating a variable, it first exchanges rows to move the entry with the largest absolute value
Absolute value
In mathematics, the absolute value |a| of a real number a is the numerical value of a without regard to its sign. So, for example, the absolute value of 3 is 3, and the absolute value of -3 is also 3...
to the "pivot position". Such "partial pivoting" improves the numerical stability of the algorithm; some variants are also in use.
Upon completion of this procedure the augmented matrix will be in row-echelon form and may be solved by back-substitution.
With the increasing popularity of multi-core processors, programmers now exploit thread-level parallel Gaussian elimination algorithms to increase the speed of computing. The shared-memory programming model (as opposed to the message exchange model) pseudocode is listed below.
// Note this code sample has been mangled (missing attr init/scope? bad gauss indentation?)...
// What is original source? Can we get valid C or else simplified pseudo code instead of this hybrid?
void parallel(int num_threads,int matrix_dimension)
int i;
for(i=0;i
pthread_attr_destroy(&attr); // Free attribute and wait for the other threads
for(i=0;i
void *gauss(int thread_id)
int i,k,j;
for(k=0;k
(k%num_thread)) //interleaved-row work distribution
for(j=k+1;j
M[k][k]=1;
barrier(num_thread,&mybarrier); //wait for other thread finishing this round
for(i=k+1;i
thread_id)for(j=k+1;j
M[i][k]=0;}
barrier(num_thread,&mybarrier);
return NULL;
void barrier(int num_thread, barrier_t * mybarrier)
pthread_mutex_lock(&(mybarrier->barrier_mutex));
mybarrier->cur_count++;
if(mybarrier->cur_count!=num_thread)
pthread_cond_wait(&(mybarrier->barrier_cond),&(mybarrier->barrier_mutex));
else
mybarrier->cur_count=0;
pthread_cond_broadcast(&(mybarrier->barrier_cond));
pthread_mutex_unlock(&(mybarrier->barrier_mutex));
See also
- LU decompositionLU decompositionIn linear algebra, LU decomposition is a matrix decomposition which writes a matrix as the product of a lower triangular matrix and an upper triangular matrix. The product sometimes includes a permutation matrix as well. This decomposition is used in numerical analysis to solve systems of linear...
can be viewed as a matrix form of Gaussian elimination - Frontal solverFrontal solverA frontal solver, due to Irons is an approach to solving sparse linear systems which is used extensively in finite element analysis. It is a variant of Gauss elimination that automatically avoids a large number of operations involving zero terms....
, an improvement of Gaussian elimination for the case of sparse systems - Gauss–Jordan eliminationGauss–Jordan eliminationIn linear algebra, Gauss–Jordan elimination is an algorithm for getting matrices in reduced row echelon form using elementary row operations. It is a variation of Gaussian elimination. Gaussian elimination places zeros below each pivot in the matrix, starting with the top row and working downwards....
, a procedure that further reduces the coefficient matrix to the point where back-substitution is no longer necessary - Grassmann's algorithm is a numerically stable variant of Gaussian elimination for systems over the real or complex numbers. The algorithm avoids subtractions so is less sensitive to errors caused by subtraction of two roughly equal numbers.
External links
- WebApp descriptively solving systems of linear equations with Gaussian Elimination
- A program that performs Gaussian elimination similarly to a human working on paper Exact solutions to systems with rational coefficients.
- Gaussian elimination www.math-linux.com.
- Gaussian elimination as java applet at some local site. Only takes natural coefficients.
- Gaussian elimination at Holistic Numerical Methods Institute
- LinearEquations.c Gaussian elimination implemented using C language
- Gauss–Jordan elimination Step by step solution of 3 equations with 3 unknowns using the All-Integer Echelon Method provides a very clear, elementary presentation of the method of row reduction.