Matrices, linear maps, and linear systems

Scaling horizontally and then shearing an image can give a different result from shearing first and then scaling. The same issue appears in an AI layer: a table of feature rows, a weight matrix, and a bias must meet a specific shape contract. Writing a multiplication sign without stating which map acts first or which axis holds features is insufficient to define the calculation.

This module connects matrix entries to maps, equations, and executable operations. You will derive products from composition, preserve solution sets during elimination, classify three kinds of system, and separate invertibility from numerical reliability. The final section reconciles column-vector mathematics with batched row-oriented code.

Retrieval check: compute a dot product and a linear combination, distinguish a point from a displacement, and predict a broadcast shape. Review Module 09 and its NumPy preparation. Optional algorithm-cost comparisons use Module 06; they are not required for the algebraic core. Matrices here have finite real entries and positive row and column counts.

An m×n matrix A is a rectangular array of entries aᵢⱼ, indexed by row i and column j. Its shape records the order: m output coordinates, n input coordinates when used as a column-vector map. A matrix can also store a table of m observations and n features. The same rectangular layout supports both roles, but an application must identify which interpretation is intended. Rows of a data table are not automatically output coordinates of a particular physical transformation.

For a column vector x∈Rⁿ, define y=Ax∈Rᵐ by $y_i=\sum_{j=1}^n a_{ij}x_j$. Every output entry is the dot product of a row of A with x. Matching n input entries to n row entries is necessary; m may differ from n. A two-by-three matrix maps three coordinates to two, while a three-by-two matrix has the opposite input/output sizes. Neither is interchangeable with the other merely because both store six numbers.

Alternatively Ax is a linear combination of A’s columns: $Ax=\sum_jx_jA_{:j}$. In particular Aeⱼ is column j, because the coordinate unit vector selects that column. Thus choosing the images of all coordinate unit vectors determines the map on every coordinate vector. This gives a geometric meaning to the stored columns without requiring a separate example for every possible input.

A map L is linear when L(u+v)=L(u)+L(v) and L(αu)=αL(u) for all inputs and real α. Matrix maps satisfy these laws by distributing each finite row sum. A linear map necessarily sends zero to zero: L(0)=L(0)+L(0) implies L(0)=0. Checking only a few vectors does not prove linearity, but a formula constructed as Ax has a general coordinate proof.

Conversely, any linear map from Rⁿ to Rᵐ has a matrix in these coordinate conventions. Write x=∑xⱼeⱼ and use linearity to obtain L(x)=∑xⱼL(eⱼ); choose those images as columns. Coordinate choice matters: a different basis can give a different array for the same abstract map, a topic developed in Module 11. The matrix’s relationship to its coordinates should remain explicit.

An affine map x↦Ax+b adds a fixed output offset. It preserves line and affine-combination structure but is linear only when b=0. For example x↦2x+1 on the real line fails the zero-to-zero requirement. In an AI layer the bias supplies this offset; calling every layer linear is common shorthand, but the mathematical distinction affects proofs about addition, origins, and fixed points.

Matrix addition and scalar multiplication are entrywise and require equal shapes. A+B represents the sum of two compatible linear maps with the same domain and codomain; αA represents output scaling. Addition of a two-by-three matrix and a three-by-two matrix is not defined by this rule. Array broadcasting can define a different numerical operation, so distinguish matrix algebra’s shape rules from every operation allowed by a general array library.

Let B be n×p and A be m×n. Acting on x∈Rᵖ, first B produces n coordinates and then A produces m. Their composition is represented by the m×p matrix AB with $(AB)_{ik}=\sum_{j=1}^n a_{ij}b_{jk}$. The matching inner dimension n is the number of intermediate coordinates. The outer dimensions m,p determine the output matrix shape. On column vectors, the rightmost map acts first.

Derive the formula rather than memorising a diagram. The ith coordinate of A(Bx) is ∑ⱼaᵢⱼ∑ₖbⱼₖxₖ. Reordering finite sums gives ∑ₖ(∑ⱼaᵢⱼbⱼₖ)xₖ, exactly (AB)x. This also explains why matrix multiplication differs from entrywise multiplication: intermediate coordinates are summed to form one composite coefficient. An entrywise product is useful in other tasks, but it does not implement this composition.

