Basic Linear Algebra Operations

The previous lesson covered operations that match two arrays element by element: broadcasting brings together two values at the same position. The tariff calculation itself does not work like that. When a subscriber’s three-period reading is distributed over three tariff tiers, every result value is the joint calculation of a whole row and a whole column. This lesson measures these two operations — matrix multiplication and the solve problem derived from it.

The subject of the measurement is not algebra; the theory is not taught here. What is measured is how many operations two writings that produce the same result take and which decision calling the abstraction takes on. The sharpest numbers of this topic come from here, because in the solve problem the gain of abstraction is not the line count but the choice of algorithm.

• NC9. In matrix multiplication, abstraction does not change the operation count. What it changes is the number of lines at the call site and the disappearance of hand-written index computation.
• NC10. In solving, calling the abstraction is choosing an algorithm. The operation counts of the definition-based way and of elimination are not of the same order of magnitude, and the gap grows quickly with size.
• NC11. The same operation count can give two different results. The pivoting decision determines the size of the residual; input sensitivity, on the other hand, is independent of the calculation, and in an ill-conditioned system a difference the size of the measurement resolution moves the result by half.

The Same Operation, Another Writing

# matmul.py — MODEL: matrix multiplication is written first by hand with three loops, then
# as a single call. Reading table (1260, 3), tariff matrix (3, 3); the result is each
# subscriber's share in the three tariff tiers. The dataset and the tariff are fictional.
import inspect

def rng(seed):
    s = seed & 0xFFFFFFFF
    def next_value():
        nonlocal s
        s = (1664525 * s + 1013904223) & 0xFFFFFFFF
        return s / 4294967296
    return next_value

SUBSCRIBERS, PERIODS = 1260, 3
r = rng(20260218)
RAW = [round(4 + r() * 56, 2) for _ in range(SUBSCRIBERS * PERIODS)]
r2 = rng(9110723)
p = [r2() for _ in range(SUBSCRIBERS * PERIODS)]
EMPTY = set(sorted(range(len(p)), key=lambda j: p[j])[:581])
READINGS = [0.0 if j in EMPTY else RAW[j] for j in range(SUBSCRIBERS * PERIODS)]
# Tariff matrix: the distribution ratio of each period's reading over the three tiers (10/25/40 m3).
TARIFF = [0.62, 0.28, 0.10,
          0.55, 0.31, 0.14,
          0.48, 0.34, 0.18]

OPS = 0     # one multiplication and one addition
INDEX = 0   # number of hand-written index computations

def multiply_by_hand(A, na, ma, B, mb, nb):
    global OPS, INDEX
    C = [0.0] * (na * nb)
    for i in range(na):
        for j in range(nb):
            t = 0.0
            for k in range(ma):
                t += A[i * ma + k] * B[k * nb + j]
                OPS += 1
                INDEX += 2
            C[i * nb + j] = t
            INDEX += 1
    return C

def multiply(A, ba, B, bb):
    return multiply_by_hand(A, ba[0], ba[1], B, bb[0], bb[1])

OPS, INDEX = 0, 0
C1 = multiply_by_hand(READINGS, SUBSCRIBERS, PERIODS, TARIFF, PERIODS, PERIODS)
i1, d1 = OPS, INDEX
OPS, INDEX = 0, 0
C2 = multiply(READINGS, (SUBSCRIBERS, PERIODS), TARIFF, (PERIODS, PERIODS))
i2, d2 = OPS, INDEX
lines = lambda f: len(inspect.getsource(f).strip().splitlines())
row_fmt = lambda w: lambda h: "".join(str(v).ljust(w[0]) if j == 0 else str(v).rjust(w[j])
                                      for j, v in enumerate(h))
y = row_fmt([28, 20, 16, 10])
deviating = sum(1 for a, b in zip(C1, C2) if abs(a - b) > 1e-9)
print(f"(1260, 3) x (3, 3) -> (1260, 3); result {len(C1)} elements, {deviating} elements deviate "
      f"above the 1e-9 threshold\n")
print(y(["measure", "three hand loops", "single call", "ratio"]))
print(y(["lines at call site", lines(multiply_by_hand), lines(multiply), f"{lines(multiply_by_hand)/lines(multiply):.1f}x"]))
print(y(["multiply-add count", i1, i2, f"{i1/i2:.1f}x"]))
print(y(["hand index computations", d1, 0, "-"]))
print(f"\nmultiplication count {SUBSCRIBERS} x {PERIODS} x {PERIODS} = {SUBSCRIBERS*PERIODS*PERIODS}; abstraction does not "
      f"change this number, it only removes {d1} index computations from the call site")

