LU Decomposition - Linear Algebra

SECTION 1

Historical Context & Motivation

Imagine you have a big system of equations — maybe five, ten, or even a hundred equations with the same number of unknowns. Solving them one at a time would take forever. Mathematicians realized centuries ago that organizing numbers into rectangular grids called matrices could speed things up. Over time, they discovered clever ways to split a matrix into simpler pieces, making the whole process much faster.

LU Decomposition is one of these clever tricks. The letters L and U stand for "Lower" and "Upper" — two special triangle-shaped matrices. The idea grew out of a method you might already know: Gaussian elimination (named after the famous mathematician Carl Friedrich Gauss). LU Decomposition packages Gaussian elimination in a reusable form so you can solve many problems without starting over each time.

1810s • Gauss & Systematic Elimination

Carl Friedrich Gauss develops a systematic method for eliminating variables from systems of equations. This method — now called Gaussian elimination — lays the groundwork for LU Decomposition.

1938 • Banachiewicz Formalizes LU

Polish mathematician Tadeusz Banachiewicz formally describes the decomposition of a matrix into lower and upper triangular factors, giving the technique a clear mathematical framework.

1948 • Turing & Doolittle Refine the Method

Alan Turing and others publish stable algorithms for computing LU Decomposition on early computers. The Doolittle and Crout algorithms become standard approaches.

1970s–Today • Backbone of Scientific Computing

LU Decomposition becomes a core routine in software libraries like LAPACK. It powers everything from weather simulations to 3-D video game physics.

The central question LU Decomposition answers is: Can we break a matrix into two simpler matrices so that solving equations becomes as easy as working backward through a list? The answer is yes — and it turns out to be one of the most useful ideas in all of applied mathematics.

SECTION 2

Core Principles & Definitions

Before diving into the math, let's nail down a few key ideas. LU Decomposition rests on the concept that every "well-behaved" square matrix can be written as the product of two triangular matrices. Understanding what "triangular" means is the first step.

Lower Triangular Matrix (L)

A square matrix where all entries above the main diagonal are zero. Think of a staircase that only has steps on the bottom-left side. The diagonal entries are usually set to 1.

Upper Triangular Matrix (U)

A square matrix where all entries below the main diagonal are zero. It's the mirror image of L — steps only appear on the upper-right side.

Factorization (A = L × U)

LU Decomposition means finding L and U so that when you multiply them together, you get back the original matrix A. This is called a factorization — just like factoring 12 into 3 × 4.

Forward & Back Substitution

Once you have L and U, solving A·x = b splits into two easy steps: first solve L·y = b going forward, then solve U·x = y going backward. Each step is simple because one matrix is triangular.

✦

KEY TAKEAWAY

Think of LU Decomposition like taking apart a LEGO set. The original model (matrix A) is complex, but you split it into two simpler sub-assemblies (L and U). Each sub-assembly is easy to work with on its own. When you need to solve a new problem with the same structure, you don't have to rebuild from scratch — you just swap in new pieces.

SECTION 3

Visual Explanation

The diagram below shows how a 3×3 matrix A is split into a lower triangular matrix L and an upper triangular matrix U. Notice where the zeros sit in each triangle and how the non-zero entries form a staircase pattern.

The original 3×3 matrix A (left) equals the product of a lower triangular matrix L (center) and an upper triangular matrix U (right). The shaded triangles show where the non-zero values live.

In the diagram, the cyan-shaded triangle in L contains the multipliers that you would use during Gaussian elimination. The violet-shaded triangle in U is the row-echelon form (the result after you finish eliminating). By saving both pieces, you avoid re-doing elimination if you need to solve the same system with a different right-hand side.

SECTION 4

Mathematical Framework

Let's put the pieces together with equations. Suppose you have a system of linear equations written in matrix form: A · x = b, where A is the coefficient matrix, x is the vector of unknowns, and b is the vector of constants on the right side.

THE FACTORIZATION

A = L × U

