lin

Linear algebra — machines/lin.shoddy

Summary

lin is everything that asks a question of a matrix — a grid of numbers in rows and columns — or takes one apart:

determinants — the one number that sums up a square matrix — and inverses
linear systems: sets of equations solved together, like A x = b
echelon forms and rank
the LU, QR and Cholesky factorisations — ways of splitting a matrix into simpler pieces whose product rebuilds it
characteristic polynomials and eigenvalues — an eigenvalue is a number that says how much the matrix stretches things, and the characteristic polynomial is the equation whose roots they are
the vector words that are geometry rather than arithmetic

Every word is pure and every word is prefixed Lin.

The type lives in matrix, not here. Matrix, Mat, MatGet, MatMul, Transp and Ident are that machine's. A program that uses lin normally includes both, because it has to name Matrix to call anything on this page:

Include "matrix.shoddy"
Include "lin.shoddy"

A Brief History of Elimination

The oldest algorithm in this machine is also the central one, and it is about two thousand years old. Chapter eight of the Nine Chapters on the Mathematical Art, compiled in Han-dynasty China, solves systems of linear equations by laying the coefficients out as a rectangular array of counting rods. Then it does exactly what LinSolve does: it scales and subtracts rows until the answers fall out. Europe reinvented the method piecemeal. It picked up Gauss's name because he systematised it for the least-squares work that recovered the lost asteroid Ceres in 1801 — a naming so unhistorical that Gauss himself called the method ordinary. The modern chapter came in 1948, when Alan Turing framed elimination as a factorisation — the L and U this machine computes. He also analysed how its rounding errors behave, founding the numerical linear algebra that decides details like pivoting. Rods to Ceres to rounding analysis, it has always been the same move: make the rectangle triangular, then read the answers off the edge.

Why It's Useful

Why it is a separate machine. matrix owns a data structure and its arithmetic. lin owns the algorithms over it. Split that way, a program that only multiplies matrices does not carry an eigenvalue solver, and an algorithm added here never disturbs the type. It is the same line eng draws against stats, and the same one json draws against dict.

Partial pivoting throughout. A pivot is the entry an elimination step divides by. Every elimination here takes the largest available entry in the column as its pivot. Without that, a perfectly ordinary matrix divides by a number that is only small by accident, and the answer is noise that still looks like an answer.

Every factorisation returns one record. LinLuOf hands back LFactor, UFactor, Perm and PSign together, because P A = L U is only true of all four at once. Two words returning two halves would invite a caller to keep one and drop the other.

Solvers refuse rather than mislead — the same rule eng follows. LinSolve checks the determinant before it starts eliminating. A singular system — one whose determinant is zero — has either no solution or infinitely many. Handing back one vector from either case, without saying which, is worse than stopping.

User's Guide

Solving a system

Include "matrix.shoddy"
Include "lin.shoddy"

Def Main()
    ' x + 2y = 5, 3x + 4y = 11
    Let a = Mat(2, 2, { 1, 2, 3, 4 })
    Let b = MatRow(Mat(1, 2, { 5, 11 }), 1)

    Let x = LinSolve(a, b)
    Print(Nth(x, 1))          ' 1
    Print(Nth(x, 2))          ' 2

    Print(LinDet(a))          ' -2
    Print(MatGet(LinInv(a), 1, 1))   ' -2

LinSolve is Gauss-Jordan on the augmented matrix, not LinInv followed by a multiply. Inverting to solve is both slower and less accurate, and it is the habit this word exists to remove. LinInv is still here for when the inverse itself is the answer.

Rank and echelon forms

Let s = Mat(3, 3, { 1, 2, 3, 2, 4, 6, 1, 1, 1 })   ' row 2 is twice row 1
Print(LinRank(s))                  ' 2
Print(LinDet(s))                   ' 0
Print(LinRank(Mat(2, 3, { 1, 2, 3, 2, 4, 6 })))    ' 1 — works on any shape

LinRef gives row echelon form and LinRref the reduced form, where every pivot is 1 and alone in its column. Pivot columns are skipped rather than forced, so both work on any shape and any rank. And LinRank, LinInv and LinSolve are all LinRref underneath.

