Matrix Operations

By the end of this lesson you will be able to solve a system of linear equations Ax = b the way real numerical libraries do: factor A once into PA = LU, then finish with two cheap triangular sweeps. You will trace lupSolve, luDecomposition and lupDecomposition line by line on small worked matrices (with the pivot swaps visible), see exactly why a zero pivot breaks plain LU and a singular matrix breaks everything, invert a matrix as n solves, and fit a curve through noisy data with the normal equations. Every algorithm is implemented in idiomatic Dart and checked against residuals A·x − b and random matrices.

0. Linear systems, triangular systems, and a note on rounding

A system of linear equations is a set of "recipes" that must all hold at once: 1·x₁ + 2·x₂ = 3 and 3·x₁ + 4·x₂ = 7, say. Writing the numbers in a grid A (rows = equations), the unknowns as a column x and the right-hand sides as a column b, the whole set is the single sentence A x = b. If A is nonsingular (its rows are genuinely independent, so no row is a combination of the others) there is exactly one x, namely x = A⁻¹b. If A is singular you get either no solution or infinitely many.

Small example. Take x₁ + 2x₂ = 3 and 3x₁ + 4x₂ = 7. Subtract 3 × (equation 1) from equation 2: (4 − 6)x₂ = 7 − 9, i.e. −2x₂ = −2, so x₂ = 1; put that back into equation 1: x₁ = 3 − 2 = 1. The answer is x = (1, 1). Everything in this chapter is a careful, computer-friendly way of doing exactly this — clearing entries below a diagonal and then reading off the unknowns from the bottom up.

Words used in this lesson. A matrix is a rectangular grid of numbers; a vector is a single column of numbers. The transpose Aᵀ flips a matrix over its diagonal (row i becomes column i). The identity matrix I has 1s on the diagonal and 0s elsewhere (multiplying by it changes nothing). Lower-triangular means every entry above the diagonal is 0, upper-triangular every entry below is 0, and unit lower-triangular adds “1s on the diagonal”. A permutation matrix P only reorders rows, so PA is A with its rows shuffled. Rows and columns are numbered from 0, as in Dart, so a[i][j] is row i, column j. Σ (“sigma”) means “add up”: Σ(j < i) l[i][j]·y[j] = l[i][0]·y[0] + … + l[i][i−1]·y[i−1], and it is an empty sum (= 0) when i = 0. Θ(n²) and Θ(n³) describe how the number of arithmetic steps grows with the size n of the matrix (see the growth-of-functions lesson); the nested loops used throughout are the ones from insertion sort. Matrix multiplication and its divide-and-conquer speed-ups are met again in the divide-and-conquer lesson. The Schur complement and “symmetric positive-definite” are defined where they first matter (the LU and SPD sections).

Computing A⁻¹ and multiplying is a bad way to solve a system (slower and less accurate). The trick is to notice that some systems are trivially easy: in a lower-triangular matrix (zeros above the diagonal) equation 0 has one unknown, equation 1 then has one new unknown, and so on — you solve top to bottom (forward substitution). In an upper-triangular matrix you solve bottom to top (back substitution). Each costs only Θ(n²).

LUP decomposition. Find a unit lower-triangular L (1s on the diagonal), an upper-triangular U and a permutation P (a reshuffling of rows) with PA = LU. Then Ax = b becomes PAx = Pb, i.e. L(Ux) = Pb. Call y = Ux: solve Ly = Pb by forward substitution, then Ux = y by back substitution. P is stored as an array π (called pi in code) whose values are row numbers 0 … n−1: row i of PA is row π[i] of A, so (Pb)[i] = b[π[i]].