A = original n × n matrix | L = lower triangular (1s on diagonal) | U = upper triangular

Because A = L × U, the equation A · x = b becomes (L × U) · x = b. We introduce a helper variable y and split the problem into two smaller problems.

STEP 1 — FORWARD SUBSTITUTION

L · y = b → solve for y

Start from the first row and work downward. Each row gives you the next value of y because L is lower triangular.

STEP 2 — BACK SUBSTITUTION

U · x = y → solve for x

Start from the last row and work upward. Each row gives you the next value of x because U is upper triangular.

Why is this better than just doing elimination every time? Imagine you need to solve A · x = b for ten different b vectors (ten different sets of constants). Without LU, you'd do Gaussian elimination ten times. With LU, you factor A once, then run the two quick substitution steps for each new b. The factorization is the expensive part; the substitutions are fast.

COMPUTING MULTIPLIERS

l_ij = a_ij / u_jj (for i > j)

Each entry l_ij in L is the multiplier you use to eliminate the entry in row i, column j during Gaussian elimination.

SECTION 5

Step-by-Step Breakdown of the Algorithm

Let's walk through exactly how you build L and U from a matrix A. The process mirrors Gaussian elimination, but you carefully record every multiplier along the way.

The flowchart traces the LU algorithm: for each column, compute the multipliers (stored in L), then eliminate entries below the pivot. After all columns, U is the resulting upper triangular matrix.

Pick the pivot: The pivot is the diagonal entry in the current column. For column j, the pivot is the entry at row j, column j.
Compute multipliers: For each row i below the pivot (i > j), divide the entry at (i, j) by the pivot. This fraction goes into L at position (i, j).
Eliminate: Subtract the multiplier times the pivot row from row i. This zeros out the entry below the pivot.
Move to the next column: Repeat the process for the next column until the matrix is upper triangular — that's U.

⚠

Watch Out for Zero Pivots!

If a pivot entry is zero, you can't divide by it. In practice, you swap rows first (called partial pivoting) to move a non-zero value into the pivot position. This gives a slightly modified decomposition: PA = LU, where P is a permutation matrix that records the row swaps.

SECTION 6

Worked Example

Let's decompose the following 3×3 matrix A into L and U, then use the result to solve a system.

GIVEN MATRIX

A = [ [2, 1, 1], [4, 3, 3], [8, 7, 9] ]

We also want to solve A · x = b where b = [4, 10, 24]^T.

LU Decomposition & Solve

1
Step 1 — Eliminate Column 1

The pivot is A(1,1) = 2. For Row 2, the multiplier is 4 ÷ 2 = 2. Subtract 2 × Row 1 from Row 2: [4, 3, 3] − 2 × [2, 1, 1] = [0, 1, 1]. For Row 3, the multiplier is 8 ÷ 2 = 4. Subtract 4 × Row 1 from Row 3: [8, 7, 9] − 4 × [2, 1, 1] = [0, 3, 5]. Store multipliers 2 and 4 in column 1 of L.

After Step 1: matrix becomes [ [2, 1, 1], [0, 1, 1], [0, 3, 5] ]. L column 1: l₂₁ = 2, l₃₁ = 4
2
Step 2 — Eliminate Column 2

The new pivot is the (2,2) entry = 1. For Row 3, the multiplier is 3 ÷ 1 = 3. Subtract 3 × Row 2 from Row 3: [0, 3, 5] − 3 × [0, 1, 1] = [0, 0, 2]. Store multiplier 3 in L.

U = [ [2, 1, 1], [0, 1, 1], [0, 0, 2] ]. L = [ [1, 0, 0], [2, 1, 0], [4, 3, 1] ]
3
Step 3 — Forward Substitution (L · y = b)

Solve L · y = [4, 10, 24]. Row 1: y₁ = 4. Row 2: 2(4) + y₂ = 10, so y₂ = 2. Row 3: 4(4) + 3(2) + y₃ = 24, so y₃ = 24 − 16 − 6 = 2.