(1260, 3) x (3, 3) -> (1260, 3); result 3780 elements, 0 elements deviate above the 1e-9 threshold

measure                         three hand loops     single call     ratio
lines at call site                            13               2      6.5x
multiply-add count                         11340           11340      1.0x
hand index computations                    26460               0         -

multiplication count 1260 x 3 x 3 = 11340; abstraction does not change this number, it only removes 26460 index computations from the call site

The number that emerges when the gain of abstraction in matrix multiplication is measured is a warning. The multiplication count is 11,340 in both writings — this number comes from the shapes, not from the writing, and no call lowers it. The only thing that changes is where it is written: the lines at the call site fall from 13 to 2, and 26,460 hand-written index computations disappear. Each of the expressions A[i * ma + k] and B[k * nb + j] is the hand-written form of reading flat memory as two dimensions, and swapping the factors in one of these two expressions produces no error, only a different number.

This is the real name of the gain: abstraction does not speed up the calculation, it lowers the number of places where a mistake can be written. When the scale grows this does not change; the multiplication count is the product of the shapes, and a library can change it only by choosing a different algorithm. The next section shows exactly where that is done.

The Algorithm the Abstraction Chooses

# solve.py — MODEL: the same system of equations is solved in two ways. The first is written
# from the definition (the determinant by cofactor expansion), the second by elimination. An operation
# counter counts both; no library is called, both algorithms are written out in the open.
COF = 0
ELIM = 0

def det_cofactor(M, n):
    global COF
    if n == 1:
        return M[0]
    t = 0.0
    for j in range(n):
        minor = [M[i*n+k] for i in range(1, n) for k in range(n) if k != j]
        COF += 1
        t += (1 if j % 2 == 0 else -1) * M[j] * det_cofactor(minor, n - 1)
    return t

def solve_elimination(M, b, n, pivot=True):
    global ELIM
    A = [list(M[i*n:(i+1)*n]) + [b[i]] for i in range(n)]
    for k in range(n):
        if pivot:
            best = max(range(k, n), key=lambda i: abs(A[i][k]))
            A[k], A[best] = A[best], A[k]
        for i in range(k + 1, n):
            c = A[i][k] / A[k][k]
            for j in range(k, n + 1):
                A[i][j] -= c * A[k][j]
                ELIM += 1
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        t = A[i][n]
        for j in range(i + 1, n):
            t -= A[i][j] * x[j]
            ELIM += 1
        x[i] = t / A[i][i]
    return x

cofactor_ops = lambda n: 0 if n <= 1 else n * (1 + cofactor_ops(n - 1))
rk = rng(70117)
def square_system(n):
    M = [round(1 + rk() * 9, 3) for _ in range(n * n)]
    b = [round(1 + rk() * 9, 3) for _ in range(n)]
    return M, b

w = row_fmt([6, 18, 14, 16, 12])
print(w(["n", "cofactor ops", "formula", "elimination ops", "ratio"]))
for n in (2, 3, 4, 5, 6, 7, 8):
    M, b = square_system(n)
    COF, ELIM = 0, 0
    det_cofactor(M, n)
    k = COF
    solve_elimination(M, b, n)
    print(w([n, k, cofactor_ops(n), ELIM, f"{k/max(ELIM,1):.0f}x"]))
for n in (12,):
    M, b = square_system(n)
    ELIM = 0
    x = solve_elimination(M, b, n)
    residual = max(abs(sum(M[i*n+j] * x[j] for j in range(n)) - b[i]) for i in range(n))
    print(w([n, "-", cofactor_ops(n), ELIM, f"{cofactor_ops(n)/ELIM:.0f}x"]))
    print(f"\nfor n = {n} the cofactor expansion asks for {cofactor_ops(n)} operations; elimination "
          f"finishes in {ELIM} operations, a ratio of {cofactor_ops(n)/ELIM:.0f} times")
    print(f"the solution of elimination satisfies the equations: largest residual {residual:.2e}, "
          f"threshold 1e-9 (coefficients of magnitude 1e1, resolution of the 64-bit float type 1e-16)")

