Often optimization procedures result in a set of N linear equations with N unknown variable values. Gaussian elimination is one method for determining the variable values.
There are diverse ways to present a system of linear equations, but these are simply notational variations. One is a classic reveal of each equation. Equation set (C1) is an example with three unknowns, x, y, and z. The coefficients a, b, and c, and the right‐hand side elements, r, all have known values. Note that the unknowns x, y, and z each appear linearly (to the first power and independent of other unknowns) in the equation set. Remove the third equation and the cz elements, and it is a system of two equations and two unknowns. In general it could be a system of N equations linear in N unknowns. If a variable does not appear in a particular equation, insert it with a coefficient value of zero:
Note that the structure of the abovementioned equation set has all of the unknown variables and coefficients on the left‐hand side, in the same order, and the constant on the right‐hand side. This makes it convenient to represent in matrix–vector notation:
Although now there is only one equal sign, for the LHS matrix–vector product (a vector) to equal the RHS vector, each element on the vector must be identical, indicating the same three equations. Note that the matrix is square, the number of rows represents the number of equations, and the number of columns represents the coefficients on each unknown.
A shorthand matrix–vector representation is
There is no need to have separate symbol for each unknown. Instead of x, y, and z, they can be x1, x2, and x3. Similarly the matrix of coefficients could use the symbol mi,j for each element with i representing the row and j representing the column (counting from the upper left element):
Gaussian elimination is a two‐part procedure that solves for the x‐values of equation set (C3). The result of the first part of the procedure is to convert the system to one in which the elements of the matrix are unity along the main diagonal and zero in the lower triangle.
With this structure the value of x3 is easily obtained from the third equation. . With x3 evaluated, the second to bottom equation permits a direct calculation for x2. This is the second part of the procedure.
The first part of the procedure is to convert the equation set into the desired form.
This is the desired structure that permits a simple step‐by‐step solution.
Note: If there were N equations, this first part would have N stages.
Solve for each variable in reverse order:
That was easy! However, there could be a problem if a leading coefficient is zero. Then you have a divide by zero in the normalization (A‐steps) of each stage. Even if the original coefficient value is not zero, this could happen because it becomes zero after the stages. The solution is to resort the columns to make leading coefficient not zero. Better yet, prior to each A‐step, sort the columns so that coefficient with the maximum absolute value is the lead, then this minimizes the effects of truncation error in calculations. This is termed pivoting.
For example, an original equation set might be
The A‐step would direct that all elements in the first equation are to be divided by the x‐coefficient value of “+3.” But the z‐coefficient value of “–10” has a larger absolute value. Accordingly, rearrange the equation set so that z is the first unknown variable in the list. This does not change the solution to the equations:
Note: The elements in the unknown vector have shifted order from (x, y, z) to (z, y, x). In matrix–vector notation,
Switching coefficients to make the leading one the maximum requires switching all column elements in the coefficients matrix and row elements in solution column vector.
For simplicity, there is no need to normalize every successive row, just the one you’re leaving, and then subtract that row times the leading coefficient:
If a divide by zero occurs, it is because the solution is indeterminate (if two or more equations are linearly dependent). So, if a leading coefficient is zero (or equivalently zero numerically, for instance, its absolute value is less than 0.00000001), then send an “indeterminate” message and stop execution.
Since the rows operated on are at the stage number or higher, this procedure can be expressed in pseudocode.
Since pivoting scrambles the solution position in the solution vector, the answer needs to be unscrambled. In the example, in Equation (C19) the solution {x, y, z} became reordered as {z, y, x}. To keep track of the order of elements in the solution vector, I use a place vector, in which place holds the original element index.
For example, place is initialized as
After the pivot, x and z are switched so the place vector becomes
There is no need to have an M matrix and an M′ matrix and an M″ matrix. Since each matrix and RHS vector operations update elements, just update the m element.
For example,
The “:=” acknowledges that this is a computer assignment statement, not the algebraic “equals.”
This column switching in the coefficient matrix to make the leading coefficient of the top working row be the coefficient with largest magnitude (absolute value) is called partial pivoting. Full pivoting also switches working rows so that the divide by value in the upper left of the working part of the matrix is the largest magnitude of all remaining in the working matrix.
' Gauss Elimination with Pivoting
' R. Russell Rhinehart
' 29 June 2009, 23 Aug 2012
' Gauss Elimination of a system of N (up to 10) linear equations with pivoting on lead coefficient
Option Explicit
Dim H(10, 10) As Double ' Coefficient Matrix - Hessian in Newton-Raphson
Dim V(10) As Double ' RHS vector - Jacobean or Gradient in Newton-Raphson
Dim X(10) As Double ' Vector of Unknowns - Change in Decision Variables in Optimization
Dim XS(10) As Double ' Scrambled Vector of Unkowns
Dim Place(10) As Integer ' Position of solution vector
Dim TC As Double ' Temporary holder for coefficient value for switching
Dim TI As Integer ' Temporary value of column index when switching or unscrambling
Dim Nvar As Integer ' Number of variables, length of solution vector
Dim i As Integer ' index counter for row
Dim k As Integer ' index counter for row
Dim j As Integer ' index counter for column
Dim MaxCoefficient As Double ' value for largest coefficient in the row
Dim Index As Integer ' k-value for largest coefficient
Dim Coefficient As Double ' value of rearragned equation lead coefficient
Dim sum As Double ' Sum of values in reconstructing the solution
Sub Gauss_Elimination_w_Pivot()
Nvar = Cells(17, 13)
' Initialize sequence of solution vector elements
For i = 1 To Nvar
Place(i) = i
Next i
' Read Data
For i = 1 To Nvar
For j = 1 To Nvar
H(i, j) = Cells(i + 5, j + 1)
Next j
Next i
For i = 1 To Nvar
V(i) = Cells(i + 5, 13)
Next i
' Perform row operations to get upper-right diagonal matrix
For k = 1 To Nvar - 1
' Search for largest coefficient in row k
MaxCoefficient = Abs(H(k, k))
Index = k
For j = k + 1 To Nvar
If Abs(H(k, j)) > MaxCoefficient Then
MaxCoefficient = Abs(H(k, j))
Index = j
End If
Next j
' switch Index and k columns
For i = 1 To Nvar
TC = H(i, Index)
H(i, Index) = H(i, k)
H(i, k) = TC
Next i
TI = Place(Index)
Place(Index) = Place(k)
Place(k) = TI
' normalize lead coefficient of row k to unity - make H(k,k) = 1
Coefficient = H(k, k)
If Abs(Coefficient) < 0.0000000001 Then 'test for linearly dependent equations.
For i = 1 To Nvar
Cells(i + 5, 17) = "Indeterminate"
Next i
End
End If
For j = 1 To Nvar
H(k, j) = H(k, j) / Coefficient
Next j
V(k) = V(k) / Coefficient
' subtract row-k-times-lead-coefficient in each remaining row (i>k) - make M(i,k) = 0
For i = k + 1 To Nvar
Coefficient = H(i, k)
For j = 1 To Nvar
H(i, j) = H(i, j) - H(k, j) * Coefficient
Next j
V(i) = V(i) - V(k) * Coefficient
Next i
Next k
' normalize last row - row i=Nvar
Coefficient = H(Nvar, Nvar)
If Abs(Coefficient) < 0.0000000001 Then
For i = 1 To Nvar
Cells(i + 5, 17) = "Indeterminate"
Next i
End
End If
For j = 1 To Nvar
H(Nvar, j) = H(Nvar, j) / Coefficient
Next j
V(Nvar) = V(Nvar) / Coefficient
' solve for scrambled solution values
XS(Nvar) = V(Nvar)
For j = Nvar - 1 To 1 Step -1
sum = 0
For k = j + 1 To Nvar
sum = sum + XS(k) * H(j, k)
Next k
XS(j) = V(j) - sum
Next j
' Unscramble solution
For i = 1 To Nvar
TI = Place(i)
X(TI) = XS(i)
Next i
' Display answer
For i = 1 To Nvar
Cells(i + 5, 17) = X(i)
Next i
' check if solution vector returns RHS vector
' Re-read original matrix coefficients
For i = 1 To Nvar
For j = 1 To Nvar
H(i, j) = Cells(i + 5, j + 1)
Next j
Next i
' Calculate and display RHS
For i = 1 To Nvar
sum = 0
For j = 1 To Nvar
sum = sum + X(j) * H(i, j)
Next j
Cells(i + 5, 14) = sum
Next i
End Sub
Sub GenerateRHS()
Nvar = Cells(17, 13)
For i = 1 To Nvar
For j = 1 To Nvar
H(i, j) = Cells(i + 5, j + 1)
Next j
X(i) = Cells(i + 5, 18)
Next i
' Calculate and display RHS
For i = 1 To Nvar
sum = 0
For j = 1 To Nvar
sum = sum + X(j) * H(i, j)
Next j
Cells(i + 5, 13) = sum
Next i
End Sub