LU Decomposition

So while we were poking at the Black-Scholes PDE numerically, this kind of matrix operation kept showing up:

What we want is $x$, so all we had to do was find the inverse of $A$. But finding that inverse is hard, and coding it up is even harder. So, to make life… a little easier, we learn LU decomposition — or so they say~

To understand LU decomposition better, I’m gonna start from Gaussian elimination.

Gaussian elimination, middle-school style

Let me explain the idea at a middle-school level. You’ve all solved a problem like this in elementary school:

There are 2-wheeled bikes and 4-wheeled bikes, 100 of them total. But somebody ducked down — whoosh!!! — and for some reason counted only the wheels. It’s already weird that they got an accurate wheel count, but apparently the wheels added up to 354.

How many 2-wheeled bikes? How many 4-wheeled bikes?

(In the time it took to count the wheels, you could’ve just counted the bikes lol)

Anyway. When we were kids, we’d start by guessing 50/50, then think about the wheel count — “should there be fewer 4-wheelers? fewer 2-wheelers?” — and inch our way to the answer.

But Korea’s dreaded 8th graders fear nothing. They just slam equations down. The 8th grader slams them down like this:

Set up the equations. Then — call it elimination, call it substitution, whatever — they double the top equation and line things up like this:

And then write this:

And then they just barrel through to the next step like an automated machine:

Now this 8th grader grows up. Goes to college. Takes linear algebra (the required class in basically every STEM major). And now they’d solve it like this. (Of course nobody actually solves it this way lol)

Rewritten as:

Multiply the top row by 2:

Multiply the top row by $(-1)$:

Replace row 2 with (row 2 + the top row):

↓

And replace the top row with (top row + row 2):

↓

Multiply both rows by $1/2$ (and for row 1, also multiply by $-1$):

──── tangent warning ────

If you’re seeing this for the first time, you might be staring at the red arrows above going “…wait, why is that allowed?”

Because I didn’t get it the first time either, heh. The reason it’s allowed comes down to properties of equality: if you have an equality $A = B$, you can add or subtract the same thing $C$ on both sides.

So in this equation

think of it like this:

Now let’s clearly add the same number to both sides of the bottom equality. We’ll add $-200$ to both the left side and the right side:

But that $-200$ on the left of the bottom equation? Let’s rewrite it as a different number in disguise:

That’s how you can explain the red arrows: adding or subtracting $k$ times one row from another is allowed. lol

──── tangent over ────

The compressed version of why all this is allowed: because it’s an equality! That’s it. But OK, let me say it the fancier way.

This matrix finally turned into the form

But — what if (just let this wash over you) it had been

this matrix??? Through the same kind of “elimination”, it would’ve had to become:

Let me describe what the “form we’re trying to reach” looks like:

And of course this shape, just like above, is reached using only the properties of equality. There’s more detail in Systems of Linear Equations and Gaussian Elimination (Linear Algebra I Studied #9).

A matrix with those features is called a row echelon form (when every leading 1 also has only zeros above it, like the one above, it’s the reduced row echelon form, RREF). And to turn a matrix into that form using the properties of equality, you get exactly three moves:

1. Multiply a row by a (nonzero) constant.
2. Swap two rows (or swap them after scaling).
3. Add a constant multiple of one row to another row.

That’s it. Those are the only moves. (These three are called the transformation rules.) Yep yep yep — using only these to bring a matrix into row echelon form is what we call Gaussian elimination.

Don’t lose the thread, by the way: we started from Gaussian elimination in order to learn LU decomposition.

EROMs — the matrices that do the elimination for you

Whoa!!! But!!! Here’s the thing about those transformation rules.

Swapping two rows:

Doesn’t matter, right?? Even doing it like this??? But you can also say that this swap means multiplying by the matrix $\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$~~~:

Picking one row of a matrix and multiplying it by a constant:

Doing that doesn’t break the system of equations, right??? But that scaling is also the same as multiplying by the matrix $\begin{pmatrix} r & 0 \\ 0 & 1 \end{pmatrix}$~~~~:

Finally, adding a constant multiple of one row to another row:

Doing this doesn’t break the system either, right?!? But the move is also the same as multiplying by $\begin{pmatrix} 1 & 0 \\ r & 1 \end{pmatrix}$~~~:

These transformation rules — even for a big~~~ matrix — are really “pick two rows out of many and do that thing.” That’s why I drew them above as $2 \times n$ matrices with two rows. The point I’m making here is not “eyeball it and figure it out.” The point is:

The act of left-multiplying by a certain matrix is itself the act of applying a transformation rule.

That’s what I want to land.

cf. You might be thinking “OK but does this work for big matrices?” To help, here’s a 4×4 example. To swap row 2 and row 3 in a 4×4 matrix:

— same as multiplying by this matrix. To multiply row 3 by $r$:

— same as multiplying by this matrix. To add $r \times$ row 1 to row 4:

— same as multiplying by this matrix. Keep that in mind for matrices with more than 2 rows, too.

In short, a special matrix that carries out an elementary row operation like this is called an Elementary Row Operation Matrix (EROM), .

From row echelon form to upper-triangular to LU

Ahhhhh, so. If you keeeeeeep multiplying EROMs onto some wild messy not-pretty-at-all matrix (= keep applying transformation rules), it turns into a pretty row echelon form.

But!!! What if that row echelon form is a square $n \times n$ matrix???

Then can’t we call it an upper triangular matrix?!?!

In an equation made of matrices, , say matrix $A$ reaches row echelon form after $k$ transformation rules. That is — say it becomes upper triangular after $k$ rules. So: — multiplying these EROMs a finite number of times,

— like this — turns $A$ (if it’s an $n \times n$ square matrix) into an upper triangular matrix. We just confirmed that.

So let’s take this guy, , as one big lump and call it . Then we’ve got .

Now — each of those individual $R$’s is, no matter what, an invertible matrix.

(Why each one is unconditionally invertible — explaining that would take forever T_T. I’ll leave that to linear algebra… probably you’d want to start from my “Linear Algebra I Studied #8” onward. But to follow #8 you’d want the stuff before #8 too, which is why I’m hesitating to drop a link… so I won’t… heh heh.)

Anyway, so:

what this means is — this is possible.

For LU, stop after forward elimination, once the matrix is upper triangular. Do not keep eliminating above the pivots to reach RREF. For now, assume the pivots are nonzero and no row swaps are needed.

Woohoo…… OK now — since becomes a lower triangular matrix, is the conclusion, so it looks like all I have left to explain is why becomes a lower triangular matrix.

First — the row-scaling and below-pivot row-addition EROMs we use for forward elimination are lower triangular. → Scroll back up and check!! heh heh. General row operations need not have this shape.

Then two facts:

1. The inverse of a row swap is that same row swap.

   

2. The inverse of a nonsingular lower triangular matrix is also lower triangular.

OK so here’s how it really has to be phrased. In , the operations can include row swaps. Keeping track of those swaps gives $PA=LU$, so we solve $LUx=Pb$. We cannot simply move all swaps to the front without adjusting the other factors.

Ah but — if we rename the suitably permuted system , then this renamed $A$ has an $LU$ decomposition when its elimination pivots are nonzero! The original matrix need not have an LU decomposition without row swaps.

Yayyyy~~~ $A$ → decomposable into $LU$!!!!

Why we wanted this in the first place

Quick reminder of why we’re learning LU decomposition in the first place. In this matrix equation, , using the inverse, — this is hard to code.

So instead, in our goal is just to find vector $x$ anyway, so if we change that equation into the form , we can find $x$ without ever coding an inverse matrix.

Now the actually hardcore part

OK NOW it gets really hellish lol

What I mean is: given a matrix $A$, we need to write code that computes the elements of $L$ and $U$ in the LU decomposition. First —

we saw this matrix being LU-decomposed like this:

(Because if you do Gaussian elimination, the leading 1s land like that in the $U$ matrix?? That is — since I’m assuming the diagonal of $U$ is 1, it’d be right to call this the Crout method.)

Let’s look at how to code $L$.

What this is saying right now is

— so we can express $l_{ij}$ in terms of $u_{ij}$ and $a_{ij}$!!!

I mean, if you brute-force Gaussian elimination on $Ax = b$, it becomes $Ux = (RR \dots R)b$, and at that point we already know both $u$ and $a$, right? So we don’t know exactly what $(RRR \dots R)$ looks like, and in these trying times finding its inverse, the lower triangular matrix, is a pain — what this means is, we can just slam $L$ in using only $u_{ij}$ and $a_{ij}$.

Ah but, that would mean we’d have to do Gaussian elimination on $A$ by hand first and then code, which is also kinda not right.

So from here on, I’m gonna show you what I’ll call the swoosh-swoosh-swoosh calculation method. Instead of finding everything in one big BANG!!!~~~ — one swoosh from $L$, one swoosh from $U$, then one more swoosh from $L$, one more swoosh from $U$ — if we go sequentially like that, we can take matrix $A$ and rrrrrrrrrroll right through it, defining every element of $L$ and $U$.

Might be hard to grok at first, but look at it carefully once and it’ll click lol. OK let’s go. I’ll keep it easy with a 3×3 like above, then generalize to $n \times n$.

Looking at this equation first,

— these relations fall right out. Putting them all together, let’s call it .

So now what’s left is the $(2,2)$, $(3,2)$, $(3,3)$ entries.

I organized things for $L$; let’s organize $U$ too.

Looking at the figure again:

As you can see, $l$ and $u$ are tangled into each other. So when coding, “BANG! define all of $L$, then BANG!!! define all of $U$” just won’t cut it for solving $Ax = b$ with LU decomposition.

So what does it mean to use the swoosh-swoosh-swoosh sequential method? We already know we can rewrite it as $LUx = b$. And while we don’t know $L$ and $U$ exactly yet,

— this much, we do know.

So what’s the answer;;;;

When you fill in these blanks — don’t fill in all of $L$~~~ and then all of $U$ — don’t do that.

One from $L$, one from $U$. One from $L$, one from $U$.

Find them in that direction, and by rrrrrrrrrolling through it~~~~~ that way, you fill in all of $L$ and $U$.

As this equation shows:

• Column 1 of $L$ can be found from $A$ alone.
• Row 1 of $U$ can be found entirely from the $(1,1)$ element of $L$ and $A$.
• Column 2 of $L$ needs row 1 of $U$. But row 1 of $U$ was already conquered, so we can use it.
• Look at row 2 of $U$ too — the same way, the materials for it were prepped just before. heh heh.

Kyahaha! Aha-ha-ha.

So now it makes sense that you can code your way from $Ax = b$ all the way to $LUx = b$. Alright — all that’s left is heading toward the solution, vector $x$. Let’s cover just that and then go straight to how the real coding goes! heh heh.

Forward sub, then back sub

I mentioned this briefly above.

──── copy-paste warning ────

Quick reminder of why we’re learning LU decomposition in the first place. In this matrix equation, , using the inverse, — this is hard to code.

So instead, in our goal is to find vector $x$ anyway, so if we change that equation into the form $LUx = b$, we can find $x$ without ever coding an inverse matrix.

──── copy-paste over ────

Like I said, I mentioned it briefly before. Now that we’ve finished coding $LU$, let’s use the logic above and head toward the solution, vector $x$.

So first let’s look at it as $Ux = z$!!! That makes it $Lz = b$, right?!?!

We can see it in this form, and since we now know both $L$ and $U$,

Looking at row 2:

Oh nice — for row 2, since we could find $z_1$ just before, we can find $z_2$ too~~~ I’ll do just row 3 and then generalize.

Now generalizing to $z_i$,

we can generalize it like this!!! From $z_1$ all the way to $z_n$, sequentially, all found!!!~~~

And now let’s head to $Ux = z$.

We can see it looks like this!!! This time, unlike before, it’s an upper triangular matrix — that is, the bottom row has the most zero terms!!! So we’ll go the way of finding things from the bottom row up, one at a time. (We currently have all of $z_1$ through $z_n$.)

Ooh, $x_n$ we can pin down!!! Let’s try the row above lol:

Ahh~~ since we have $x_n$, this is a form where we can find $x_{n-1}$ too!!!

But right now I can’t say it in concrete numbers — I have to keep saying $n$-this and $n$-that, which is a real headache. And generalizing from right here is a bit much, so let me do just one more, OK lol hoo hoo hoo (T_T). Next up is the row above that, row $n-2$:

This is so painfully annoying lol hoo hoo… now quick, let’s generalize. Generalize!

Whoa — the principle is finally done???

The actual code

So now… it’s time to code. But honestly… could I ever have come up with this code myself….? (T_T) I think this is gonna be more of an “understanding the code you’re given” kind of session!!!

This is a function procedure: toss in $A$ and $b$ from $Ax = b$, and out pops here’s your solution $x$~~!!! Here it is!

This teaching routine does not do pivoting: it needs nonzero, reasonably sized pivots. A nonsingular matrix can still fail here if a pivot is zero; use pivoted LU for general problems. The VBA version also assumes arrays start at index 0.

Function Solve_SE(matA() As Double, matb() As Double) As Double()

Dim rowA As Long, rowB As Long, colA As Long, n As Long, _
    j As Long, i As Long, k As Long, No As Long
Dim L() As Double, U() As Double, Lsum As Double, Usum As Double, _
    modB() As Double, sumofB As Double, x() As Double, sum As Double, matx As Double

rowA = UBound(matA(), 1)
colA = UBound(matA(), 2)
rowB = UBound(matb(), 1)

If (rowA <> colA) Or (rowB <> colA) Then Exit Function

ReDim L(rowA, colA) As Double
ReDim U(rowA, colA) As Double

For j = 0 To UBound(L, 1)
    For i = 0 To UBound(L, 1)
        If j = 0 Then
           L(i, j) = matA(i, j)
           U(j, i) = matA(j, i) / L(0, 0)

        Else
           If i < j Then
              L(i, j) = 0
           Else
              Lsum = 0
              For k = 0 To j - 1
                  Lsum = Lsum + L(i, k) * U(k, j)
              Next
              L(i, j) = matA(i, j) - Lsum
           End If

          If i < j Then
             U(j, i) = 0
          ElseIf i = j Then
             U(j, i) = 1
          Else
             Usum = 0

             For k = 0 To i - 1
                 Usum = Usum + U(k, i) * L(j, k)
             Next
             U(j, i) = (matA(j, i) - Usum) / L(j, j)
          End If
       End If
    Next
Next

ReDim modB(UBound(matb)) As Double
modB(0) = matb(0) / L(0, 0)

For i = 1 To UBound(modB)
    sumofB = 0
    For k = 0 To i - 1
        sumofB = sumofB + L(i, k) * modB(k)
    Next
    modB(i) = (matb(i) - sumofB) / L(i, i)
Next

ReDim x(UBound(modB)) As Double

For No = UBound(x) To 0 Step -1
    sum = 0
    For i = 0 To UBound(x)
        sum = sum + U(No, i) * x(i)
    Next
    x(No) = modB(No) - sum
Next

Solve_SE = x()

End Function

There’s no way to fit this in a single screenshot (T_T), and for it to look like code it has to be a screenshot, so here it is torn into pieces:

Honestly, I think studying code is something you do on your own, heave-ho, heave-ho — but please just take this as a reference for how I studied this code… T_T

Continuing with the code:

It works in Python too, with the same result as in VBA.

import numpy as np

def Solve_SE(matA, matb):
    rowA = len(matA); colA = len(matA.T);     rowB = len(matb)
    if (rowA != colA) | (rowB != colA):
        raise NotImplementedError
    else:
        pass
    L = np.zeros((rowA, colA))
    U = np.zeros((rowA, colA))
    for j in range(0, len(L.T)):
        for i in range(0, len(L)):
            if j==0:
                L[i, j] = matA[i, j]
                U[j, i] = matA[j, i] / L[0, 0]
            else:
                if i < j:
                    L[i, j] = 0
                else:
                    Lsum = 0
                    for k in range(0, j):
                        Lsum = Lsum + L[i, k]*U[k, j]
                    L[i, j] = matA[i, j] - Lsum
                if i < j:
                    U[j, i] = 0
                elif i == j:
                    U[j, i] = 1
                else:
                    Usum = 0
                    for k in range(0, i):
                        Usum = Usum + L[j, k]*U[k, i]
                    U[j, i] = (matA[j, i] - Usum) /L[j, j]
    modB = np.zeros((len(matb), 1))
    modB[0] = matb[0] / L[0, 0]
    for i in range(0, len(modB)):
        sumofB = 0
        for k in range(0, i):
            sumofB = sumofB + L[i, k]*modB[k]
        modB[i] = (matb[i] - sumofB) / L[i, i]
    x = np.zeros((len(modB), 1))
    for no in range(len(x)-1, 0-1, -1):
        sum=0
        for i in range(0, len(x)):
            sum = sum + U[no, i] * x[i]
        x[no] = modB[no] - sum
    return x

matA=np.zeros((3, 3))
matb=np.zeros((3, 1))
matA[0, 0]=4; matA[0, 1]=-1; matA[0, 2]=0
matA[1, 0]=-1; matA[1, 1]=4; matA[1, 2]=-1
matA[2, 0]=0; matA[2, 1]=-1; matA[2, 2]=4
matb[0]=3
matb[1]=2
matb[2]=3
matx = Solve_SE(matA, matb)
print(matx)

Heads-up. The figures here are the original notes, slips and all. A few things to watch for:

• In the handwritten notes on the VBA code, three index slips: the L loop adds L(i,k) * U(k,j) (not L(i,j) * U(k,j)), and its branch runs for i >= j, the diagonal included; the U loop adds L(j,k) * U(k,i) and divides by L(j,j) (not L(j,i)); and the back substitution adds U(No,i) * x(i). The typed VBA in the post is right.
• In the Python screenshot, the U loop should add L[j, k]*U[k, i], not L[k, i]*U[j, k]. The printed [1, 1, 1] comes out right only because the test matrix is tridiagonal; the code block in the post has the corrected line.
• In the repeated copy-paste block, the system we change to is LUx = b, not Ax = b.
• For a general permutation matrix, P inverse is P transpose. A single row swap is its own inverse, as used here.