Python Linear Equations - Gaussian Elimination

Question

Goal

Given a set of points, I'm trying to find the coefficients of the linear equation that satisfies all the points provided.

For example, if I wanted to find the linear equation (ax + by + c = z):

3x + 2y + 2 = z

I would need a minimum of three, three-dimensional points:

(2, 2, 12)
(3, 4, 19)
(4, 5, 24)

Given enough points with coordinates (x, y, z), I should be able to find (a, b, c) using Gaussian Elimination.

However, I think I'm having issues back solving the matrix in special cases. You can view my first stab at a python implementation here: https://gist.github.com/anonymous/8188272

Let's take a look at a few examples...

Data Set 1

Use the following "hand made" points (x, y, z):

(2, 2, 12)
(3, 4, 19)
(4, 5, 24)

Perform LU decomposition on the following matrix:

[[  2.   2.   1.  12.]
 [  3.   4.   1.  19.]
 [  4.   5.   1.  24.]]

Back solve the U matrix:

[[  4.    5.    1.   24. ]
 [  0.   -0.5   0.5   0. ]
 [  0.    0.    0.5   1. ]]

The returned result (a, b, c):

[3.0, 2.0, 2.0]

Correct! Everything seems fine...

Data Set 2

Use the following "hand made" points (x, y, z):

(3, 4, 19)
(4, 5, 24)
(5, 6, 29)

Perform LU decomposition on the following matrix:

[[  3.   4.   1.  19.]
 [  4.   5.   1.  24.]
 [  5.   6.   1.  29.]]

Back solve the U matrix:

[[  5.00000000e+00   6.00000000e+00   1.00000000e+00   2.90000000e+01]
 [  0.00000000e+00   4.00000000e-01   4.00000000e-01   1.60000000e+00]
 [  0.00000000e+00   0.00000000e+00   4.44089210e-16   0.00000000e+00]]

The returned result (a, b, c):

[1.0, 4.0, 0.0]

While technically it is a solution, not what I was looking for!

Data Set 3

Use the following "hand made" points (x, y, z):

(5, 6, 29)
(6, 7, 34)
(7, 8, 39)

Perform LU decomposition on the following matrix:

[[  5.   6.   1.  29.]
 [  6.   7.   1.  34.]
 [  7.   8.   1.  39.]]

Back solve the U matrix:

[[  7.00000000e+00   8.00000000e+00   1.00000000e+00   3.90000000e+01]
 [  0.00000000e+00   2.85714286e-01   2.85714286e-01   1.14285714e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   3.55271368e-15]]

Implementation crashes...

Thoughts

In Data Set 2 and 3 the last row and second to last row are "special". The second to last row has the same value for "b" and "c" (which is true in my special example!). Unfortunately, I lack the mathematical knowledge to make heads or tails from it.

Is there something special case I need to handle when the last row is all zeros and the row above it has equal values across?

Thanks in advance!

Answer 1

Check out numpy.linalg, scipy.linalg and scipy.optimize.

For example using Data set 1

(2, 2, 12)
(3, 4, 19)
(4, 5, 24)

for 3x + 2y + 2 = z use numpy.linalg.solve

>>> a = np.array([[2, 2, 1],
                  [3, 4, 1],
                  [4, 5, 1]])
>>> b = np.array([12, 19, 24])
>>> np.linalg.solve(a, b)
array([ 3.,  2.,  2.])

For example using Data set 2

(3, 4, 19)
(4, 5, 24)
(5, 6, 29)

I get ...

>>> a = np.array([[3, 4, 1], [4, 5, 1], [5, 6, 1]])
>>> b = np.array([19, 24, 29])
>>> np.linalg.solve(a, b)
array([ 0.73333333,  4.26666667, -0.26666667])

Was this the answer you were looking for? This is one valid answer, but since a is rank 2, [2., 3., 1.] is also a valid answer. See below ...

For example using Data set 3

(5, 6, 29)
(6, 7, 34)
(7, 8, 39)

repeat process ...

>>> a = np.array([[5, 6, 1], [6, 7, 1], [7, 8, 1]])
>>> b = np.array([29, 34, 39])
>>> np.linalg.solve(a, b)
LinAlgError: Singular matrix

This means that the determinant of your coefficients is zero and therefore the inverse of the matrix is infinite, or in general terms, the system of equations has an infinite number of solutions or no solution.

>>> np.linalg.det(a)
0.0
>>> 1. / np.linalg.det(a)
inf

Remember you are solving a system ax = b by assuming that x = inv(a)b therefore inv(a) must exist and be finite. If a is singular, then inv(a) is infinite.

>>> np.linalg.inv(a)
LinAlgError: Singular matrix

So how about trying a least squares method to find the best solution.

>>> np.linalg.lstsq(a,b)
(array([ 2.,  3.,  1.]),       # solution "x"
 array([], dtype=float64),     # residuals, empty if rank > a.shape[0] or < a.shape[1]
 2,                            # rank
 array([  1.61842911e+01,   2.62145599e-01,   2.17200830e-16]))    # singular values of "a"

So it has found the best solution as [2., 3., 1.], which lucky you, is actually a solution to your conditions! The residuals are returned as empty, because the as @wim said, a is rank deficient, EG: not equal to the dimensions of square matrix a, or full rank.

Answer 2

Yes, it is a special case you need to handle differently. In cases 2 and 3 you have a rank deficient matrix. Generally it can mean there are infinitely many solutions, or no solutions.

You can identify whether these cases are going to occur by checking the determinant of the matrix you make by stacking those 3-vectors.

>>> import numpy as np
>>> from scipy.linalg import det
>>> data1 = np.array([(2, 2, 12), (3, 4, 19), (4, 5, 24)])
>>> data2 = np.array([(3, 4, 19), (4, 5, 24), (5, 6, 29)])
>>> data3 = np.array([(5, 6, 29), (6, 7, 34), (7, 8, 39)])
>>> det(data1)
-1.9999999999999982
>>> det(data2)
5.551115123125788e-17
>>> det(data3)
8.881784197001213e-16

Example 1 was a full rank matrix, which geometrically tells you that the 3 points are linearly independent.

Examples 2 and 3 make matrices with zero determinant, which tells you that the points are linearly dependent.

Answer 3

What you are looking for is the plane that contains all 3 points. The reason your solutions seem to act strangely for datasets 2 & 3 is that the three points are collinear, so no unique solution exists (i.e. there are an infinite number of planes that contain any given line). This is reflected in your LU decomposition with the appearance of a 'zero' row, since your matrix is rank 2.

Assuming you stick to finding the plane that contains all 3 points, you need to ensure that your matrix is actually rank 3. If it is, then you have a solution. If it is rank 2, then any plane which contains the common line is a valid solution.

Note: If you try to find a common plane for 4 points, then you may find that no such solution exists.