n           cofactor ops       formula elimination ops       ratio
2                      2             2               4          0x
3                      9             9              14          1x
4                     40            40              32          1x
5                    205           205              60          3x
6                   1236          1236             100         12x
7                   8659          8659             154         56x
8                  69280         69280             224        309x
12                     -     823059744             704    1169119x

for n = 12 the cofactor expansion asks for 823059744 operations; elimination finishes in 704 operations, a ratio of 1169119 times
the solution of elimination satisfies the equations: largest residual 2.66e-15, threshold 1e-9 (coefficients of magnitude 1e1, resolution of the 64-bit float type 1e-16)

The question is this: the unit prices of the three tariff tiers are unknown, but the total bill of several subscriber groups and each group’s consumption in the tiers are known. This is a system of equations and it is solved in two ways. The way written from the definition computes the determinant by cofactor expansion; the elimination way subtracts rows from one another and isolates the unknowns one by one.

The table shows that the two ways are not of the same order of magnitude. Up to n = 3 elimination does more work; the direction turns at n = 4, and the cofactor expansion rises above the 32 operations of elimination with 40 operations. At n = 6 the ratio is 12 times, at n = 8 it is 309 times. The measured numbers match the formula exactly — the operation count of the cofactor expansion follows the recursion $T(n) = n \cdot (1 + T(n-1))$, and this is a growth in which the multiplier itself grows as n grows. For n = 12 it comes to 823,059,744 operations; elimination finishes the same system in 704 operations. The ratio is 1,169,119 times.

This is the form of the measure that carries over to the course: calling a solve function is not a shortcut, it is a choice of algorithm. The gain is not in the line count but in the way the operation count grows. For the same reason, a solve function called without knowing which algorithm was chosen is a decision made without knowing which operation count was agreed to.

Numerical Sensitivity

# sensitivity.py — MODEL: the unit price of three tariff tiers is solved from the total bill of
# three subscriber groups. The coefficients are each group's m3 amount in the tiers;
# they are fictional. The consumption profiles of the first two groups are very close to each other.
matvec = lambda M, x, n: [sum(M[i*n+j] * x[j] for j in range(n)) for i in range(n)]
TRUE_PRICES = [12.50, 18.75, 27.40]     # unit price, fictional
COEF = [10.00, 15.00, 5.00,
        10.02, 15.04, 5.01,
        8.00, 20.00, 12.00]
BILL = matvec(COEF, TRUE_PRICES, 3)
x0 = solve_elimination(COEF, BILL, 3)
COEF2 = list(COEF)
COEF2[0] = COEF[0] + 0.01               # a single change the size of the meter resolution
x1 = solve_elimination(COEF2, BILL, 3)
u = row_fmt([16, 16, 18, 18])
print(f"determinant of the coefficient matrix {det_cofactor(COEF, 3):.4f}\n")
print(u(["tier", "solved price", "after 0.01 m3", "deviation percent"]))
for i, label in enumerate(["10 m3", "25 m3", "40 m3"]):
    print(u([label, f"{x0[i]:.6f}", f"{x1[i]:.6f}",
             f"{100*abs(x1[i]-TRUE_PRICES[i])/TRUE_PRICES[i]:.2f}"]))
print(f"\na single 0.01 m3 change in the input moves the price of the first tier by "
      f"{100*abs(x1[0]-TRUE_PRICES[0])/TRUE_PRICES[0]:.2f} percent; there is no error in the calculation")

# The pivoting decision. As the first coefficient shrinks, elimination without pivoting breaks down;
# the value here is not a real reading, it is a probe that measures the behavior of the algorithm.
v = row_fmt([18, 20, 20, 16])
print("\n" + v(["first coefficient", "pivoted residual", "pivotless residual", "operation count"]))
for eps in (1e-2, 1e-6, 1e-10, 1e-14, 0.0):
    T = [eps, 15.00, 5.00, 10.02, 15.04, 5.01, 8.00, 20.00, 12.00]
    bt = matvec(T, TRUE_PRICES, 3)
    residual_of = lambda xx: max(abs(sum(T[i*3+j] * xx[j] for j in range(3)) - bt[i])
                                 for i in range(3))
    ELIM = 0
    xp = solve_elimination(T, bt, 3, True)
    n_ops = ELIM
    try:
        xn = solve_elimination(T, bt, 3, False)
        pr = f"{residual_of(xn):.3e}"
    except ZeroDivisionError:
        pr = "no solution"
    print(v([f"{eps:g}", f"{residual_of(xp):.3e}", pr, n_ops]))
