15,850,298 members
Articles / Programming Languages / Visual Basic

# High precision native Gaussian Elimination

Rate me:
26 Aug 2013CPOL12 min read 66.1K   2.3K   33   20
Solving linear equations and finding the inverse matrix using Gaussian elimination.

## 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:

VB
```Public Function GaussElimination(ByVal ResultingVector As Matrix) As Matrix

' Check that both the matrix dimensions are correct
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:

VB
```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':

VB
`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:

VB
```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
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:

VB
```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:

VB
```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)).

VB
```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:

VB
```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:

VB.NET
```  Private Shared DeterminantSign As Integer = 1

Public Function Det() As Double

DeterminantSign = 1

'Check that both the matrix dimensions are correct
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)

'Order the original equation set from
'the highest to lowest according to their unity value
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

'Forward elimination
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

'Check for zeroes here as there is no need to
If FractionsMatrix.SortByValue(a, k).ToDouble = 0 Then
HitZero = True
Exit For
End If
Next

' No zeroes found; resume to backward substitution
If HitZero = False Then
Exit For
End If
Next

Dim result As New Fraction(1, 1)

'Calculate the diagonal product
For i As Integer = 0 To Me.Cols - 1
result *= a.Item(i, i)
Next

'Return the diagonal product multiplyed by the changing sign given by row swaps
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.

## 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

Written By
Chief Technology Officer
Norway
This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.

 First Prev Next
 excellent Southmountain8-Feb-20 8:37 Southmountain 8-Feb-20 8:37
 Error in code? Member 1340131118-Sep-17 22:49 Member 13401311 18-Sep-17 22:49
 Re: Error in code? Member 1340131129-Sep-17 22:41 Member 13401311 29-Sep-17 22:41
 Re: Error in code? Kenneth Haugland11-Oct-17 8:50 Kenneth Haugland 11-Oct-17 8:50
 GetContinuedFraction Member 15656279-Sep-17 4:43 Member 1565627 9-Sep-17 4:43
 Big Num - and don't work Ilya Pronyashin31-Dec-14 9:15 Ilya Pronyashin 31-Dec-14 9:15
 Re: Big Num - and don't work Kenneth Haugland31-Dec-14 19:50 Kenneth Haugland 31-Dec-14 19:50
 This is what you always hope to find on the intrenet. jesneed13-Feb-14 6:32 jesneed 13-Feb-14 6:32
 Re: This is what you always hope to find on the intrenet. Kenneth Haugland14-Feb-14 7:28 Kenneth Haugland 14-Feb-14 7:28
 My vote of 5 Anshul Mehra18-Jul-13 17:40 Anshul Mehra 18-Jul-13 17:40
 Re: My vote of 5 Kenneth Haugland18-Jul-13 18:19 Kenneth Haugland 18-Jul-13 18:19
 My vote of 2 YvesDaoust16-Jul-13 4:29 YvesDaoust 16-Jul-13 4:29
 Re: My vote of 2 Kenneth Haugland16-Jul-13 10:42 Kenneth Haugland 16-Jul-13 10:42
 My vote of 5 Sergio Andrés Gutiérrez Rojas7-Jul-13 17:31 Sergio Andrés Gutiérrez Rojas 7-Jul-13 17:31
 Re: My vote of 5 Kenneth Haugland8-Jul-13 2:57 Kenneth Haugland 8-Jul-13 2:57
 Well done, a 5 Espen Harlinn7-Jul-13 15:01 Espen Harlinn 7-Jul-13 15:01
 Re: Well done, a 5 Kenneth Haugland8-Jul-13 2:56 Kenneth Haugland 8-Jul-13 2:56
 Last Visit: 31-Dec-99 19:00     Last Update: 2-Mar-24 20:11 Refresh 1