y = [4, 2, 2]
4
Step 4 — Back Substitution (U · x = y)

Solve U · x = [4, 2, 2]. Row 3: 2x₃ = 2, so x₃ = 1. Row 2: x₂ + 1(1) = 2, so x₂ = 1. Row 1: 2x₁ + 1(1) + 1(1) = 4, so 2x₁ = 2 and x₁ = 1.

x = [1, 1, 1]

💡

Quick Check

Plug x = [1, 1, 1] back into A · x: [2(1)+1(1)+1(1), 4(1)+3(1)+3(1), 8(1)+7(1)+9(1)] = [4, 10, 24] = b. ✓ It checks out!

SECTION 7

Strengths & Limitations

LU Decomposition is powerful, but like any tool it has situations where it shines and situations where other methods might be better. The table below summarizes the key trade-offs.

Aspect	Strength ✓	Limitation ✗
Speed for multiple solves	Factor once, solve many times. Each new right-hand side only needs two substitution passes.	Initial factorization has the same cost as Gaussian elimination, so for a single solve there's no speed advantage.
Determinant & Inverse	The determinant of A is simply the product of the diagonal entries of U. Finding the inverse is also streamlined.	Computing a full inverse is still expensive for very large matrices.
Numerical stability	With partial pivoting (PA = LU), the method is stable for most practical problems.	Without pivoting, small rounding errors can blow up, especially for ill-conditioned matrices.
Matrix type	Works on any square matrix that doesn't need row swaps (or with pivoting, any non-singular square matrix).	Does not directly apply to non-square (rectangular) matrices. QR or SVD are used instead.

Strengths and limitations of LU Decomposition

✦

KEY TAKEAWAY

Think of LU Decomposition like pre-cooking a sauce base. The base takes time to prepare, but once it's done, you can quickly make many different dishes by adding different toppings. If you only need one dish, you might as well cook from scratch. But when you need to serve ten dishes with the same base, the upfront effort pays off enormously.

SECTION 8

Connection to Advanced Methods

LU Decomposition is just one member of a family of matrix factorizations. As you move into more advanced math and science, you'll encounter related decompositions designed for specific types of matrices or problems.

Decomposition	Form	When You'd Use It
LU	A = L × U	General square systems, finding determinants.
Cholesky	A = L × L^T	Symmetric positive-definite matrices (common in statistics). About twice as fast as LU.
QR	A = Q × R	Least-squares problems, non-square matrices, eigenvalue algorithms.
SVD	A = U × Σ × V^T	Data compression, image processing, machine learning. The most versatile decomposition.

How LU relates to other matrix factorizations

Mastering LU Decomposition gives you a strong foundation for all of these. The core idea — break a complicated matrix into simpler factors — runs through all of advanced linear algebra. In college-level courses and engineering applications, you'll see LU and its cousins working behind the scenes in simulation software, search engines, and even the recommendation algorithms on streaming services.

SECTION 9

Practice Problems

Question 1: CONCEPTUAL

In an LU Decomposition A = L × U, what special property do the diagonal entries of the lower triangular matrix L typically have, and why is this choice made?

Show answer

The diagonal entries of L are all set to 1. This convention is called the Doolittle form. We choose 1s on L's diagonal because it makes the decomposition unique — without this rule, there would be infinitely many ways to split A into L and U. Setting L's diagonal to 1 means all the "scaling" information lives in U's diagonal, keeping everything organized. It also simplifies the forward substitution step because multiplying by 1 doesn't change anything.

Question 2: BASIC CALCULATION

Find the LU Decomposition of the 2×2 matrix A = [ [3, 6], [1, 4] ]. Write out L and U.

Show answer

Start with the pivot A(1,1) = 3. The multiplier for Row 2 is 1 ÷ 3 = 1/3. Subtract (1/3) × Row 1 from Row 2: [1, 4] − (1/3)[3, 6] = [0, 2]. So U = [ [3, 6], [0, 2] ] and L = [ [1, 0], [1/3, 1] ]. You can verify: L × U = [ [1×3+0×0, 1×6+0×2], [(1/3)×3+1×0, (1/3)×6+1×2] ] = [ [3, 6], [1, 4] ] = A. ✓

