LU Decomposition

The problem, from first principles

10-gaussian-elimination solves Ax = b by row-reducing A down to an upper-triangular form. But that same row-reduction work is thrown away the moment you need to solve a different b with the same A — a common situation (simulating many right-hand sides against one fixed system matrix). LU decomposition captures the row-reduction work itself as two reusable matrices, so solving for a new b afterward is nearly free.

From theory to code

Implement lu_decompose(A), returning (L, U) such that A = LU, using the exact row-elimination steps from 10-gaussian-elimination but recording each elimination multiplier into L instead of discarding it. The signature and docstring are already in the editor.

Constraints

• A is a square n × n NumPy array that does not require row swaps (a partial pivoting-free case — real solvers handle the swap case too, but that's a refinement on top of this question's core idea, not covered here).
• L is lower-triangular with 1s on its diagonal; U is upper-triangular.
• L @ U must reconstruct A exactly (up to floating-point precision).

Examples

Example 1 (2 × 2)

A = [[4, 3],
     [6, 3]]

Zero out A[1, 0] using pivot row 0: the multiplier is 6 / 4 = 1.5. Subtract 1.5 * row 0 from row 1, giving [6 - 6, 3 - 4.5] = [0, -1.5].

L = [[1.0, 0.0],      U = [[4.0,  3.0],
     [1.5, 1.0]]           [0.0, -1.5]]

Check: L @ U = [[4, 3], [1.5 * 4, 1.5 * 3 - 1.5]] = [[4, 3], [6, 3]].

Example 2 (3 × 3)

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

• Column 0, pivot 2: multipliers 4 / 2 = 2 for row 1 and 8 / 2 = 4 for row 2. Row 1 becomes [0, 1, 1], row 2 becomes [0, 3, 5].
• Column 1, pivot 1: multiplier 3 / 1 = 3 for row 2. Row 2 becomes [0, 0, 2].

L = [[1, 0, 0],       U = [[2, 1, 1],
     [2, 1, 0],            [0, 1, 1],
     [4, 3, 1]]            [0, 0, 2]]

Each multiplier sits in L at the position of the entry it zeroed out: L[1, 0] = 2, L[2, 0] = 4, L[2, 1] = 3.

Hints