Compatible matrix multiplication is associative: (AB)C=A(BC). Expanding either entry gives the same finite double sum ∑ⱼ∑ₖaᵢⱼbⱼₖcₖℓ. This allows grouping a fixed ordered composition in different ways. It does not allow swapping the order. Grouping can change intermediate array sizes and numerical rounding or work, even though the real-number algebraic result is the same.

The identity matrix Iₙ has diagonal ones and all other entries zero. It maps each n-coordinate input to itself. For A of shape m×n, IₘA=A and AIₙ=A, using differently sized identities on the two sides. A diagonal matrix scales each coordinate separately. Zero diagonal entries collapse directions, whereas negative diagonal entries reverse them. Neither identity nor diagonal terminology implies a general matrix is invertible.

The transpose Aᵀ interchanges rows and columns, with (Aᵀ)ⱼᵢ=aᵢⱼ. Its shape is n×m. The product law is (AB)ᵀ=BᵀAᵀ: entry (k,i) on the left is ∑ⱼaᵢⱼbⱼₖ, while the right gives ∑ⱼbⱼₖaᵢⱼ, equal by scalar commutativity. The reversed order also makes the dimensions fit. Transposing each factor without reversing them may be undefined or represent a different map.

A symmetric square matrix satisfies A=Aᵀ. This describes entry symmetry, not automatically positivity, invertibility, or a particular data interpretation. The identity is symmetric and invertible; the all-zero square matrix is symmetric and singular. Later modules use symmetry for inner products, spectral geometry, and covariance. Here it is one checkable matrix property, with its own definition and scope.

Block matrices partition rows and columns into groups. Block products follow the same sum-of-compatible-products rule, provided intermediate partitions align. For instance [A B] times a vertically stacked vector [u;v] equals Au+Bv when u and v match the respective column groups. A block is still a rectangular map with dimensions; treating every block as a scalar can conceal an invalid product. The notation helps organise a larger system without changing the underlying entry formula.

A linear system Ax=b specifies m equations in n unknown coordinates. Row i states ∑ⱼaᵢⱼxⱼ=bᵢ. The right-hand side b has one entry per equation. The augmented matrix [A|b] stores both coefficients and constants; the separator distinguishes unknown coefficients from the fixed target. Changing only A while leaving b untouched can change the solution set, even when the coefficient operation looks algebraically familiar.

There are three elementary row operations: swap two equations; multiply an equation by a nonzero scalar; add a scalar multiple of another equation to one equation. Apply each to the entire augmented row, including the right-hand side. Each operation has an inverse of the same type, so it preserves the solution set in both directions. Multiplication by zero is excluded because it can erase a constraint and cannot be reversed.

To prove the row-addition operation, suppose an assignment satisfies both original rows. It also satisfies their indicated linear combination, so it satisfies the new row. Conversely, the unchanged source row and the new row recover the old target by subtracting the same multiple. This reverse implication is what establishes equality of solution sets. A one-way statement that all old solutions remain would not rule out extra new solutions.

A pivot is a selected nonzero coefficient used to eliminate a variable from other rows. If a candidate diagonal entry is zero but a lower row has a nonzero entry in that column, swap rows. If every remaining entry in that column is zero, skip it and continue to later columns rather than divide by zero. Pivot columns identify constrained coordinates; a skipped column can correspond to a free variable once consistency is established.

For exact arithmetic any nonzero suitable pivot permits a legal algebraic step. For floating-point work, choosing the largest available absolute entry in the column is partial pivoting, which avoids dividing by a needlessly small available pivot. This practical rule does not make every system well conditioned or prove every numerical error is small. The lab uses exact rational Fractions for its trace, separating algebraic row equivalence from floating rounding.

Row operations change equations, not the coordinate identities of the unknowns. Swapping rows reorders constraints; swapping columns reorders variables and requires corresponding changes to the interpretation of x. A solver may use column permutations deliberately, but a learner must record them. Confusing a row interchange with an unknown-variable interchange can produce a numerically plausible answer assigned to the wrong feature or quantity.

The elimination invariant is that the current augmented system has exactly the original solution set. At each step the row-operation proof preserves that fact, while the pattern of zeros makes subsequent solution easier. Termination comes from advancing the pivot row or column through finite bounds. As with earlier algorithm proofs, the invariant alone does not establish that progress occurs; the finite traversal of candidate positions supplies that argument.

