Click here to Skip to main content
Click here to Skip to main content

High precision native Gaussian Elimination

, 26 Aug 2013 CPOL
Rate this:
Please Sign up or sign in to vote.
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:

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:  

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| = |AT|
  • 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

        '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. 

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

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

Kenneth Haugland
Engineer
Norway Norway
I hope that you like the stuff I have created and if you do wish to say thank you then a donation is always appreciated.
You can donate here[^].

Comments and Discussions

 
GeneralThis is what you always hope to find on the intrenet. Pinmemberjesneed13-Feb-14 6:32 
GeneralRe: This is what you always hope to find on the intrenet. PinmemberKenneth Haugland14-Feb-14 7:28 
GeneralMy vote of 5 PinmemberAnshul Mehra18-Jul-13 17:40 
AnswerRe: My vote of 5 PinprofessionalKenneth Haugland18-Jul-13 18:19 
GeneralMy vote of 2 PinmemberYvesDaoust16-Jul-13 4:29 
GeneralRe: My vote of 2 PinprofessionalKenneth Haugland16-Jul-13 10:42 
GeneralMy vote of 5 PinmemberSergio Andrés Gutiérrez Rojas7-Jul-13 17:31 
GeneralRe: My vote of 5 PinprofessionalKenneth Haugland8-Jul-13 2:57 
GeneralWell done, a 5 PinmvpEspen Harlinn7-Jul-13 15:01 
GeneralRe: Well done, a 5 PinprofessionalKenneth Haugland8-Jul-13 2:56 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Terms of Use | Mobile
Web01 | 2.8.141220.1 | Last Updated 26 Aug 2013
Article Copyright 2013 by Kenneth Haugland
Everything else Copyright © CodeProject, 1999-2014
Layout: fixed | fluid