Gaussian elimination
The Gaussian elimination method or simply Gaussian process (after Carl Friedrich Gauss ) is an algorithm of the mathematical sciences of linear algebra and numerical analysis. It is an important method for solving systems of linear equations and based on the fact that although elementary transformations change the system of equations, but get the solution. This allows to make any clearly solvable system of equations step shape of the light determines the solution by successive elimination of the unknowns or the solution set can be read.
The number of required operations is in a matrix of the order. In its basic form, the algorithm is sensitive to rounding errors, but with small modifications ( Pivoting ) it provides for general systems of linear equations, the standard solution method represents and is part of all major libraries for numerical linear algebra such as NAG, IMSL and LAPACK.
- 3.1 forward substitution
- 3.2 backward substitution
- 3.3 algorithm in pseudocode
- 3.4 Incomplete decompositions
- 4.1 computational complexity and memory requirements
- 4.2 accuracy 4.2.1 iterative refinement
- 5.1 statements on the solvability of the linear system
- 5.2 determinant
- 5.3 Calculation of the inverse
Explanation
A linear system of three equations and three unknowns, and the right side of the form:
The algorithm for the calculation of the variable, and can be divided into two stages:
In the first step, the equation system is placed on step-like shape. Stepped shape means that each line, at least one variable occurs less, that at least one variable is eliminated. In the above equations would be, and eliminate, in the third line is then only the variable:
To achieve the levels form elementary row operations are used, with the help of which the system of equations is transformed into a new, but which has the same solution set. Sufficient are two types of elementary row operations:
The method then consists, starting in the first column to do with transformations of the first kind given by skillfully adding the first row, all entries up to the first zero. This is then placed in the thus modified second column continues with time multiple of the second line are added to the following row and so on. This step works only if the diagonal element of the current column is not zero. In such a case, the second type of the line forming is necessary, since a non-zero entry on the diagonal can be produced by a row interchange. With the help of these two kinds of transformations it is possible to bring any linear system of equations on step shape.
Another type of elementary transformation the swapping of columns. It is not required to implement the algorithm, but sometimes employed in the computer program, for stability reasons. The position of the variable is changed in the system of equations. In computing per capita is sometimes still multiplying a row by a number useful to avoid about to complicated fractures. This causes additional computational effort, and is therefore not an option in computer programs and further changes the determinant of the coefficient matrix, which brings theoretical disadvantages.
In the second step of the process, the reverse insertion can be based on the last line, in which only a variable appears, the variable is calculated and used in the overlying line.
An alternative to this is the Gauss-Jordan algorithm, wherein not only the lower portions are eliminated, but also the top, thus creating a diagonal shape in which then eliminates the above-mentioned second step.
Example
For clarity, the coefficients of the linear system of equations is written in the with extended coefficient matrix:
Now is reshaped so that and are zero, by adding suitable multiples of the first equation to the second and third equations. The corresponding multiplier is obtained by adding the element to be eliminated ( first) divided by the pivot element (here and ). Since the two elements and zero to be, the two multipliers are respectively multiplied.
Skip to second row so the times and the third line shows the times the first row is added. So that zero is a multiple of the second line is added to the third row, in this case the times:
If the number divided by to calculate the multiplier (here for the first two lines of the figure, the third time the number) is zero, this row interchanged with a lower down. The last line means
This equation is easy to solve and supplies. Thus results for the second line
And on. This means that all variables are calculated:
Control by the row sum
The transformation may be monitored by calculating the row sum.
Here, the sum of all elements of each row was added in the last column. For the first row is the row-sum. Since no transformations are performed on the first line, its row sum does not change. In the first transformation of this equation system to the second line, the times will be the first added. If one does, it also applies to the row total. This result is the row sum of the formed second line. To check the bills so you can perform the transformations at the row total. Are all invoices correctly, the row sum of the transformed line must arise.
Pivoting
The Gaussian elimination method is feasible in general not without row permutations. Is replaced by the above example, the algorithm without row interchange can not start. To remedy this, one chooses an element of the first column of the coefficient matrix, called the pivot element, which is equal to 0.
Then we interchange the first row with the pivot row:
For the calculation by hand, it is helpful to have a 1 or minus one to choose as pivot element, so that in the further course of the process there are no breaks. For the calculation with the aid of a computer, it makes sense to choose the amount largest element in order to obtain a stable algorithm as possible. If we choose the pivot element in the current column, it is called Spaltenpivotisierung. Alternatively, you can select the pivot in the current line.
In this case the columns are exchanged accordingly.
In the backward substitution is to be observed that the variables have changed their position in the system of equations. If we choose as pivot the element largest amount of total residual matrix, then one speaks of complete Pivoting or full pivoting. For both row and Spaltenvertauschungen are necessary in general.
Pivoting is feasible without significant additional effort, if not the entries of the matrix and the right-hand side are reversed, but the permutations are stored in an index vector.
LU decomposition
The delivery of solving a quadratic equation with a unique solution system as a computer program, it makes sense (also known as LU decomposition or triangular factorization ) the Gaußalgorithmus as LU decomposition to interpret. This is a regular separation of the matrix into the product of a lower left triangular matrix (left or Engl. "Lower " ), and an upper-right triangular matrix (to the right, also referred to by Engl. "Upper "). The following example illustrates this:
It has the above-mentioned step-like shape, and the matrix used to store the required forming steps corresponding multiplications Frobeniusmatrizen. This shows the existence of the decomposition. To achieve uniqueness, the diagonal elements of the matrix are defined as 1. Saving the transformation steps has the advantage that different "right side " the system of equations can be solved efficiently by forward and backward substitution.
The row permutations generally required can be described by a permutation matrix.
Forward substitution
In the forward substitution to compute a solution of the linear system, or in the case of bill with Pivoting. This is by the equation with the solution of the equation system in the original relationship.
Written out, the equation system has the following form:
Then the following formula holds for the components:
Starting with succession can be calculated, respectively by used already known.
Backward substitution
The backward substitution to calculate the solution of the original system of equations by dissolving similar to the forward insertion. The difference being that one starts with, and then successively calculates the values of. The corresponding formula is
Algorithm in pseudocode
The following algorithm performs a LU decomposition of the matrix A without Pivoting by simultaneously L and R generated on the basis of A:
Input: matrix A / / Initialization R: = A L: = E_n / / N-1 iterations for i: = 1 to n-1 / / Rows of the residual matrix are run through for k: = i 1 to n / / Calculation of L L ( k, i): = R ( k, i ) / R (i, i) / / Note: before checking for null values necessary / / Columns of the residual matrix are run through for j: = i to n / / Calculation of R R ( k, j ) = R ( k, j) - L ( k, i) * R (i, j) Output: matrix L, the matrix R Final describe L and R just the factorization of A in the desired form, ie A = LR for a left lower triangular matrix L and a right upper triangular matrix R.
Alternatively, ( from possible interest of memory efficiency), a simultaneous development of the L and R in the A directly possible to be described by the following algorithm:
Input: matrix A / / N-1 iterations for i: = 1 to n-1 / / Rows of the residual matrix are run through for k: = i 1 to n / / Calculation of L A ( k, i): = A ( k, i ) / A (i, i) / / Note: before checking for null values necessary / / Columns of the residual matrix are run through for j: = i 1 to n / / Calculation of R A (k, j ) = A ( k, j) - A ( k, i) * a ( i, j) Output: Matrix A ( in modified form ) Incomplete decompositions
The LU decomposition has the drawback that it is often fully occupied even with sparse matrices. Then instead of all entries, only those calculated in a given occupation pattern, it is called an incomplete LU factorization. This provides a convenient approximation to the matrix and can thus be used as a preconditioner in the iterative solution of linear systems of equations. In the case of symmetric positive definite matrices is called an incomplete Cholesky decomposition.
Features of the method
Computational complexity and memory requirements
The number of arithmetic operations for the LU decomposition is about at a matrix. The cost of the forward and backward substitution is square () and therefore negligible overall. However, since the cost increases with the cube of the dimension of the matrix, one can estimate the calculation time for a different array on the basis of computation time for a matrix. What is needed therefore the algorithm on a specific computer for a matrix of dimension about 10 seconds, so he needs for a matrix of dimension on the same machine about seconds about 3 hours. Thus, the Gaussian elimination method is a fast direct method for solving linear systems of equations, for a QR decomposition requires at least twice as many arithmetic operations. Nevertheless, the algorithm should be used only for small systems of equations to medium dimension (up to about ). For matrices of higher dimension, however, the computing time increases fast and iterative methods are often better. This approach to the solution gradually and need in each step for a fully staffed matrix arithmetic operations. The convergence rate of such methods is highly dependent on the properties of the matrix and can be the computing time required concrete is difficult to predict.
The calculation can be carried out in the memory of the matrix, so that in addition to the storage of even no additional memory requirement arises. For a fully staffed matrix of dimension would have to save a million coefficients. This corresponds to the IEEE 754 format double to about 8 megabytes. In iterative methods that work with matrix-vector multiplications, but an explicit storage of itself may not be necessary, so that these processes are possibly preferable.
For special cases, time and storage can be significantly reduced by using specific properties of the matrix and its LU decomposition can be exploited. As required for the Cholesky decomposition of symmetric positive definite matrices, only half the amount of arithmetic operations and memory. Another example are banded matrices with fixed bandwidth, since the LU decomposition obtains the band structure and thus reduces the effort on. For a few special sparse matrices, it is possible to exploit the cast structure, so that the LU decomposition is also sparse. Both associated with a reduced memory requirement.
Accuracy
Thus, the calculation of x is sufficiently accurate, may on the one hand the condition of the matrix is not too bad and the machine precision used is not too low. On the other hand it requires a solution process, which is sufficiently stable. A good algorithm is thus characterized by a high stability.
Generally, the process without Pivoting is unstable. Therefore, most Spaltenpivotisierung is used for solution. Thus, the method for most matrices is stable feasible, especially as it became clear by the work of James H. Wilkinson after the Second World War. It can, however, specify matrices for which the stability constant grows exponentially with the dimension of the matrix. With complete Pivoting the stability can be improved further, but then also increases the cost of the pivot search, so this is rarely used. Generally have better stability QR factorization, but these are also expensive to compute.
In strictly diagonally dominant or positive definite matrices (see also Cholesky decomposition ), the Gaussian method can be performed stably without Pivoting, occur so no zeros on the diagonal on.
Iterative refinement
A practical approach to compensate for these inaccuracies of calculation consists of an iterative refinement by splitting procedure, there stands on the LU decomposition is a good approximation of the matrix A are available, which is easily invertible. For this, start with the computed solution and calculated in each step the residual
After that calculated using the LU decomposition, the solution of the equation system
And sets
As it usually goes only to small corrections, often rich few iteration steps. In general, a higher accuracy is, however, necessary for the calculation of the residual. Reaches the iterative refinement is not sufficient to attain the desired accuracy, remains only the choice of another procedure or a transformation of the problem to obtain a more favorable matrix, such as a smaller condition.
The iterative refinement is used for example in the LAPACK routine DSGESV. In this routine, the LU decomposition is determined in single-precision and double precision achieved by iterative refinement of the solution with exactly twice the calculated residual.
The Gaussian process as a theoretical tool
The Gaussian process is in addition to its importance for the numerical treatment of linear systems of equations with a unique solution also an important tool in theoretical linear algebra.
Statements about the solvability of the linear system
A system of linear equations may have one, several, or no solution. When using the full Pivoting Gauss method brings each matrix to echelon form triangle in which all entries below a certain row are zero and the diagonal appear no zero entries. The rank of the matrix is then the number of non- zero rows of the matrix. The solubility is then obtained from the interaction with the right side: Belong to the zero -line non-zero entries of the right-hand side, the system of equations is unsolvable, otherwise solvable, the dimension of the solution set of the number of unknowns equal to minus the rank.
Since the second equation is a multiple of the first equation, the equation system has infinite number of solutions. At the elimination of X in the second equation to disappear completely, leaving only the first equation. Solving this for x, you can specify the amount of solution as a function of y:
Determinant
Further, the Gaussian method provides a possibility for calculating the determinant of a matrix. Since the elementary row transformations have determinant 1, up to row permutations whose determinant is -1 ( but this only changes the sign and can easily correct therefore ), the resulting upper triangular matrix has the same determinant as the original matrix, but can much easier be calculated: it is the product of the diagonal elements.
Calculation of the inverse
Another possibility of using the Gaussian method is to calculate the inverse of the matrix. To this end, the algorithm is applied to a right- extended by a unit matrix scheme and continued after said first phase, to the left, a unit matrix has been reached. In the right part is the inverse matrix. This method is numerically not recommended and the explicit calculation of the inverse can usually be circumvented.
History
Already in the Chinese mathematics book Jiu Zhang SuanShu (Eng. Nine Books arithmetic technique ), which was written between AD 200 and before 100, there is an example, but a clear demonstration of the algorithm based on the solution of a system with three unknowns. 263 Liu Hui published a comprehensive review of the book, which was received subsequently in the text corpus. The Jiuzhang SuanShu was to the 16th century a significant source of mathematical education in China and surrounding countries.
In Europe, a process was only in 1759 by Joseph -Louis Lagrange published, containing the basic elements. Carl Friedrich Gauss employed as part of its development and application of the method of least squares with linear systems of equations, the normal equations occurring there. His first publication on the subject dates from 1810 ( Disquisitio de elementis ellipticis Palladis ), but he mentioned in 1798 in his diaries cryptically that he had solved the problem of elimination. It is certain that he used the method of calculating the orbit of the asteroid Pallas 1803-1809. In the 1820s, he described the first time something like an LU decomposition. The elimination procedure was used in the following years mainly in geodesy (see in Gauss ' services ), and then the second namesake of the Gauss-Jordan method is not about the mathematician Camille Jordan, but the surveyor William Jordan.
During and after the Second World War, the investigation of numerical methods gained in importance and the Gauss method has now been applied increasingly to problems regardless of the method of least squares. John von Neumann and Alan Turing defined the LU decomposition in the usual form today and examined the phenomenon of rounding errors. These issues have been satisfactorily resolved until the 1960s by James Hardy Wilkinson, who showed that the method with Pivoting is backward stable.