An augmented system is in row-echelon form when nonzero coefficient rows appear before zero coefficient rows and each successive pivot lies strictly to the right of the preceding pivot. Entries below pivots are zero. A contradiction row may have all zero coefficients and a nonzero right-hand side; it should be recognised explicitly rather than counted as an ordinary unknown-variable pivot. We separate coefficient pivots from the augmented constants.

Gaussian elimination produces such a form by choosing a pivot, optionally normalising its row, eliminating entries below it, and continuing with the remaining rows and columns. Normalising a pivot is not required for echelon form, but the lab does so to simplify its back substitution. Eliminating above pivots as well would produce a reduced row-echelon form; that is a further transformation, not necessary merely to solve a unique triangular system.

Back substitution starts at the last pivot row. After free-variable values are chosen, each pivot equation solves its pivot variable from already known later variables, dividing by the nonzero pivot when the row is not normalised. Proceed upward until all pivot variables are assigned. The structure ensures later unknowns are known at the required moment. Solving top-down without that structure could require a value not yet determined.

Three cases describe exact real solution sets. If a row reads 0=c with c≠0, the system is inconsistent and has no solution. If no contradiction exists and every unknown column has a pivot, there is a unique solution. If no contradiction exists and at least one unknown is free, the system is underdetermined and has infinitely many real solutions, because each real choice of a free coordinate can be extended by back substitution.

For rows [1,1|3],[2,2|7], eliminating gives a contradictory zero-coefficient row with constant one, so no solution exists. Changing the second constant to six instead gives a redundant equation and x+y=3, whose solutions are (3−t,t) for real t. The coefficient matrix is identical in both cases. Singularity is a matrix property, whereas consistency also depends on b; reporting only “singular” does not classify the system’s solution set.

Counting equations alone is inadequate. A square system can be inconsistent or have infinitely many solutions if its equations are dependent. A tall system can have a unique solution when every unknown has a pivot and all extra equations are compatible. A wide system cannot constrain every unknown with coefficient pivots, but it can still be inconsistent. Always inspect the actual transformed equations and their constants instead of using shape as a complete classification.

The classification also gives a completeness argument for the parametrisation. Every assignment satisfying the echelon equations has some free-coordinate values; upward pivot equations then force exactly the reported dependent coordinates. Conversely, any allowed free-coordinate choice satisfies all pivot equations after back substitution and all remaining zero rows when no contradiction exists. Row equivalence transfers both directions to the original system, so the displayed family neither misses nor adds solutions.

Exact zero and near-zero are different in numerical computation. The rational lab can decide a coefficient is exactly zero; floating calculations need a tolerance and scale-aware interpretation. A tiny residual or a tiny computed determinant is not, by itself, an exact symbolic classification. Module 11 introduces numerical rank and Module 28 treats floating-point reliability. Here the exact small examples establish the algebraic categories before those numerical questions are added.

A square n×n matrix A is invertible when there is an n×n matrix A⁻¹ satisfying both A⁻¹A=Iₙ and AA⁻¹=Iₙ. The inverse undoes the map and gives x=A⁻¹b for each target b. If B and C are both inverses, B=B(AC)=(BA)C=C, so the inverse is unique. Square shape is a necessary part of this definition; rectangular maps can have one-sided inverses under other conditions, but those are different claims.

For a square matrix, elimination has a pivot in every column exactly when it is invertible. Full pivots make every target system consistent and uniquely solvable. Solving Aeⱼ for each coordinate unit target builds the inverse’s columns. Conversely, an invertible map cannot have a nonzero vector z with Az=0, because applying A⁻¹ would force z=0. A missing coefficient pivot gives a nonzero homogeneous solution by choosing a free variable, so it contradicts invertibility. This connects the map, equation, and pivot interpretations.

For the constructed matrix B, the column equations give AB=I. Also BAx and x both solve Ay=Ax; uniqueness for that target gives BAx=x for every x, hence BA=I. This verifies both inverse obligations rather than quietly treating a one-sided identity as the full definition.

For a two-by-two matrix with rows (a,b),(c,d), the determinant is ad−bc. Its columns form the two edges of an image parallelogram, and this expression is its signed area relative to the unit square. Absolute determinant measures area scaling; a negative sign indicates reversed orientation. Zero means the parallelogram collapses to a line or point. The determinant is not an individual output length or the number of stored entries.