Rounding used on this page. All computation happens in ordinary 64-bit floating point (the same as Dart's double), but every number displayed is rounded to at most 3 decimal places with trailing zeros dropped (so 0.6 and 5 stay short, and 0.214285… shows as 0.214). Nothing is rounded between steps — only the printing is. The worked matrices on this page are chosen to be exact in decimals, so they display exactly.

Reading π. The array pi says which original row sits in each row of PA: row i of PA is row pi[i] of A, and the same reordering is applied to b. Dart matrices are List<List<double>>. The two small helpers below show how pi is applied.

Every check below compares vectors and matrices with a few small helper functions; this is their complete code (plain Dart lists, no packages):

Forward substitution: solving L y = Pb

Picture a queue of people passing a number along. Person 1 can answer at once. Person 2's answer needs only person 1's. Person 3 needs the answers of persons 1 and 2, and so on. Nobody ever waits for someone behind them in the line, so you walk down the queue exactly once, top to bottom. That is forward substitution.

Small example. L = [[1, 0], [2, 1]] and Pb = (3, 8). Row 0 reads y[0] = 3. Row 1 reads 2·y[0] + y[1] = 8, so y[1] = 8 − 2·3 = 2. The answer is y = (3, 2) — each row used only values that were already found.

This is the first half of lupSolve. The panel below gives it as a procedure of its own so it can be traced alone. Because L is unit lower-triangular (1s on its diagonal) no division is needed: y[i] = b[π[i]] − Σ(j < i) l[i][j]·y[j]. Watch three things: π chooses which entry of b starts each row; the terms l[i][j]·y[j] use only y-values that are already final; and the last frame multiplies L·y to prove that it equals Pb.

The second example is the degenerate case: with L = I every product is 0 and y is just b reordered by π. The third lets you type any system; lupDecomposition is run first to supply L and π.

Starting from b instead of Pb. With the first example below (pi = [1, 2, 0], b = (9, 4, 2)) the first row of the permuted system uses b[1] = 4 (because pi[0] = 1), not b[0] = 9. Forgetting the permutation gives a perfectly smooth-looking but wrong y, and the mistake only shows up when you check A·x = b at the end.

Running time and correctness. Row i (counting from 0) does i multiplications, so the total is 0 + 1 + … + (n−1) = n(n−1)/2 multiplications and as many subtractions: Θ(n²) (n = 3: 3 multiplications; n = 100: 4,950). Invariant: at the start of iteration i, y[0] … y[i−1] already satisfy the first i equations of Ly = Pb; since row i of L has zeros beyond column i and l[i][i] = 1, equation i contains only the unknown y[i], and line 4 solves it exactly. For a general (non-unit) lower-triangular matrix you would also divide by l[i][i].

Input size → what is feasible: n ≤ 2000 → n(n−1)/2 ≈ 2·10⁶ multiplications (instant); n = 10⁵ → 5·10⁹ multiplications and a dense 10¹⁰-entry matrix (impossible), so a structured matrix (banded, see the tridiagonal question) is needed.

Forward substitution solves a lower-triangular system top to bottom in Θ(n²). Row i needs only the answers above it. Always start from (Pb)[i] = b[π[i]].

Dart implementation (it expects the permuted right-hand side, so call it as forwardSubstitution(l, [for (final i in pi) b[i]]), because pi[i] is the original row that row i of PA came from):

Back substitution: solving U x = y

Now the queue is reversed: the person at the END of the line has only one unknown and can answer immediately; the person before them needs just that one answer plus their own; and so on back to the front. You walk up the queue once, from the bottom row to the top.

Small example. U = [[2, 1], [0, 3]] and y = (5, 6). Row 1 reads 3·x[1] = 6, so x[1] = 2. Row 0 reads 2·x[0] + 1·x[1] = 5, so x[0] = (5 − 1·2)/2 = 1.5. The answer is x = (1.5, 2). This time each row divides by its diagonal entry u[i][i].

This is the second half of lupSolve: x[i] = (y[i] − Σ(j > i) u[i][j]·x[j]) / u[i][i], computed for i = n − 1 down to 0. The division is where a zero on the diagonal of U becomes fatal: it means A was singular.

Going the wrong way. If you loop i = 0 … n − 1 upwards the sum needs x-values that do not exist yet. The loop must run downwards, and the divisor must be the diagonal entry of the same row (u[i][i]), not of the row below.

Running time and correctness. Row i does n − 1 − i multiplications and one division, so the total is n(n−1)/2 multiplications + n divisions = Θ(n²) (n = 3: 3 multiplications, 3 divisions; n = 100: 4,950 and 100). Invariant: at the start of iteration i, x[i+1] … x[n−1] already solve the last n − 1 − i equations; equation i has zeros left of the diagonal in U's row, so only x[i] is new, and one division isolates it. Together with forward substitution: n(n−1) multiplications + n divisions for lupSolve.

Input size → what is feasible: n ≤ 2000 → n(n−1)/2 ≈ 2·10⁶ multiplications and 2,000 divisions (instant); n = 10⁵ → 5·10⁹ steps and 10¹⁰ stored entries (impossible for a dense U).

Back substitution solves an upper-triangular system bottom to top in Θ(n²). It divides by u_ii, so a zero diagonal entry means “singular”.

Dart implementation (0-indexed; y is the vector produced by forward substitution):

lupSolve(l, u, pi, b)

Imagine unpacking a chain of dominoes: L is a staircase where each equation only depends on answers you already have (top to bottom), so you collect y one value at a time. U is the same staircase upside-down, so you collect x from the last unknown up. π simply says "which of the original equations sits in this row after the rows were shuffled to help numerical stability".

Small example. A = [[1, 2], [3, 4]], b = (5, 11). The biggest first-column entry is 3 (row 1), so lupDecomposition swaps: pi = [1, 0], L = [[1, 0], [1/3, 1]], U = [[3, 4], [0, 2/3]]. Forward: y[0] = b[1] = 11, then y[1] = b[0] − (1/3)·11 = 4/3. Back: x[1] = (4/3)/(2/3) = 2, then x[0] = (11 − 4·2)/3 = 1. So x = (1, 2), and indeed 1 + 2·2 = 5 and 3 + 4·2 = 11.

Inputs are an already-computed L, U, π (from the sections below) and any right-hand side b. Below, the system A = [[3,6,6],[4,4,2],[3,9,4]], b = (9,4,2) is solved. Watch three things: π picks which entry of b starts each forward step; each y[i] subtracts the terms already known; and the residual check at the end.

Solving with the wrong b. lupSolve receives the ORIGINAL b and applies π itself. Passing an already-permuted b (or forgetting that L and U describe PA, not A) yields a solution to a different system. The residual check A·x = b on the last frame is the quick way to catch this.

Running time. Forward substitution multiplies l_ij·y_j for every j < i: 0 + 1 + … + (n−1) = n(n−1)/2 multiplications. Back substitution does the same plus n divisions. Total: n(n−1) multiplications + n divisions = Θ(n²) (n = 3: 6 multiplications and 3 divisions; n = 100: 9,900 and 100). This is why factoring once and re-solving for many b's is cheap. Correctness: LUx = Pb by construction, and each substitution step solves one equation exactly using values that are already final, so no step can be wrong unless a division by u_ii = 0 happens — precisely the singular case the edge example shows. The table in the inversion section counts these operations by running the loops.

Input size → what is feasible: n ≤ 2000 → n(n−1) ≈ 4·10⁶ multiplications per solve (instant), so 1,000 right-hand sides cost 4·10⁹ — still better than re-factoring each time; n = 10⁴ → 10⁸ per solve but each factor holds 10⁸ doubles = 800 MB.

lupSolve = forward substitution (Ly = Pb) then back substitution (Ux = y): Θ(n²) once PA = LU is known. Factor once, solve for as many b as you like.

Dart implementation (0-indexed; pi[i] holds a 0-based row number):

luDecomposition(a)

This is Gaussian elimination with a receipt. To clear the entries below the first pivot you subtract multiples of the top row from the rows beneath; the multiples you used are exactly the numbers in L's first column, and the top row itself is U's first row. What remains is a smaller (n−1)×(n−1) problem (the Schur complement S = A′ − v wᵀ / a₀₀), and you repeat. When everything is cleared, L (multipliers) times U (final rows) rebuilds A.

Small example. A = [[2, 1], [4, 5]]. The pivot is u[0][0] = 2 and the multiplier is l[1][0] = 4/2 = 2, which clears the 4. The Schur complement is the single number 5 − 2·1 = 3, so u[1][1] = 3. Result: L = [[1, 0], [2, 1]], U = [[2, 1], [0, 3]], and L·U = [[2, 1], [4, 5]] = A.

Precisely: split A = [[a₀₀, wᵀ], [v, A′]]. Then A = [[1, 0], [v/a₀₀, I]]·[[a₀₀, wᵀ], [0, S]] with S = A′ − v wᵀ/a₀₀. The pseudocode is this recursion with the tail call turned into the loop for k. The loop invariant: before iteration k (counting from 0), rows and columns 0 … k−1 of L and U are final, and the working array's lower-right (n−k)×(n−k) block is the current Schur complement.

Plain LU has no way to swap rows. If a pivot u[k][k] is 0, line 7 divides by zero — even when A is perfectly nonsingular (example 2). Tiny pivots are nearly as bad (they give huge multipliers and lose digits). LU is only safe for special matrices such as symmetric positive-definite or diagonally dominant ones (see the SPD section); everything else uses LUP.

Running time. The update on line 10 executes (n−1−k)² times in iteration k, so the total number of multiply-subtracts is Σ_k=0ⁿ⁻¹(n−1−k)² = (n−1)n(2n−1)/6 = Θ(n³) (≈ n³/3). Counting confirmed by the verify program:

n	2	3	4	5	10
multiply-subtracts (line 10)	1	5	14	30	285

Doubling n multiplies the work by about 8.

Input size → what is feasible: n ≤ 500 → n³/3 ≈ 4·10⁷ multiply-subtracts (instant); n = 2000 → 2.7·10⁹ (tens of seconds in Dart); n = 10⁴ → 3.3·10¹¹ (hopeless). Use LU only for SPD or diagonally dominant matrices.

luDecomposition records each elimination step: multipliers go into L, pivot rows into U, and the shrinking Schur complement is what is left to factor. It costs about n³/3 multiply-subtracts and works only when no pivot is zero.

Dart implementation (returns a record (l, u); the throw is an added guard for the zero-pivot case that the pseudocode leaves out):

lupDecomposition(a)

Before each elimination step, look down the current column and pick the row whose entry is biggest in absolute value — like choosing the sturdiest plank to stand on — then swap that row to the top. Dividing by a big pivot keeps every multiplier between −1 and 1, so rounding errors cannot blow up; and if the biggest entry is 0 the whole column is empty, which proves A is singular. That is partial pivoting.

Small example. A = [[0, 1], [1, 1]]. Plain LU would divide by a[0][0] = 0. Partial pivoting scans column 0 = (0, 1), finds the larger |entry| in row 1 and swaps the rows (pi = [1, 0]): the matrix becomes [[1, 1], [0, 1]], the multiplier is 0/1 = 0, and the Schur complement is 1 − 0·1 = 1. So L = I and U = [[1, 1], [0, 1]], with PA = LU.

The procedure works in place: when it finishes, the array a holds L's multipliers below the diagonal (L's unit diagonal is implied) and U on and above the diagonal, and pi records the row order. Note that line 12 swaps whole rows, including the multipliers already stored to the left, so earlier steps stay consistent with the new order.

