Solving a System of Linear Equations

Problem

You have a number of linear equations (that is, equations that look like "2x + 10y + 8z = 54"), and you want to figure out the solution: the values of x, y, and z. You have as many equations as you have variables, so you can be certain of a unique solution.

Solution

Create two Matrix objects. The first Matrix should contain the coefficients of your equations (the 2, 10, and 8 of "2x + 10y + 8z = 54"), and the second should contain the constant results (the 54 of the same equation). The numbers in both matrices should be represented as floating-point numbers, rational numbers, or BigDecimal objects: anything other than plain Ruby integers.

Then invert the coefficient matrix with Matrix#inverse, and multiply the result by the matrix full of constants. The result will be a third Matrix containing the solutions to your equations.

For instance, consider these three linear equations in three variables:

	2x + 10y + 8z = 54 
	7y + 4z = 30 
	5x + 5y + 5z = 35

To solve these equations, create the two matrices:

	require 'matrix'
	require 'rational' 
	coefficients = [[2, 10, 8], [0, 7, 4], [5, 5, 5]].collect! do |row|
	 row.collect! { |x| Rational(x) }
	end
	coefficients = Matrix[*coefficients]
	# => Matrix[[Rational(2, 1), Rational(10, 1), Rational(8, 1)],
	# => [Rational(0, 1), Rational(7, 1), Rational(4, 1)],
	# => [Rational(5, 1), Rational(5, 1), Rational(5, 1)]]

	constants = Matrix[[Rational(54)], [Rational(30)], [Rational(35)]]

Take the inverse of the coefficient matrix, and multiply it by the results matrix. The result will be a matrix containing the values for your variables.

	solutions = coefficients.inverse * constants
	# => Matrix[[Rational(1, 1)], [Rational(2, 1)], [Rational(4, 1)]]

This means that, in terms of the original equations, x=1, y=2, and z=4.

Discussion

This may seem like magic, but it's analagous to how you might use algebra to solve a single equation in a single variable. Such an equation looks something like Ax = B: for instance, 6x = 18. To solve for x, you divide both sides by the coefficient:

The sixes on the left side of the equation cancel out, and you can show that x is 18/6, or 3.

In that case there's only one coefficient and one constant. With n equations in n variables, you have n2 coefficients and n constants, but by packing them into matrices you can solve the problem in the same way.

Here's a side-by-side comparision of the set of equations from the Solution, and the corresponding matrices created in order to solve the system of equations.

	2x + 10y + 8z = 54 | [ 2 10 8] [x] [54]
	x + 7y + 4z = 31 | [ 1 7 4] [y] = [31]
	5x + 5y + 5z = 35 | [ 5 5 5] [z] [35]

If you think of each matrix as a single value, this looks exactly like an equation in a single variable. It's Ax = B, only this time A, x, and B are matrices. Again you can solve the problem by dividing both sides by A: x = B/A. This time, you'll use matrix division instead of scalar division, and your result will be a matrix of solutions instead of a single solution.

For numbers, dividing B by A is equivalent to multiplying B by the inverse of A. For instance, 9/3 equals 9 * 1/3. The same is true of matrices. To divide a matrix B by another matrix A, you multiply B by the inverse of A.

The Matrix class overloads the division operator to do multiplication by the inverse, so you might wonder why we don't just use that. The problem is that Matrix#/ calculates B/A as B*A.inverse, and what we want is A.inverse*B. Matrix multiplication isn't commutative, and so neither is division. The developers of the Matrix class had to pick an order to do the multiplication, and they chose the one that won't work for solving a system of equations.

For the most accurate results, you should use Rational or BigDecimal numbers to represent your coefficients and values. You should never use integers. Calling Matrix#inverse on a matrix full of integers will do the inversion using integer division. The result will be totally inaccurate, and you won't get the right solutions to your equations.

Here's a demonstration of the problem. Multiplying a matrix by its inverse should get you an identity matrix, full of zeros but with ones going down the right diagonal. This is analagous to the way multiplying 3 by 1/3 gets you 1.

When the matrix is full of rational numbers, this works fine:

	matrix = Matrix[[Rational(1), Rational(2)], [Rational(2), Rational(1)]]
	matrix.inverse
	# => Matrix[[Rational(-1, 3), Rational(2, 3)],
	# => [Rational(2, 3), Rational(-1, 3)]]

	matrix * matrix.inverse 
	# => Matrix[[Rational(1, 1), Rational(0, 1)],
	# => [Rational(0, 1), Rational(1, 1)]]

But if the matrix is full of integers, multiplying it by its inverse will give you a matrix that looks nothing like an identity matrix.

	matrix = Matrix[[1, 2], [2, 1]]
	matrix.inverse
	# => Matrix[[-1, 1],
	# => [0, -1]]

	matrix * matrix.inverse 
	# => Matrix[[-1, -1],
	# => [-2, 1]]

Inverting a matrix that contains floating-point numbers is a lesser mistake: Matrix#inverse tends to magnify the inevitable floating-point rounding errors. Multiplying a matrix full of floating-point numbers by its inverse will get you a matrix that's almost, but not quite, an identity matrix.

	float_matrix = Matrix[[1.0, 2.0], [2.0, 1.0]]
	float_matrix.inverse
	# => Matrix[[-0.333333333333333, 0.666666666666667],
	# => [0.666666666666667, -0.333333333333333]]

	float_matrix * float_matrix.inverse 
	# => Matrix[[1.0, 0.0], 
	# => [1.11022302462516e-16, 1.0]]