For general two-by-two A, direct multiplication verifies the inverse formula $(ad-bc)^{-1}\begin{bmatrix}d&-b\\-c&a\end{bmatrix}$ when the determinant is nonzero. If ad−bc=0, the columns are dependent and the map collapses a nonzero direction, so it has no two-sided inverse. This small proof supports the area interpretation and inverse criterion without relying on a successful numerical library call.

In higher dimensions the determinant generalises signed volume scaling, and a square real matrix is invertible exactly when its determinant is nonzero. We use this general theorem with its scope stated; a full construction of the determinant from alternating multilinear volume is a follow-on topic. Under elimination, row swaps change its sign, nonzero row scaling multiplies it by that scale, and row addition preserves it. A triangular matrix’s determinant is the product of its diagonal entries, tying the criterion to its pivots.

Algebraically x=A⁻¹b is correct when an inverse exists. Numerically, use a direct system solver for the requested right-hand side rather than explicitly constructing all inverse entries. A solver can use a factorisation and substitution without forming a complete inverse, and reusable factorisations help with several targets. The lab’s well-behaved example makes both approaches agree; that does not prove inverse code always fails or that every direct solve is insensitive to data error.

Cramer’s rule expresses each unknown as a ratio of determinants, replacing one coefficient column with b in the numerator. It is an exact formula for nonsingular square systems, but evaluating many determinants separately repeats work and does not cure sensitivity in a ratio. It is useful for a small symbolic derivation, while a direct factorisation-based solve is the numerical workflow used here. An algebraically valid expression need not be the preferred computational procedure.

A residual is r=A x̂−b for a proposed solution x̂. A zero residual establishes that the proposed vector solves the stored equations in exact arithmetic. A small numerical residual measures near-satisfaction of those equations; it does not establish that the original data were accurate or that the inferred vector is close to a solution for slightly different data. Its scale and the system’s sensitivity matter.

For A=diag(1,10⁻¹²), the target b=(1,10⁻¹²) gives x=(1,1). Changing only the second target by 10⁻⁹ gives x′=(1,1001). The matrix remains invertible and x′ has zero residual for the changed target, yet the solution changed substantially. This is inherent sensitivity of the equations, not automatically a solver bug. Spectral and condition-number analysis later quantify it; this example already refutes the claim that invertibility or residual alone guarantees reliable coefficients.

For an input column x with d coordinates and output column y with o coordinates, write y=Mx+b with M of shape o×d and b of shape o×1. Every row of M selects one output’s weighted feature sum. For a batch of B input records stored as rows, use X of shape B×d and W=Mᵀ of shape d×o. Then Y=XW+bᵀ has shape B×o, with the same affine transformation applied to each row.

The transpose relationship reconciles conventions; it does not make multiplication commutative. Row-oriented XW and column-oriented Mx use the same coefficients under explicitly different storage orientations. A library’s weight layout may choose either convention. Read its definition and annotate the axes before copying a formula. The numbers B,d,o stand for observation, feature, and output counts, not interchangeable generic sizes.

In NumPy, X @ W is this matrix product and X * W is entrywise multiplication or broadcasting. When shapes happen to match, the wrong operator can run and give plausible numbers. Likewise a flat bias (o,) broadcasts over B rows, while a column bias (o,1) can fail or add different row offsets when B=o. Equality of two axis lengths can conceal a semantic mistake, so shape compatibility alone is not sufficient.

For the first two rows of the worked example, a wrong bias column adds 0.5 to both outputs of the first row and −0.5 to both outputs of the second. It gives rows (1.5,2.5),(−1.5,2.5), different from the intended output-specific bias. Both results have shape (2,2). Label the output and batch axes, and check values on a deliberately nonuniform small example as well as asserting the shape.

The elementary product algorithm uses mnp scalar multiplications for an m×n matrix times an n×p matrix. Summing each output’s n products requires n−1 additions when the first product initializes the sum, giving mp(n−1) additions; an implementation starting from zero can count n additions instead. State the charged operations before comparing counts. These formulas come directly from the finite index domains, and optional asymptotic interpretation is available in Module 06.

Affine bias addition contributes one further addition per output entry. A different parenthesisation of a longer product can change intermediate dimensions and work even though associativity preserves the exact map. Numerical libraries exploit structure and hardware, so the simple scalar-operation count is a transparent educational model, not a direct prediction of seconds on a particular machine. Keep algebraic validity, operation counts, and measured time distinct.