Read every animation the same way: first the pivot choice (scan the column, remember the biggest absolute value), then the row swap (π and whole rows move together), then the multipliers and the update, and finally a frame where the lower-right block is highlighted as the Schur complement — the smaller matrix of the same kind that the next iteration factors. The Schur complement of the pivot a[k][k] is S = A′ − v·wᵀ / a[k][k], where v is the column below the pivot, wᵀ the row to its right and A′ the remaining block; its entries are exactly the numbers written by lines 16–17.

Biggest number is not biggest size. The pivot is the entry of largest absolute value: in the column (2, −9, 4) the pivot is −9, not 4. Comparing signed values (or swapping only part of a row) breaks the guarantee |l_ij| ≤ 1 and, in the second case, breaks PA = LU itself.

Running time. Pivot search is Θ(n²) in total, swaps are Θ(n²), and the update loop is the same (n−1−k)² per step as plain LU: Θ(n³). Pivoting costs only a constant factor extra. Correctness: PA = LU follows from the same Schur-complement argument, applied after permuting the pivot row to the top; nonsingularity of A implies the Schur complement is nonsingular, so a nonzero pivot always exists until the end. Partial pivoting also guarantees |l_ij| ≤ 1 (checked on 300 random matrices in the verify program).

Input size → what is feasible: n ≤ 500 → n³/3 ≈ 4·10⁷ multiply-subtracts plus n(n+1)/2 = 125,250 pivot comparisons (instant); n = 2000 → 2.7·10⁹ (tens of seconds); n ≥ 10⁴ → use a library (BLAS/LAPACK) or exploit structure.