Question 3: INTERMEDIATE

Given A = [ [1, 2, 4], [3, 8, 14], [2, 6, 13] ], perform LU Decomposition to find L and U.

Show answer

Column 1: Pivot = 1. Multiplier for Row 2: 3/1 = 3. Row 2 becomes [3,8,14] − 3×[1,2,4] = [0,2,2]. Multiplier for Row 3: 2/1 = 2. Row 3 becomes [2,6,13] − 2×[1,2,4] = [0,2,5]. Column 2: New pivot = 2 (in the updated Row 2). Multiplier for Row 3: 2/2 = 1. Row 3 becomes [0,2,5] − 1×[0,2,2] = [0,0,3]. Therefore U = [ [1,2,4], [0,2,2], [0,0,3] ] and L = [ [1,0,0], [3,1,0], [2,1,1] ]. Checking: L×U reproduces A. Each multiplier was placed in L at the position of the entry it helped eliminate.

Question 4: APPLIED

A video game physics engine needs to solve A · x = b for the same 3×3 matrix A = [ [2,1,1], [4,3,3], [8,7,9] ] every frame, but with a different b each frame. In Frame 1, b = [1, 1, 1]. Using the LU factorization from the worked example (L = [ [1,0,0],[2,1,0],[4,3,1] ] and U = [ [2,1,1],[0,1,1],[0,0,2] ]), find x for Frame 1.

Show answer

Forward substitution (L·y = b): y₁ = 1. Row 2: 2(1) + y₂ = 1, so y₂ = −1. Row 3: 4(1) + 3(−1) + y₃ = 1, so y₃ = 1 − 4 + 3 = 0. Thus y = [1, −1, 0]. Back substitution (U·x = y): Row 3: 2x₃ = 0, so x₃ = 0. Row 2: x₂ + 1(0) = −1, so x₂ = −1. Row 1: 2x₁ + 1(−1) + 1(0) = 1, so 2x₁ = 2 and x₁ = 1. The solution is x = [1, −1, 0]. This process took only the two substitution passes — no re-doing elimination — which is why game engines rely on LU for speed.

Question 5: CRITICAL THINKING

Consider the matrix B = [ [0, 1], [2, 3] ]. Can you perform LU Decomposition on B without row swaps? Explain why or why not, and describe what modification would fix the problem.

Show answer

No, you cannot decompose B into LU without row swaps. The first pivot B(1,1) = 0, and you cannot divide by zero to compute multipliers. This doesn't mean B is singular (its determinant is 0×3 − 1×2 = −2 ≠ 0), just that the standard LU algorithm breaks down. The fix is partial pivoting: swap Row 1 and Row 2 to get B' = [ [2,3], [0,1] ]. Now the pivot is 2 and the multiplier is 0/2 = 0, giving L = [ [1,0], [0,1] ] (the identity) and U = B'. We record the swap in a permutation matrix P = [ [0,1], [1,0] ], so the full decomposition is PA = LU. This shows that partial pivoting extends LU to handle any non-singular matrix.

SUMMARY

Summary

LU Decomposition splits a square matrix A into a lower triangular matrix L (with 1s on the diagonal) and an upper triangular matrix U so that A = L × U. This factorization captures the work of Gaussian elimination in a reusable form. The multipliers used during elimination fill the entries of L, and the resulting row-echelon matrix becomes U.

To solve a system A · x = b, you first perform forward substitution on L · y = b, then back substitution on U · x = y. The big advantage is that once L and U are found, solving for new right-hand sides is fast. When a zero pivot appears, partial pivoting (row swaps recorded in a permutation matrix P) keeps the method stable. LU is the foundation for more advanced factorizations like Cholesky, QR, and SVD — all built on the same idea of breaking matrices into simpler, more manageable pieces.