An implementation contract should record input/output shapes, axis meanings, numeric representation, and the intended exceptional cases. np.linalg.solve(A,b) in our tested environment expects a square nonsingular coefficient matrix for a unique exact-system solve. It does not return the family of an underdetermined system or classify a contradictory singular system. Least-squares approximation and minimum-norm solutions require their own definitions, which arrive in Module 12.

Verification has several layers. Check dimensions before multiplication, validate a small hand example, and verify residuals or inverse identities as appropriate. A known exact solution can reveal a transposed convention or misplaced bias. Numerical tolerance is chosen for the data and computation; it does not transform finite checks into a general proof. The row-operation and composition proofs explain why the formulas work beyond the particular arrays tested.

Use the NumPy preparation and its pinned numerical environment. The elimination trace additionally uses standard-library Fractions for exact arithmetic. Downloaded scripts generated the captured outputs, including the actual residuals.

Lab A Matrix products and an affine batch

Predict: SH,HS and their images of (1,1); annotate X,W,bias shapes. Run: compare the loop product with NumPy and inspect the affine output. Explain: why the rightmost column map acts first and why the final bias is output-specific. Change: replace a weight matrix with a rectangular one and update the bias/output contract accordingly. The script rejects ragged matrices and incompatible inner dimensions.

Download lab1_products.py

"""Shape-checked matrix multiplication and an explicitly batched affine map."""
import numpy as np

def shape(matrix):
    if not matrix or not matrix[0] or any(len(row) != len(matrix[0]) for row in matrix):
        raise ValueError("Require a nonempty rectangular matrix")
    return len(matrix), len(matrix[0])

def multiply(a, b):
    rows, inner = shape(a)
    other_inner, cols = shape(b)
    if inner != other_inner:
        raise ValueError("Inner dimensions must agree")
    return [[sum(a[i][k] * b[k][j] for k in range(inner)) for j in range(cols)] for i in range(rows)]

scale, shear = [[2, 0], [0, 1]], [[1, 1], [0, 1]]
for name, a, b in [("S H (shear first)", scale, shear), ("H S (scale first)", shear, scale)]:
    product = multiply(a, b)
    assert np.array_equal(np.array(product), np.array(a) @ np.array(b))
    print(name + ":", product, "; image of (1,1):", multiply(product, [[1], [1]]))
x, weights, bias = [[1, 0], [0, 1], [2, -1]], [[1, 2], [-1, 3]], [0.5, -0.5]
linear = multiply(x, weights)
affine = [[value + bias[j] for j, value in enumerate(row)] for row in linear]
reference = np.array(x) @ np.array(weights) + np.array(bias)
assert reference.shape == (3, 2) and np.array_equal(reference, np.array(affine))
print("Shapes X / W / bias / output:", shape(x), shape(weights), (len(bias),), shape(affine))
print("XW:", linear)
print("XW + bias:", affine)
for a, b in [([[1, 2]], [[1, 2]]), ([[1, 2], [3]], [[1], [2]])]:
    try:
        multiply(a, b)
    except ValueError as error:
        print("Rejected:", error)

S H (shear first): [[2, 2], [0, 1]] ; image of (1,1): [[4], [1]]
H S (scale first): [[2, 1], [0, 1]] ; image of (1,1): [[3], [1]]
Shapes X / W / bias / output: (3, 2) (2, 2) (2,) (3, 2)
XW: [[1, 2], [-1, 3], [3, 1]]
XW + bias: [[1.5, 1.5], [-0.5, 2.5], [3.5, 0.5]]
Rejected: Inner dimensions must agree
Rejected: Require a nonempty rectangular matrix

Lab B Exact pivoted elimination

Predict: which row supplies the first pivot and whether each system is unique, inconsistent, or underdetermined. Run: inspect the exact augmented rows, pivots and free columns. Explain: a zero constant row versus a contradiction and verify the unique solution in the original equations. Change: alter the duplicate-row right-hand side and predict the new classification before running. Printed column indices are zero-based; mathematical x,y,z names follow their original order.

Download lab2_elimination.py

"""Exact rational elimination with pivot selection and all three system classes."""
from fractions import Fraction