lupDecomposition = elimination that always brings the largest-magnitude entry of the column to the diagonal first. It never divides by zero on a nonsingular matrix, keeps every multiplier within [−1, 1], and costs Θ(n³). It reports “singular matrix” when a whole pivot column is zero.

Dart implementation (the row-swap step becomes a swap of whole row references; the singular check uses the tolerance eps = 1e-12 instead of exact zero because floating-point subtraction rarely gives exactly 0. The animations on this page use the scale-free version of the same test, pivot ≤ 10⁻¹⁰·max|a_ij|, so that a singular matrix with entries near 10⁶ is still caught; the fixed eps is fine for the modest entries used in the Dart examples):

Why pivoting matters: zero and tiny pivots

Dividing by a tiny number is like looking at a rounding error through a powerful magnifying glass. A computer stores only about 16 significant digits, so every number carries a hidden error of roughly one part in 10¹⁶. Divide by 10⁻²⁰ and that invisible error is multiplied by 10²⁰ — suddenly it is bigger than the answer. Dividing by exactly 0 is the extreme case where nothing works at all. Choosing a large pivot means dividing by a big number, which shrinks errors instead.

Small example. A = [[0, 1], [1, 1]] is invertible (det = −1), yet its first pivot is 0 and plain elimination stops immediately (the previous section's pitfall). Exchanging the two rows fixes it. Now imagine the 0 is replaced by 10⁻²⁰: the matrix is still fine mathematically, but the computer now divides by a number that is almost zero.

The system below is A = [[e, 1], [1, 1]] with right-hand side b. Its exact solution has a closed form (x[0] = (b[0] − b[1])/(e − 1), x[1] = b[1] − x[0]), which serves as an unrounded reference. Each player runs the same five arithmetic lines twice on real 64-bit doubles: once with the row exchange switched off, once with partial pivoting on. Only the numbers differ.

“Not zero” is not the same as “safe”. Any check that only avoids exact zero (if (pivot != 0)) still lets tiny pivots through, and tiny pivots destroy accuracy long before they reach zero. The rule is to pick the largest candidate, which is why lupDecomposition compares absolute values instead of testing for zero. (Pivoting protects against unstable elimination; it cannot rescue a matrix that is itself extremely sensitive — a different issue called ill-conditioning.)

Why it works. With partial pivoting every multiplier satisfies |l_ij| ≤ 1, so a step can add at most the size of the pivot row to another row: the entries of the Schur complements cannot grow by more than a factor of 2 per step. In theory the growth could reach 2ⁿ⁻¹ for contrived matrices, but such cases are known to be extremely rare in practice, which is why partial pivoting is the standard default. Without pivoting there is no bound at all: the multiplier is a[1][0]/a[0][0] and can be arbitrarily large, as the first example shows (l[1][0] = 10²⁰). For special matrices — symmetric positive-definite (see the SPD section) or diagonally dominant — the pivots are automatically safe and pivoting can be skipped.

Input size → what is feasible: the 2×2 experiment is O(1); on an n × n matrix pivoting adds only n(n+1)/2 comparisons and at most n row swaps to the Θ(n³) work, so n ≤ 500 stays instant — there is never a reason to skip pivoting on a general matrix.

A zero pivot stops elimination; a tiny pivot silently ruins it. Partial pivoting (largest absolute value in the column, rows swapped) fixes both at a cost of a few comparisons and Θ(n²) swaps.

Dart version of the two-pass experiment (the same operations, in the same order, as the animation; the verify program checks the numbers quoted above):

Inverting a matrix with LUP

A⁻¹ is the matrix X with AX = I. Read that column by column: A times X's j-th column must equal the j-th column of the identity, e_j. That is n separate linear systems with the same A — so factor A once and call lupSolve n times.

Small example. A = [[1, 2], [3, 4]] has the factors from the lupSolve example. Column 0 of A⁻¹ solves Ax = (1, 0): forward gives y = (0, 1), back gives x[1] = 1/(2/3) = 1.5 and x[0] = (0 − 4·1.5)/3 = −2, so column 0 is (−2, 1.5). Column 1 solves Ax = (0, 1) and gives (1, −0.5). Hence A⁻¹ = [[−2, 1], [1.5, −0.5]].

The panel below writes the method as steps. Cost: one Θ(n³) decomposition plus n solves of Θ(n²) each: Θ(n³) overall. (Inversion and matrix multiplication have the same asymptotic difficulty: multiplying reduces to inverting a block matrix, and inverting reduces to multiplying via the Schur complement — both reductions are coded and checked in the panel below.)

Counting the work. Each number below comes from running the same loops as the algorithms (instrumented counters), and the verify program repeats the counts and checks them against the closed forms n(n−1) and (n−1)n(2n−1)/6. The inversion column is one decomposition plus n solves. At n = 100 a single solve costs 9,900 multiplications, while a full inversion costs about 1.3 million — roughly 133 times more.

Inverting is almost never needed. To solve Ax = b, lupSolve is faster (Θ(n²) after factoring, no extra n solves) and more accurate than computing A⁻¹b. Use an explicit inverse only when you truly need every entry of A⁻¹.

A⁻¹ is n solves of A x = e_j that share one factorization: one Θ(n³) decomposition plus n Θ(n²) solves = Θ(n³). Compute it only when you truly need every entry.

Input size → what is feasible: n ≤ 300 → (4/3)·n³ ≈ 3.6·10⁷ multiplications (instant); n = 1000 → 1.3·10⁹ (about 15 s); a single system Ax = b never needs it — lupSolve costs one factorization plus n(n−1).

Dart implementation (with the small helper that splits the in-place result into L and U):

Least squares and the normal equations

You measured 5 points that almost lie on a smooth curve, but noise means no curve passes through all of them. Least squares picks the curve whose vertical misses ε_i = F(x_i) − y_i have the smallest total squared size ‖ε‖². Squaring makes big misses count much more and makes the minimum a clean calculus problem.

Small example. Points (0, 1), (1, 3), (2, 2) and a line F(x) = c[0] + c[1]·x. Then A = [[1,0],[1,1],[1,2]], AᵀA = [[3, 3], [3, 5]] and Aᵀy = (6, 7). Solving [[3,3],[3,5]]·c = (6,7) gives c = (1.5, 0.5), the line 1.5 + 0.5x. Its misses are ε = (0.5, −1, 0.5): they add up to 0 and are perpendicular to both columns of A, exactly what the normal equations demand.

Model F(x) = c[0] + c[1]·x + … + c[n−1]·xⁿ⁻¹. Stack the m data points into the m × n matrix A with A[i][j] = x[i]^j, so the prediction vector is Ac and the miss vector is ε = Ac − y (m > n, so this is overdetermined). Setting the derivatives of ‖ε‖² to zero gives the normal equations AᵀA c = Aᵀy. If A has full column rank, AᵀA is symmetric positive-definite (SPD), hence nonsingular, and c = (AᵀA)⁻¹Aᵀy = A⁺y where A⁺ is the pseudoinverse. The panel below is this lesson’s own step-by-step summary of the method. In practice you do not form the inverse: multiply out AᵀA and Aᵀy, decompose the n × n system and call lupSolve. (For SPD matrices plain LU never meets a zero pivot, as the SPD section shows, but LUP is harmless and safer, so the panel uses it.)

The last frame of each example also draws the result: orange dots are the data, the blue curve is the fitted F, and red vertical segments are the misses ε (zero-length when the fit is exact).

More coefficients is not better. With n = m the curve passes through every point (ε = 0) but “fits the noise”; and with fewer distinct x values than coefficients the columns of A are dependent, AᵀA is singular and the fit is not unique (custom input shows the “singular matrix” error). Choose n well below m.

Running time. Forming AᵀA costs Θ(mn²), Aᵀy costs Θ(mn), the decomposition Θ(n³) and the solve Θ(n²): Θ(mn² + n³) — for fixed low degree, linear in the number of data points. Correctness: the minimum of a convex quadratic is where its gradient vanishes; Aᵀ(Ac − y) = 0 says the error vector is perpendicular to every column of A, which the final frame of each example verifies numerically. Warning: forming AᵀA squares the condition number, so for badly-conditioned fits professionals use QR factorization; the normal equations are the simple route.

Input size → what is feasible: m ≤ 10⁵ points with n ≤ 10 coefficients → m·n² ≈ 10⁷ multiply-adds (instant); n = 1000 coefficients → n³ = 10⁹, and AᵀA squares the condition number, so use QR instead.

Least squares = build A, form AᵀA c = Aᵀy (an n × n system that is symmetric positive-definite), and hand it to lupDecomposition + lupSolve. The error vector ends up perpendicular to every column of A.

Dart implementation (needs import 'dart:math' as math;; the powers use math.pow(x, j) for x^j):

Symmetric positive-definite matrices

Imagine a perfectly smooth bowl. Wherever you stand and whichever direction you walk away from the very bottom, the ground rises. A symmetric positive-definite (SPD) matrix is the algebraic version of such a bowl: the quantity xᵀAx (a single number for every direction x) is strictly positive for every x ≠ 0. A matrix whose “bowl” has a saddle or a flat valley has some direction in which the number is zero or negative, and it is not positive-definite.

Two words to define. A matrix is symmetric if A = Aᵀ (its entries mirror across the diagonal: a[i][j] = a[j][i]). It is positive-definite if xᵀAx > 0 for every nonzero vector x. SPD matrices are the ones you meet in least squares (AᵀA with independent columns), in physics (energy) and statistics (covariance). Four facts about them: they are nonsingular, every leading submatrix is SPD, the Schur complement of an SPD matrix is SPD, and therefore luDecomposition never divides by zero on them — every pivot is strictly positive.

Small example. A = [[2, 1], [1, 2]]. The first pivot is 2 > 0; the Schur complement is 2 − 1·1/2 = 1.5, so the second pivot is 1.5 > 0. Both pivots are positive, so A is SPD, and det A = 2·1.5 = 3. Contrast [[1, 2], [2, 1]]: its Schur complement is 1 − 2·2/1 = −3 < 0, so it is symmetric but not positive-definite. The players below trace exactly this test: symmetric check, then pivot, column v, and the new Schur complement, step by step.

The procedure in the code panel turns those facts into a test. It uses the fact that a symmetric matrix is positive-definite exactly when plain LU produces only positive pivots, and it never needs row swaps — which is why LU (and even the square-root variant, the Cholesky factorization A = LLᵀ, whose diagonal is the square roots of the pivots) is safe here. The least-squares matrix AᵀA of the previous section is SPD whenever A has independent columns, which is why the normal equations can always be solved.

Symmetric is not enough. [[1, 2], [2, 1]] is perfectly symmetric but has a negative second pivot; [[1, 1], [1, 1]] is symmetric with a zero second pivot (singular). Also note the test is on the pivots of the Schur complements, not on the entries of A — A = [[4, 5], [5, 4]] has all positive entries and still fails.

Why the Schur complement stays SPD. Write A = [[A_k, Bᵀ], [B, C]] and S = C − B·A_k⁻¹·Bᵀ. For any x = (y, z), xᵀAx = (y + A_k⁻¹Bᵀz)ᵀ A_k (y + A_k⁻¹Bᵀz) + zᵀSz (“completing the square”). Choosing y = −A_k⁻¹Bᵀz kills the first term, leaving zᵀSz = xᵀAx > 0 for every z ≠ 0, so S is positive-definite. By induction every pivot is positive. Cost: the test is one LU without pivoting, ≈ n³/3 multiply-subtracts, versus ≈ n³/6 for Cholesky; the check of symmetry is Θ(n²).

Input size → what is feasible: n ≤ 500 → the LU pivot test is n³/3 ≈ 4·10⁷ steps and Cholesky (n−1)n(n+1)/6 ≈ 2·10⁷ (both instant); n = 5000 → 4·10¹⁰ / 2·10¹⁰ (minutes), so use sparse or iterative methods for huge SPD systems.

SPD = symmetric with xᵀAx > 0 for all x ≠ 0. Equivalent test: symmetric and every LU pivot strictly positive. SPD systems need no pivoting, and AᵀA (independent columns) is always SPD.

Dart implementation of the test (same idea, run on a copy of A; eps is the same 1e-12 tolerance used for singular pivots):

Quiz

Interview questions

Cheat sheet

Procedure	Time	Space	Needs pivoting?	Fails when	When to use
Forward substitution (Ly = Pb)	Θ(n²): n(n−1)/2 mult	O(n)	n/a	never (unit diagonal)	Lower-triangular systems; first half of lupSolve
Back substitution (Ux = y)	Θ(n²): n(n−1)/2 mult + n div	O(n)	n/a	u_ii = 0 (singular)	Upper-triangular systems; second half of lupSolve
lupSolve	Θ(n²)	O(n) (x, y)	n/a (uses π)	u_ii = 0 (singular)	Every solve after a factorization; many right-hand sides
luDecomposition	Θ(n³) ≈ n³/3 mult-subs	Θ(n²) (L, U)	No (cannot swap)	Any zero pivot, even if A nonsingular	SPD or diagonally dominant matrices
lupDecomposition	Θ(n³)	Θ(n²) in place + π	Yes (partial)	Whole pivot column is 0 → "singular matrix"	Default for general square A
Inversion via LUP	Θ(n³)	Θ(n²)	Yes	Singular A	Only when all of A⁻¹ is needed
Pivoting (partial)	O(n²) extra comparisons and swaps	O(n) for π	—	whole column zero → singular	Always for general matrices; keeps \|l_ij\| ≤ 1
Least squares (normal equations)	Θ(mn² + n³)	Θ(n² + mn)	Optional (AᵀA is SPD)	Columns of A dependent (AᵀA singular)	Curve fitting with few parameters, well-conditioned data
SPD test / Cholesky	LU test ≈ n³/3, Cholesky ≈ n³/6; symmetry check Θ(n²)	Θ(n²)	No (pivots stay positive)	Nonpositive pivot → not positive-definite	Covariance / AᵀA / energy matrices; fastest factorization for SPD