Factorisations

Let f = LinLuOf(a)
' P A = L U, and the pivots times the sign are the determinant
Print(MatGet(LFactor(f), 1, 1))    ' 1 — Doolittle puts the unit diagonal on L
Print(PSign(f))                    ' -1 — a's pivot was in the second row

Let q = LinQrOf(a)                 ' A = Q R, Q orthonormal, R upper triangular
Let l = LinChol(b)                 ' A = L L', for symmetric positive definite

There is more than one valid L and U, but only one valid product. That is how a factorisation is testable at all: the suite checks P A - L U, Q R - A and L L' - A rather than any individual cell.

LinChol puts its refusal exactly where positive definiteness fails — the square root of a negative pivot. That is a fact about the matrix, not a rounding accident, and no other test is needed.

Eigenvalues

' The 2,-1 tridiagonal: eigenvalues 2-sqrt2, 2, 2+sqrt2
Let t = Mat(3, 3, { 2, 0 - 1, 0, 0 - 1, 2, 0 - 1, 0, 0 - 1, 2 })

Print(LinCharPoly(t))              ' { 1, -6, 10, -4 }
Print(Nth(LinEigVals(t), 1))       ' 0.58578...

Let e = LinEigSym(t)               ' values AND vectors, paired by position
Print(Nth(Values(e), 1))
Print(MatCol(Vectors(e), 1))       ' the eigenvector that goes with it

LinCharPoly is Faddeev–LeVerrier. It gets every coefficient from traces of matrix products and never expands a determinant symbolically. That is what makes a characteristic polynomial computable here at all. Its coefficients come back highest power first, matching eng, so they feed straight into EngPolyEval or EngCubicReal without reversing anything.

LinEigSym is cyclic Jacobi and is the complete answer for a symmetric matrix: a full set of real eigenvalues, with orthonormal eigenvectors, paired by position. Prefer it whenever the matrix is symmetric.

What It Will Not Do

LinEigVals finds the real eigenvalues only. It brackets sign changes of the characteristic polynomial, and a complex pair produces no sign change — nor does a repeated root of even multiplicity. So the length of the answer carries information. Length(LinEigVals(m)) short of Rows(m) means some eigenvalues are not real, or not simple, and were not found. A 90° rotation matrix returns an empty list, which is the honest answer: nothing it acts on keeps its direction.
n ≈ 8. Every word here is written for clarity over speed, and Shoddy is pure, so each elimination step builds a whole new matrix. LinEigVals and LinEigSym iterate and feel it first.
Cancellation in LinCharPoly. A characteristic polynomial's coefficients span magnitudes that subtraction eats, and everything built on it inherits that.
No singular value decomposition, and no least squares. LinQrOf is the piece those would be built from, and it is here.

Word Reference

Every word, and the records they hand back

The records

Word	Description
LinLu	`LFactor`, `UFactor`, `Perm`, `PSign` — the four parts of `P A = L U`. Spelled out rather than L/U/P/S because a one-letter accessor is shadowed by any local of the same letter.
LinQr	`QFactor`, `RFactor` — for the same reason.
LinEigen	`Values`, `Vectors`. The vectors are the columns, so the k-th value pairs with `MatCol(Vectors, k)`.
LinJac	`Work`, `Basis` — the Jacobi sweep's accumulator, exposed because the machine is one flat surface.

Shape and measure

Word	Description
LinIsSquare(m)	Rows equal columns.
LinIsSym(m, tol)	Symmetric to a tolerance, never to the bit: a matrix assembled by arithmetic is symmetric in intent and not in the last few digits.
LinTrace(m)	The diagonal sum.
LinNorm(m)	The Frobenius norm — the length of the matrix read as one long vector.
LinMaxAbs(m)	The largest magnitude in it.

Vectors

Word	Description
LinVNorm(v)	Euclidean length of an `Array Of Number`.
LinVUnit(v)	The same direction, length 1. Errors on the zero vector, which has no direction.
LinCross(u, v)	The cross product, three-dimensional and nothing else. A version that quietly padded to three would be answering a different question.

Building and shuffling