def echelon(coefficients, rhs):
    rows, cols = len(coefficients), len(coefficients[0])
    if len(rhs) != rows or any(len(row) != cols for row in coefficients):
        raise ValueError("Require rectangular coefficients and one RHS per row")
    augmented = [[Fraction(x) for x in row] + [Fraction(b)] for row, b in zip(coefficients, rhs)]
    pivot_row, pivots, events = 0, [], []
    for column in range(cols):
        if pivot_row == rows:
            break
        chosen = max(range(pivot_row, rows), key=lambda i: abs(augmented[i][column]))
        if augmented[chosen][column] == 0:
            continue  # No pivot in this column; it can be free.
        if chosen != pivot_row:
            augmented[chosen], augmented[pivot_row] = augmented[pivot_row], augmented[chosen]
            events.append(f"swap rows {pivot_row} and {chosen}")
        pivot = augmented[pivot_row][column]
        augmented[pivot_row] = [value / pivot for value in augmented[pivot_row]]
        for i in range(pivot_row + 1, rows):
            factor = augmented[i][column]
            augmented[i] = [a - factor * b for a, b in zip(augmented[i], augmented[pivot_row])]
        pivots.append(column)
        events.append(f"pivot row {pivot_row}, column {column}")
        pivot_row += 1
    inconsistent = any(all(value == 0 for value in row[:cols]) and row[cols] != 0 for row in augmented)
    status = "inconsistent" if inconsistent else "unique" if len(pivots) == cols else "underdetermined"
    return augmented, pivots, events, status

def solve_unique(rows, pivots):
    cols = len(rows[0]) - 1
    answer = [Fraction(0)] * cols
    for row_index in range(len(pivots) - 1, -1, -1):
        column = pivots[row_index]
        row = rows[row_index]
        answer[column] = row[-1] - sum(row[j] * answer[j] for j in range(column + 1, cols))
    return answer

cases = [
    ("pivoted unique", [[0, 2, 1], [1, -2, -3], [2, 3, 1]], [3, 0, 7]),
    ("inconsistent", [[1, 1], [2, 2]], [3, 7]),
    ("underdetermined", [[1, 1], [2, 2]], [3, 6]),
]
for name, coefficients, rhs in cases:
    rows, pivots, events, status = echelon(coefficients, rhs)
    print(name + ":", status)
    print("Events:", events)
    print("Echelon augmented rows:", [[str(value) for value in row] for row in rows])
    print("Pivot / free columns (zero based):", pivots, [j for j in range(len(coefficients[0])) if j not in pivots])
    if status == "unique":
        answer = solve_unique(rows, pivots)
        assert all(sum(Fraction(a) * x for a, x in zip(row, answer)) == b for row, b in zip(coefficients, rhs))
        assert answer == [1, 2, -1]
        print("Solution:", [str(x) for x in answer])
    assert status == ("unique" if name == "pivoted unique" else name)
print("All row operations used exact Fractions; a zero RHS row differs from a contradiction row.")

pivoted unique: unique
Events: ['swap rows 0 and 2', 'pivot row 0, column 0', 'pivot row 1, column 1', 'pivot row 2, column 2']
Echelon augmented rows: [['1', '3/2', '1/2', '7/2'], ['0', '1', '1', '1'], ['0', '0', '1', '-1']]
Pivot / free columns (zero based): [0, 1, 2] []
Solution: ['1', '2', '-1']
inconsistent: inconsistent
Events: ['swap rows 0 and 1', 'pivot row 0, column 0']
Echelon augmented rows: [['1', '1', '7/2'], ['0', '0', '-1/2']]
Pivot / free columns (zero based): [0] [1]
underdetermined: underdetermined
Events: ['swap rows 0 and 1', 'pivot row 0, column 0']
Echelon augmented rows: [['1', '1', '3'], ['0', '0', '0']]
Pivot / free columns (zero based): [0] [1]
All row operations used exact Fractions; a zero RHS row differs from a contradiction row.

Lab C Solvers singularity and sensitivity

Predict: solve [[3,1],[1,2]] at target [9,8], then inspect the two singular targets. Run: compare the direct and inverse results; both work for this example. Explain: why the singular exception cannot distinguish an infinite family from no solution, and why the diagonal system’s perturbed solution still has zero residual. Change: increase the small diagonal entry and compare the response to the same target perturbation. The example demonstrates data sensitivity rather than asserting that one library routine is always more accurate.