print("\npivoting does not change the operation count; what it changes is the size of the residual")

determinant of the coefficient matrix 0.8000

tier                solved price     after 0.01 m3 deviation percent
10 m3                  12.500000          6.239082             50.09
25 m3                  18.750000         25.001560             33.34
40 m3                  27.400000         21.154679             22.79

a single 0.01 m3 change in the input moves the price of the first tier by 50.09 percent; there is no error in the calculation

first coefficient     pivoted residual  pivotless residual operation count
0.01                         0.000e+00           8.538e-11              14
1e-06                        0.000e+00           5.687e-07              14
1e-10                        0.000e+00           4.120e-03              14
1e-14                        0.000e+00           7.237e+01              14
0                            0.000e+00         no solution              14

pivoting does not change the operation count; what it changes is the size of the residual

The upper table measures a sensitivity outside the calculation. The consumption profiles of the first two subscriber groups are very close to each other; the determinant of the coefficient matrix is a small number like 0.8. In this system, when the size of the meter resolution, that is, 0.01 m³, is added to the input, the unit price of the first tier falls from 12.50 to 6.24, a deviation of 50.09%. The other two tiers also move by 33% and 23%. There is no error in the calculation; both results are the correct solution of their own input. What moves is the input itself, and since the input came with the measurement, it carries the uncertainty of 0.01 m³.

This is the most important limit of the topic: if the result changes by half while the difference between two measurements is too small to be read from the meter, that result is not a reportable number. The value a solve function gives can be correct and meaningless at the same time; what separates the two is moving the input and looking at how far the result moves.

The lower table shows the decision that abstraction hides. Pivoting is bringing the largest coefficient to the front at every step of elimination, and it does not change the operation count at all: 14 in all five rows. What it changes is the residual. When the first coefficient is 0.01, the residual of elimination without pivoting, 8.538e-11, is below the threshold and the two ways cannot be told apart. As the coefficient shrinks the residual grows: 4.120e-03 at 1e-10, 7.237e+01 at 1e-14 — that is, the result no longer satisfies the equations at all. When the coefficient is zero, elimination without pivoting cannot even produce a solution, while elimination with pivoting gives the correct result in the same 14 operations. If it is not written whether pivoting is on when calling a solve function, nothing can be claimed about the size of the residual.

Summary

• In matrix multiplication, abstraction does not change the operation count: 11,340 multiplications are the same in both writings; the lines at the call site fall from 13 to 2 and 26,460 hand-written index computations disappear.
• In solving, the gain is a choice of algorithm: the cofactor expansion and elimination are on par up to n = 4, and they separate by 309 times at n = 8 and by 1,169,119 times at n = 12, 823,059,744 operations against 704.
• The solution of elimination satisfies the equations: at n = 12 the largest residual is 2.66e-15, far below the 1e-9 threshold.
• Pivoting does not change the operation count (fixed at 14), it changes the residual: as the first coefficient shrinks, the residual without pivoting rises from 8.538e-11 to 7.237e+01, and when the coefficient is zero no solution is produced at all.
• In an ill-conditioned tariff system (determinant 0.8), when 0.01 m³, the size of the meter resolution, is added to the input, the unit prices deviate by 50.09%, 33.34% and 22.79%; there is no error in the calculation.

Next Step

In five lessons the structure of numerical computing was built, and at every step two numbers stood side by side. When the environment was pinned, the notebook gave the same result in 17 of 17 steps, the array interface did the same job in 2 lines instead of 22 and in 30,320 bytes instead of 188,864, vectorization cut the lines to a fifth while tripling the passes, broadcasting saved 30,320 bytes with the array that is not built while silently computing 1,586,340 elements wrong with an incompatible shape, and the solve function put a 1,169,119-fold difference in operation count behind a single call.

The numbers got faster, but there is something repeated in all the code blocks of these five lessons: the only thing that says which column is what is the code. The expression READINGS[i * PERIODS + d] carries no meaning without a person who knows the second dimension is the period; that the rows of the TARIFF matrix are periods and its columns are tiers is written nowhere; which of the 581 empty cells are which cannot be read from the array. The array has no label. The next topic, dataframes, takes up exactly this gap: once labels are added, which jobs become correct on their own, and which decision does the alignment that the label takes on make invisible.