Introduction
Gaussian elimination is a technique that is often used to solve a system of linear equations, as it is a very stable method of solving them. There are many examples available
around the web that shows you how to solve them, but they are seldom explained very well, why they work and what the potential problem is, referring especially to the potential
roundoff errors. The key is to understand how the calculation of Gaussian elimination and the inverse of the matrix works, and when you do, it is rather trivial to solve the equations.
Both these methods of calculation will give you a general method in the toolbox for many important applications, as for example Spline calculations, Polynomial fits, etc.
Background
A set of linear equations can be solved using Gaussian elimination, where we need to equal the number of variables and the number of equations,
to be able to solve them, and can look like this:
However, to solve this list of equations, the approach I'm going to take is to write the equations in a matrix form:
To solve the list of equations for the unknowns we want to calculate it to get the unity matrix on the left side and the solution on the right side, as shown below:
The matrix on the most left side, with the ones forming a diagonal surrounded by zeroes, is called an
identity matrix. If you recall earlier math classes where this is discussed,
you would also notice that you have the same amount of rows and columns in the matrix, as you need X number of equations to solve for X unknowns. There is also the notion that
when the two matrices on the left hand side are multiplied together we get the following values, and the solution to the equation as:
From this information we can immediately get that the values in the diagonal, that came from the three original equations, cannot be zero. In our case this means that in the first
equation x cannot be zero, in equation two, y cannot be zero, and in equation three z cannot be zero. If they are,
we would have to rearrange the equations so that this is not the case,
for example changing equation one with equation two, etc., without violating the first condition that the diagonal must not be zero.
The equations in the example given here have no need to be rearranged as they do not violate the condition, so we begin to construct the unity matrix by
dividing the first
equation by 2, which will make the first of the ones in the diagonal matrix we seek.
The next step is to remove the x value in equation two, this is done by adding the first equation,
multiplied by 3, to the second equation:
The second equation has now zero in the x value, and this must also now be done with the third equation as well, by adding
two times the first equation:
To form the unit matrix on equation two, we would need to multiply it by 2:
To form the last zero we would need to eliminate the y in the third equation, this is simply done by withdrawing 2 times the
second equation:
Now the last equation needs to be multiplied by -1, and we get the solution for the z value z=-1:
This completes what is often called forward shifting, as the matrix has all zeroes below the diagonal ones. It is now time to start what is called backwards elimination.
This starts off with the example by eliminating the y value from the first equation, and this is done simply by withdrawing the second equation multiplied by 1/2:
Now remove the z value from the first equation by simply adding the third equation with the first:
And at last add the third equation multiplied by -1 from the second and we get:
And thereby we get the solution which is of course x=2, y=3, and z = -1. This is quite a detailed run through of the Gaussian Elimination process,
and you could comprehend that this could be done automatically, regardless of how many variables there were. It should be mentioned that
we could always do this manually.
Gaussian elimination in code
Now we proceed to get the forward elimination done in code. We will skip the first row completely, which means that the first equation remains intact in the implementation
of forward Gaussian elimination.
First we take the second equation (or following values) and divide this by the x value of the first equation and the value of x, or rather always the actual value currently
in the unity matrix. The code below describes the implementation with the matrix called a, and the right hand side vector is called b.
It is important to notice that this would set all values to zero below the unity matrix. This would actually also give the solution to the last value, in the example given
in the previous section, the value of z. From this information we could now calculate all the variables, the reason is that the last equation is solved for the last variable,
the next equation has one unknown to be determined, as we can insert the last value into it, and so on:
Public Function GaussElimination(ByVal ResultingVector As Matrix) As Matrix
If Not (Me.Cols = ResultingVector.Rows And ResultingVector.Cols = 1) Then
If Me.Cols <> ResultingVector.Rows And ResultingVector.Cols <> 1 Then
Throw New Exception("The number of columns (variables) is different " & _
"form the number rows (number of equations), " & _
"and there is too many columns in the ResultingVector matrix.")
ElseIf Me.Cols <> ResultingVector.Rows Then
Throw New Exception("The number of columns (variables) is different " & _
"form the number rows (number of equations).")
ElseIf ResultingVector.Cols <> 1 Then
Throw New Exception("There is too many columns in the ResultingVector matrix.")
End If
End If
Dim a As New Matrix(Rows, Cols)
For i As Integer = 0 To Rows - 1
For j As Integer = 0 To Cols - 1
a.Item(i, j) = Me.Item(i, j)
Next
Next
Dim result, b As New Matrix(Rows, 1)
For i As Integer = 0 To Rows - 1
b.Item(i, 0) = ResultingVector(i, 0)
Next
Dim m, Sum As Double
For k As Integer = 0 To Rows - 2
For i As Integer = k + 1 To Rows - 1
m = a.Item(i, k) / a.Item(k, k)
For j As Integer = 0 To Cols - 1
a(i, j) -= m * a.Item(k, j)
Next
b(i, 0) = b(i, 0) - m * b(k, 0)
Next
Next
For i As Integer = Rows - 1 To 0 Step -1
Sum = 0
result(i, 0) = 0
For k As Integer = i + 1 To Rows - 1
Sum += a(i, k) * result(k, 0)
Next
result(i, 0) = (b(i, 0) - Sum) / a(i, i)
Next
Return result
End Function
The values are now stored in the result vector, and voila this is the Gaussian Elimination.
The problems with Gaussian Elimination
Gaussian Elimination is almost never implemented directly as seen from the code
snippets in the previous sections, as the matrix of the equation can both have zeroes in the unity matrix
before calculation, and can also have zero in the unity matrix, when subtraction of previous equations are performed. Both of these problems are solved using pivoting,
meaning that the order of equations in the matrix is rearranged, so that the problem doesn't occur. It is important to notice that the problem of zeroes can happen in each iteration,
where we subtract one equation from another, like the example below:
We subtract the equation in the second row with a half time the first row, and get:
Actually the problem of the possible zeros in the unity matrix can be coupled with the second problem with Gaussian Elimination, the problem of accuracy
of the solution. We can demonstrate this with an example:
Even with double accuracy we do not obtain the exact value (x=1, y=1, and z=1), and Gaussian Elimination does suffer from an extreme requirement
of numbers of digits in the calculation. This could be solved using fraction calculations with the use of the class
System.Numerics.BigInteger instead of Decimal
or Double in the calculations, this would eliminate errors of calculations, even for relatively small numbers.
We should mention here that we could use either
the Gauss-Jordan or better the Gauss-Seidel iteration schema instead.
However when pivoting a matrix is performed, we should bear in mind what it really is.
We could, as suggested previously, exchange the order of the equations,
meaning exchange the row in the equation. This is referred to as a partial pivoting, as
we can also change the order of the variables. If we exchange
the rows and the variables this is referred to as full pivoting. Assuming that we have a number of equations that are put in a matrix of dimensions n * n,
we have exactly n! ways of organizing the rows, and if we exchange the variables also
we would have n! * n!.
Before we make the pivoting, there is the question of what way is most preferable of organizing them in. The usual way is to take the highest
absolute value and place that element in the diagonal in the matrix, this is done to minimize the round off errors, and it is this way I shall organize the matrix in.
This would also be the way that each element would be organized in also.
The sorting is done with the so-called bubble sorting, this is usually an inefficient way, as we only need the highest value in the first element, so one could search for the max value and just swap that as the new top element. The code is written the following way:
Public Shared Function SortByValue(ByVal A As FractionsMatrix, _
b As FractionsMatrix, ByVal index As Integer) As Fraction
Dim n, i As Integer
n = A.Rows
i = index
Dim CurrentMaximumValue As Double = 0
Dim j As Integer = index
Do
If Math.Abs(A(j, index).ToDouble) < Math.Abs(A(j + 1, index).ToDouble) Then
A.SwapRows(j, j + 1)
b.SwapRows(j, j + 1)
j = index
Else
j += 1
End If
Loop Until j = n - 1
Return A(index, index)
End Function
It does not allow for every combination of the matrix to be unfolded, but is sufficient for most practical cases.
Fraction class
The main problem in constructing a fraction class, we immediately encounter the problem of constructing a fraction from either a
Decimal or Double.
There are two choices, one to construct it from a 10 base fraction. If so we thereby need the number of significant digits after the decimal place,
and that can be done as shown below were 'Number* is the 'decimal' or 'double':
Dim SignificantDigitCount As Integer = BitConverter.GetBytes(Decimal.GetBits(Number)(3))(2)
This would return the count of all digits until the rest is zeroes, but it is not the preferred technique to use. The continued fraction is the best way of getting a fraction
out of a decimal number, as it would get the exact correct fraction of the decimal number 0.333333333 etc., and that could not be done any other way. It will also converge more quickly
for an irrational number such as the numbers pi and e. We should however be careful with the implementation as a too early conversion from decimal to integer,
as we quickly get round off errors. The code to construct the continued fraction is given below:
Sub New(ByVal Number As Decimal)
Dim ResultingFraction As New Fraction
Dim k As Double = 1
If Number < 0 Then
k = -1
End If
ResultingFraction = GetFraction(Math.Abs(Number))
ResultingFraction.Numerator *= k
Me.Numerator = ResultingFraction.Numerator
Me.Denominator = ResultingFraction.Denominator
End Sub
Private Function GetFraction(ByVal d As Decimal)
Dim Temp As Decimal = d
Dim list As New List(Of Decimal)
For i As Integer = 0 To 1000
list.Add(Math.Truncate(Temp))
Dim ff As Decimal = GetContinuedFraction(list).ToDouble
If d = ff Then
Exit For
End If
Try
Temp = 1 / (Temp - Math.Truncate(Temp))
Catch ex As Exception
Exit For
End Try
Next
Return GetContinuedFraction(list)
End Function
Private Function GetContinuedFraction(ByVal d As List(Of Decimal)) As Fraction
Dim result As New Fraction
result.Numerator = d(d.Count - 1)
result.Denominator = 1
For i As Integer = d.Count - 2 To 0 Step -1
If Not result.Denominator = 0 Then
result = New Fraction(d(i), 1) + 1 / result
End If
Next
Return result
End Function
For the fraction must contain two values, namely the Numerator and
Denominator, to perform the calculations. However the integer value could become
quite large so we need to use System.Numerics.BigInteger.
The calculation of the fractions is trivial and we can look it up in the code to see how it
is done, however we should take care when we want to convert the fraction
into a double or decimal number. The BigInteger class could have values that are extremely large, in fact too large to be converted directly into a decimal place value.
If we really want to be sure we should account for it, but the maximum value for double allows for approximately 307 digits, so if we get
a number larger than that we need
to employ a code similar to this:
Public Function ToDouble() As Double
Dim NumberOfNumeratorDigits, NumberOfDenuminatorDigits, MaximumDoubleDigits As Integer
NumberOfNumeratorDigits = GetNumberOfDigits(Me.Numerator)
NumberOfDenuminatorDigits = GetNumberOfDigits(Me.Denumerator)
MaximumDoubleDigits = CInt(Math.Log10(Double.MaxValue)) - 1
Dim ResultingDoubleValue As Double
If MaximumDoubleDigits > NumberOfDenuminatorDigits And _
MaximumDoubleDigits > NumberOfNumeratorDigits Then
ResultingDoubleValue = CDbl(Me.Numerator) / CDbl(Me.Denumerator)
Else
If MaximumDoubleDigits < NumberOfNumeratorDigits And _
MaximumDoubleDigits < NumberOfDenuminatorDigits Then
Dim Over, Below As Integer
Over = NumberOfNumeratorDigits - MaximumDoubleDigits
Below = NumberOfDenuminatorDigits - MaximumDoubleDigits
Me.Numerator /= 10 ^ Over
Me.Denumerator /= 10 ^ Below
ResultingDoubleValue = CDbl(Me.Numerator) / CDbl(Me.Denumerator)
ResultingDoubleValue *= CDbl(10 ^ (Over - Below))
ElseIf MaximumDoubleDigits < NumberOfDenuminatorDigits Then
Dim Below As Integer
Below = NumberOfDenuminatorDigits - MaximumDoubleDigits
Me.Denumerator /= 10 ^ Below
ResultingDoubleValue = CDbl(Me.Numerator) / CDbl(Me.Denumerator)
ResultingDoubleValue *= CDbl(10 ^ (Below))
ElseIf MaximumDoubleDigits < NumberOfNumeratorDigits Then
Dim Over As Integer
Over = NumberOfNumeratorDigits - MaximumDoubleDigits
Me.Numerator /= 10 ^ Over
ResultingDoubleValue = CDbl(Me.Numerator) / CDbl(Me.Denumerator)
ResultingDoubleValue *= CDbl(10 ^ (Over))
End If
End If
Return ResultingDoubleValue
End Function
The number of digits are found with the following function:
Private Function GetNumberOfDigits(ByVal d As System.Numerics.BigInteger) As Integer
Dim NumberOfDigits As Integer
If d < 0 Then
d *= -1
End If
If d = 0 Then
NumberOfDigits = 1
Else
NumberOfDigits = System.Numerics.BigInteger.Log10(d) + 1
End If
Return NumberOfDigits
End Function
The BigInteger class has a function called gcd (greatest common divisor) which can be be used to minimize the fraction. The function will give 1 as a result
if the fraction is already in its shortest form (but can also give 0 as a result if the fraction is (0/0)).
Private Shared Function MinFraction(ByVal FratNumber As Fraction) As Fraction
Dim CommonDivisor As New System.Numerics.BigInteger
Do
CommonDivisor = System.Numerics.BigInteger.GreatestCommonDivisor(_
FratNumber.Numerator, FratNumber.Denumerator)
If CommonDivisor = 0 Then Exit Do
FratNumber.Numerator /= CommonDivisor
FratNumber.Denumerator /= CommonDivisor
Loop Until CommonDivisor = 1
Return FratNumber
End Function
Inverse matrix
The inverse matrix is really the same as solving a matrix, however it demands that it is on a square form, and that the determinant is not zero.
The steps for getting the inverse matrix, provided that it really exists, is nearly the exact same as the Gaussian elimination process.
We start with the same example as in the beginning, except that the constants are not included, with the matrix written
in the form:
What we now want is to multiply, add, or subtract each row, on both sides of the equation, so that we create the unity matrix on the left and get the inverse matrix on the right:
However, as with the Gaussian elimination, zeroes in the unity matrix is still an issue and would have to be solved by pivoting. This can be done by splitting
the unity matrix into three parts, and then using normal Gaussian elimination to solve each column and rearrange them into the unity matrix afterwords:
This process could of course be automated and would include that each column of the inverse matrix could be written like this in code:
Public Function Inverse() As Matrix
Dim result As New Matrix(Rows, Cols)
For i As Integer = 0 To Rows - 1
Dim v As New Matrix(Rows, 1)
v(i, 0) = 1
Dim res As New Matrix
Dim inputs As New Matrix(Rows, Cols)
For m As Integer = 0 To Cols - 1
For l As Integer = 0 To Rows - 1
inputs(l, m) = Me(l, m)
Next
Next
res = inputs.GaussElimination(v)
For k As Integer = 0 To Rows - 1
result(k, i) = res(k, 0)
Next
Next
Return result
End Function
Determinant
Determinant is an important calculation type for any matrix library, as it is a well understood property and there are fast and numerical stable ways of calculating it. This has lead to entire books written about the use of determinants and their meaning and implications. One of the places a programmer meets determinant calculations are with geometrical tests, were most of the equations could be rewritten to utilize the calculated determinant value.
Some examples of places were it is used in geometry functions are:
- Which side of the line (2D) a point lies (I used this in this implementation, and it can also be used to calcualte a 2D convex hull)
- Which side of a triangle (3D) a point lies
- Point inside circle (2D)
- Point inside sphere (3D)
- Cramers rule of solving the system of equations
- etc.
I will use forward substitution to calculate the determinant of the matrix (one out of three options, the other two are either by formula or by the use of the minor matrix together with the cofactor.), but in order to do this one must take a look at the properties of the determinant. Some of the most important properties are:
- The determinant of a matrix remains the same if the matrix is transposed
|A| = |A<sup>T</sup>| - If one swaps two rows, or two columns, the determinant changes sign, from
|A| to -|A| - If a row or columns only has zeroes the determinant will also be zero.
The determinant is calculated using forward substitution, i.a. the same as in native Gaussian Elimination, so the rules above must be applied when calculating the determinant, especially important is the changing of sign when two columns or rows are interchanged. The last criteria, were all the values in one row, or one column is zero (meaning that the determinant will be zero), must be checked for in advance, as the calculation could give NaN as a result with this configuration.
When the forward substitution is completed one is left with a triangular matrix, and the determinant is then just the product of the main diagonal. The code becomes:
Private Shared DeterminantSign As Integer = 1
Public Function Det() As Double
DeterminantSign = 1
If Not (Me.Cols = Me.Rows) Then
Throw New Exception("The number of columns (variables) is different form the number rows (number of equations), and there is too many columns in the ResultingVector matrix.")
End If
Dim a As New FractionsMatrix(Rows, Cols)
For i As Integer = 0 To Me.Rows - 2
SortByValue(Me, i)
Next
For i As Integer = 0 To Rows - 1
For j As Integer = 0 To Cols - 1
a.Item(i, j) = Me.Item(i, j)
Next
Next
Dim m, Sum As New Fraction
Dim HitZero As Boolean = False
For f As Integer = 0 To Rows - 1
If f <> 0 Then
For i As Integer = 0 To Rows - 1
For j As Integer = 0 To Cols - 1
a.Item(i, j) = (Me.Item(i, j))
Next
Next
DeterminantSign *= -1
a.SwapRows(0, f)
HitZero = False
End If
For k As Integer = 0 To Rows - 2
For i As Integer = k + 1 To Rows - 1
m = a.Item(i, k) / a.Item(k, k)
For j As Integer = 0 To Rows - 1
a(i, j) -= m * a.Item(k, j)
Next
Next
If FractionsMatrix.SortByValue(a, k).ToDouble = 0 Then
HitZero = True
Exit For
End If
Next
If HitZero = False Then
Exit For
End If
Next
Dim result As New Fraction(1, 1)
For i As Integer = 0 To Me.Cols - 1
result *= a.Item(i, i)
Next
Return result.ToDouble * DeterminantSign
End Function The only thing missing now is to check for rows and columns that is filled with zeroes, as given in code below:
''' <summary>
''' Check for row or colums containing nothing but zeroes
''' </summary>
''' <param name="a"></param>
''' <returns></returns>
''' <remarks></remarks>
Public Function IsDetermenantZero(ByVal a As FractionsMatrix) As Boolean
Dim result As Boolean = True
For i As Integer = 0 To a.Cols - 1
result = True
For j As Integer = 0 To a.Rows - 1
If a.Item(j, i).todouble <> 0 Then
result = False
End If
Next
If result Then
Return True
End If
Next
For i As Integer = 0 To a.Rows - 1
result = True
For j As Integer = 0 To a.Cols - 1
If a.Item(i, j).todouble <> 0 Then
result = False
End If
Next
If result Then
Return True
End If
Next
Return result
End Function
This concludes the determinant calculation.
History
26.08.2013 - Added determinant calculation
References
- "Numerical simulations and case studies using visual C++ .NET" by Shaharuddin Salleh et al
- "Practical numerical methods with C#" by Jack Xu
- "Statistics and data analysis in geology, Third Edition, by John C. Davis
- "Numerical Recipes in C++" Second Edition, by William Press et. al.
- "Real time collision detection - series in interactice 3D technology", by Christer Ericson