Download lab3_solvers_and_faults.py

"""Direct solves, singular RHS distinctions, and sensitivity despite a nonzero determinant."""
import numpy as np

a = np.array([[3.0, 1.0], [1.0, 2.0]])
b = np.array([9.0, 8.0])
direct = np.linalg.solve(a, b)
via_inverse = np.linalg.inv(a) @ b
assert direct.shape == (2,) and via_inverse.shape == (2,)
assert np.allclose(direct, [2, 3], rtol=1e-12, atol=1e-12)
assert np.allclose(via_inverse, direct, rtol=1e-12, atol=1e-12)
print("Direct solve:", direct.round(6).tolist())
print("Inverse times RHS:", via_inverse.round(6).tolist())
print("Residual norms:", float(np.linalg.norm(a @ direct - b)), float(np.linalg.norm(a @ via_inverse - b)))
print("Both agree here; forming an inverse is unnecessary for solving one RHS.")

singular = np.array([[1.0, 2.0], [2.0, 4.0]])
for rhs in [np.array([3.0, 6.0]), np.array([3.0, 7.0])]:
    try:
        np.linalg.solve(singular, rhs)
    except np.linalg.LinAlgError:
        kind = "infinitely many: x+2y=3" if rhs[1] == 2 * rhs[0] else "inconsistent: second RHS is not twice first"
        print("Singular solve rejected:", rhs.tolist(), ";", kind)

sensitive = np.diag([1.0, 1e-12])
baseline_rhs = np.array([1.0, 1e-12])
changed_rhs = np.array([1.0, 1e-12 + 1e-9])
baseline, changed = np.linalg.solve(sensitive, baseline_rhs), np.linalg.solve(sensitive, changed_rhs)
assert np.allclose(baseline, [1, 1], rtol=1e-12, atol=1e-12)
assert np.allclose(changed, [1, 1001], rtol=1e-12, atol=1e-12)
print("Nonzero determinant:", float(np.linalg.det(sensitive)))
print("Baseline / perturbed solution:", baseline.round(6).tolist(), changed.round(6).tolist())
print("Small RHS change:", float(np.linalg.norm(changed_rhs - baseline_rhs)))
print("Changed-system residual:", float(np.linalg.norm(sensitive @ changed - changed_rhs)))
print("Invertibility and small residual do not guarantee insensitivity of the solution.")

Direct solve: [2.0, 3.0]
Inverse times RHS: [2.0, 3.0]
Residual norms: 0.0 0.0
Both agree here; forming an inverse is unnecessary for solving one RHS.
Singular solve rejected: [3.0, 6.0] ; infinitely many: x+2y=3
Singular solve rejected: [3.0, 7.0] ; inconsistent: second RHS is not twice first
Nonzero determinant: 1.000000000000001e-12
Baseline / perturbed solution: [1.0, 1.0] [1.0, 1001.0]
Small RHS change: 1e-09
Changed-system residual: 0.0
Invertibility and small residual do not guarantee insensitivity of the solution.

Exercises 1–12 are required; 13–14 are extensions. Use exact arithmetic for the small systems and state the solution set rather than only a solver status.

Automatic scoring covers Questions 1–9; Question 10 is written and self-reviewed.

Use linear-system, matrix-product, and determinant selections in the authors’ Mathematics for Machine Learning companion. Read the official NumPy solve contract and compare the inverse interface. These primary readings supplement the original examples and arguments here.

The library’s phrase exact-system solve refers to solving Ax=b rather than a least-squares objective; floating-point implementation is still numerical. Do not read it as a guarantee of symbolic exactness or data insensitivity.

Without notes, explain columns as coordinate-direction images, derive a product entry, and prove row-addition equivalence. Solve and classify one system in each category. State the determinant and residual claims separately, then annotate a batched affine calculation.

Exit task: take A=[[1,2],[0,1]],b=(5,2). Back substitution gives y=2,x=1; detA=1. Its inverse [[1,−2],[0,1]] recovers the input. Explain how b stored as (2,) versus (2,1) changes storage shape without changing this one mathematical target. Neither shape alone creates a batch of two independent target systems.

Ready to move on: you can interpret equations, classify their complete solution sets, and keep map and storage conventions consistent. The next lesson develops spaces, bases, rank, null directions, and parameter identifiability. See the course overview for availability.