Word	Description
LinBuild(rows, cols, f)	Give it dimensions and a rule for one cell — `[ Number Number -- Number ]` — and it makes the matrix. Everything below is written in terms of it.
LinSwapRows(m, i, j)	Exchange two rows.
LinScaleRow(m, p, k) / LinScaleCol(m, p, k)	Multiply one row or column.
LinAug(a, b)	Side by side. Errors when the row counts differ.
LinCols(m, lo, hi)	Columns lo through hi — the half of an augmented matrix that carries the answer.

Elimination

Word	Description
LinRef(m)	Row echelon form: zeros below each pivot, pivots left to right, zero rows at the bottom.
LinRref(m)	Reduced row echelon form: every pivot 1 and alone in its column. Rank, inverse and solution all fall out of this one.
LinRank(m)	The number of non-zero rows in the reduced form. Works on any shape.
LinPivot(m, col, from) / LinPivotGo(...)	The partial-pivot search: the row at or below `from` with the largest magnitude in that column.
LinElimBelow(m, p, col) / LinElimAll(m, p, col)	One elimination step, below the pivot row or on every other row.

Determinant and inverse

Word	Description
LinDet(m)	By elimination, not cofactor expansion: cofactors are n! multiplications and elimination is n³, which at n = 10 is the difference between an answer and a hung program.
LinInv(m)	Gauss-Jordan on `[A \| I]`. Errors on a singular matrix rather than returning something.
LinSolve(a, b)	Solves `A x = b` for an `Array Of Number` and returns one. Checks the determinant first.
LinPow(m, n)	Square-and-multiply, so a large exponent costs log n multiplications. A negative exponent inverts first.
LinMinor(m, i, j)	The matrix with row i and column j struck out.
LinCofactor(m, i, j)	The signed minor determinant.
LinAdj(m)	The adjugate — the transpose of the cofactor matrix. `A` times it is `LinDet(A)` times the identity.

Factorisations

Word	Description
LinLuOf(a)	Doolittle LU with partial pivoting → `LinLu`. `P A = L U`, L carries the unit diagonal, U carries the pivots, and `PSign` is the determinant of `P` — exactly the factor a determinant would otherwise lose to the row swaps.
LinQrOf(a)	Gram-Schmidt QR → `LinQr`. `A = Q R` with Q's columns orthonormal. Needs at least as many rows as columns, and refuses on linearly dependent ones.
LinChol(a)	Cholesky → the lower triangular `L` with `A = L L'`. Needs a symmetric, positive definite matrix and says which condition failed.
LinDotPart(l, u, i, j, upto)	The partial sum both LU and Cholesky are built from.

Eigenvalues

Word	Description
LinCharPoly(m)	Faddeev-LeVerrier. Highest power first, so a 3×3 gives four numbers beginning with 1 — the same convention eng's polynomial words use.
LinRootBound(cs)	Cauchy's bound: every root lies inside this radius. Generous, which is what a scan wants.
LinEigVals(m)	The real eigenvalues, ascending, by bracketing sign changes. Complex pairs and repeated roots of even multiplicity are not found, so a short answer is itself information.
LinEigSym(m)	Cyclic Jacobi → `LinEigen`. The complete answer for a symmetric matrix: every eigenvalue, with orthonormal eigenvectors, paired by position. Errors on a matrix that is not symmetric.
LinRot(n, p, q, cs, sn)	A Givens rotation, the step Jacobi is made of.
LinOffMass(m)	The total off-diagonal magnitude — what the sweep drives to zero and how it knows it has finished.

Who Uses It

No machine and no mill includes it yet — its words are already at the reckoner's prompt through its seed, and a mill that puts them to work will appear here.

The Machines It Uses

	Machine	Why
	eng	`EngSum` and `EngProdOf` for the reductions, `EngPolyEval` and `EngZeroOf` under `LinEigVals`, and `EngAngleOf` for the Jacobi rotation angle.
	matrix	The `Matrix` type itself, and all of its arithmetic — `Mat`, `MatGet`, `MatMul`, `Transp`, `Ident`, `Dot` and the vector words.
	seq	`Append` and `DropN` build the coefficient and eigenvalue lists.