# Solving Systems of Linear Equations March 25, 2012

Yet another university-related post, but I really enjoyed it so I thought I'd share: for a GUI- workshop I'm taking, we are given GUI-layout constraints as a system of linear equations, which we need to satisfy. To make life more interesting, some constraints are constant while some are parametric. There's no magic here, just some linear algebra combined with Python's overloadable nature to produce a nice and compact solver for these linear systems.

Naturally, you present your system as a *coefficient matrix*, which is Gauss-eliminated and then
"solved" for a given set of variables. I took the Gauss-Jordan elimination code from
Jarno Elonen and modified it to
support MxN matrices (not necessarily square). Here's a simple example of elimination:

```
>>> m = Matrix([1,2,4,2], [3,7,6,8], [2,2,2,9])
>>> print m.eliminate()
( 1.00 -0.00 0.00 6.00)
( 0.00 1.00 -0.00 -1.00)
( 0.00 0.00 1.00 -0.50)
```

But of course a reduced row-echelon matrix is not our final goal - we want to get the variable
*assignments*. For this, we have `solve()`

:

```
>>> sol = solve(m, ["x", "y", "z"])
>>> print sol
{'y': -1.0000000000000004, 'x': 6.0, 'z': -0.4999999999999998}
```

Modulo precision errors, we get `x = 6`

, `y = -1`

and `z = -0.5`

. If we have more equations than
variables, the "extraneous" equations must be linear combinations of previous ones, or a
contradiction will result. But what if we have less equations than variables? It means we have
some *degrees of freedom*... how would we handle that? It's actually simple - instead of being
resolved to constant values, variables will be assigned *dependent expressions*:

```
>>> m2 = Matrix([1,2,4,2], [3,7,6,8])
>>> print m2.eliminate()
( 1.00 0.00 16.00 -2.00)
( 0.00 1.00 -6.00 2.00)
>>>
>>>
>>> sol = solve(m2, ["x", "y", "z"])
>>> print sol
{'y': <BinExpr (2.0 - (-6.0 * z))>, 'x': <BinExpr (-2.0 - (16.0 * z))>,
'z': <FreeVar z>}
```

As you can see now, `z`

is a *free variable* and `x`

and `y`

are *dependent* on it. Of course more
than one variable may be free and some variables may be independent of free variables. Once a value
for `z`

is known, we can "fully evaluate" the dependent expressions:

```
>>> sol["x"].eval({"z" : 10})
-162.0
```

The code is available on my github page.
Note: all numbers are represented as `Decimals`

, to avoid loss of precision as much as possible,
and I'm using an "epsilon" value of `1E-20`

to equate numbers to each other (meaning, `x == y`

iff
`abs(x-y) <= epsilon`

).