Solving equations of upper triangular matrices

Suppose you have an equation like \(R \cdot x = b\) with \(R \in \mathbb{R}^{n \times n}\) and \(x,b \in \mathbb{R}^n\). \(b\) and \(R\) are given and you want to solve for \(x\).

Example

With \(n=5\), the problem could look like this:

This is only a shorthand for:

First step: Solve for $x_5$

First you see that \(x_5 = \frac{5}{4}\). So you divide \(b\) by the current row.

Now you replace every occurrence of \(x_5\) in the system of equations above:

Now you make the multiplications and remove the first trivial line.

Second step: update

Get the constant factors to the right side of the equations:

You're now in the same situation as in the first step. Next you will solve for \(x_4\), then for \(x_3, x_2\) and finally for \(x_1\).

This is called "back substitution".

According to Wolfram|Alpha, the solution is:

Python straightforward algorithm

I will use fractions for operations as I don't want to lose precision while dividing.

#!/usr/bin/env python
# -*- coding: utf-8 -*-

def solveUpperTriangularMatrix(R, b):
    # Convert R and b to Fraction
    from fractions import Fraction
    fR, fb = [], []
    for x, line in enumerate(R):
        fLine = []
        for y, el in enumerate(line):
            fLine.append(Fraction(el))
        fR.append(fLine)
    for el in b:
        fb.append(Fraction(el))

    # The solution will be here
    x = [Fraction(0)] * len(b)
    for step in range(len(b)-1, 0-1, -1):
        if fR[step][step] == 0:
            if fb[step] != 0:
                return "No solution"
            else:
                return "Infinity solutions"
        else:
            x[step] = fb[step] / fR[step][step]

        for row in range(step-1, 0-1, -1):
            fb[row] -= fR[row][step]*x[step]
    return x

if __name__ == "__main__":
    R = [[2, 7, 1, 8, 2],
         [0, 8, 1, 8, 2],
         [0, 0, 8, 4, 5],
         [0, 0, 0, 9, 0],
         [0, 0, 0, 0, 4]]
    b =  [3, 1, 4, 1, 5]
    x = solveUpperTriangularMatrix(R, b)
    print(x)

    # Convert x to float
    x = map(float, x)
    print(x)

A better algorithm

Just like for unipotent lower triangular matrices, we can operate directly on the given input:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

def solveUpperTriangularMatrix(R, b):
    # Convert R and b to Fraction
    from fractions import Fraction
    for x, line in enumerate(R):
        for y, el in enumerate(line):
            R[x][y] = Fraction(el)
    for x, el in enumerate(b):
        b[x] = Fraction(el)

    # The solution will be here
    for step in range(len(b)-1, 0-1, -1):
        if R[step][step] == 0:
            if b[step] != 0:
                return "No solution"
            else:
                return "Infinity solutions"
        else:
            b[step] = b[step] / R[step][step]

        for row in range(step-1, 0-1, -1):
            b[row] -= R[row][step]*b[step]

if __name__ == "__main__":
    R = [[2, 7, 1, 8, 2],
         [0, 8, 1, 8, 2],
         [0, 0, 8, 4, 5],
         [0, 0, 0, 9, 0],
         [0, 0, 0, 0, 4]]
    b =  [3, 1, 4, 1, 5]
    solveUpperTriangularMatrix(R, b)
    print(b)

Conversion to Fraction

You could think that the conversion to fraction is not necessary. But if you simply remove line 5 to 16, you will get:

[1, 0, -1, 0, 1]

because of integer arithmetic. When you convert the input to float before passing it to solveUpperTriangularMatrix, you will get

[0.8717447916666666, -0.25651041666666663, -0.3368055555555556, 0.1111111111111111, 1.25]

which is almost the same as when we calculated with Fraction and converted to float afterwards:

[0.8717447916666666, -0.2565104166666667, -0.3368055555555556, 0.1111111111111111, 1.25]

So: Using Fractions needs some computing time, but you will get better results.

Time complexity

I'll analyze the second algorithm.

The conversion of our input data is obviously in \(\mathcal{O}(n^2)\). Let's only analyse the part after the conversion.

Assume that there is exactly one solution and that line 15-21 take \(c_1\) operations and line 24 takes \(c_2\) operations.

Then we would have a total of

So the algorithms time complexity is in \(\Theta(n^2) \subsetneq \mathcal{O}(n^2)\).

Space complexity

Please note that I take advantage of Pythons dynamic typing system. I think it's difficult to see space complexity in python programs. But when you make the same in C++, you will see that you will need space in \(\mathcal{O}(n)\) when you do the conversion. Without the conversion, you're in \(\mathcal{O}(1)\).

I guess you might want to leave this choice to the user of your functions. When he wants better results, he should give the input as Fraction. When he wants to get results rather faster, he should give the input as float.

Notes

In this algorithm, we don't need anything below the diagonal.