Broadcasting

In the previous lesson all the operations were done between an array and a number, and the shape never changed. In real calculations two arrays meet, and their shapes are most often not the same: next to the reading table of 3,780 cells stand a base per period (3 values), a household count per subscriber (1,260 values) and a unit price per tariff tier (3 values). Instead of matching these one by one with loops, the array interface applies the broadcasting rule: the smaller shape is counted as repeated to fit the larger shape.

The rule is told in two lines, and it is dangerous precisely because it is only two lines. Shapes are aligned from the right; in every aligned dimension either the two numbers are equal or one of them is

1. If this condition holds, the operation is carried out. The condition holding does not mean the operation is correct, only that it is possible. What this lesson measures is the difference between the two.

• NC7. The broadcasting rule checks the compatibility of shapes, not the meaning of the operation. When the rule holds, the calculation is done and no warning is produced.
• NC8. An incompatible shape pair gives one of two results: an error or an array of unexpected size. The first stops the calculation, the second continues it and can enlarge the element count by thousands of times.

Writing the Rule

# broadcasting.py — MODEL: the broadcasting rule is written with the standard library. The rule:
# shapes are aligned from the right; in every dimension either the two numbers are equal or one is 1.
# If the rule holds, the operation is carried out. No library is called.
import math

def broadcast_shape(shape1, shape2):
    n = max(len(shape1), len(shape2))
    s1 = (1,) * (n - len(shape1)) + tuple(shape1)
    s2 = (1,) * (n - len(shape2)) + tuple(shape2)
    c = []
    for x, y2 in zip(s1, s2):
        if x == y2 or x == 1 or y2 == 1:
            c.append(max(x, y2))
        else:
            return None
    return tuple(c)

count = lambda b: 0 if b is None else math.prod(b)
PAIRS = [
    ("reading x single coefficient",         (1260, 3), (1,),        3780),
    ("reading x period base",                (1260, 3), (3,),        3780),
    ("reading x subscriber household count", (1260, 3), (1260, 1),   3780),
    ("reading x subscriber (flat form)",     (1260, 3), (1260,),     3780),
    ("reading x period (two dimensions)",   (1260, 3), (1, 3),      3780),
    ("subscriber total x household row",     (1260, 1), (1, 1260),   1260),
    ("period base x period adjustment",      (3,),      (3, 1),      3),
    ("reading x two tiers",                  (1260, 3), (2,),        3780),
]
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([38, 12, 11, 13, 11, 11, 9])
print(y(["operation", "shape 1", "shape 2", "result shape", "elements", "expected", "outcome"]))
errors = silent = correct = 0
for name, shape1, shape2, expected in PAIRS:
    c = broadcast_shape(shape1, shape2)
    if c is None:
        status, errors = "error", errors + 1
    elif count(c) != expected:
        status, silent = "silent", silent + 1
    else:
        status, correct = "correct", correct + 1
    print(y([name, str(shape1), str(shape2), str(c) if c else "-", count(c), expected, status]))
print(f"\n{len(PAIRS)} pairs: {correct} correct, {errors} errors, {silent} silently unexpected shape")

operation                                  shape 1    shape 2 result shape   elements   expected  outcome
reading x single coefficient             (1260, 3)       (1,)    (1260, 3)       3780       3780  correct
reading x period base                    (1260, 3)       (3,)    (1260, 3)       3780       3780  correct
reading x subscriber household count     (1260, 3)  (1260, 1)    (1260, 3)       3780       3780  correct
reading x subscriber (flat form)         (1260, 3)    (1260,)            -          0       3780    error
reading x period (two dimensions)        (1260, 3)     (1, 3)    (1260, 3)       3780       3780  correct
subscriber total x household row         (1260, 1)  (1, 1260) (1260, 1260)    1587600       1260   silent
period base x period adjustment               (3,)     (3, 1)       (3, 3)          9          3   silent
reading x two tiers                      (1260, 3)       (2,)            -          0       3780    error

8 pairs: 4 correct, 2 errors, 2 silently unexpected shape

The first three rows are where the rule works, and the gain is real. Subtracting a base per period takes two nested loops and an index calculation when written by hand; with broadcasting it is a single operation, and the base of 3 values behaves as if it were repeated 1,260 times — without the repeated copy being built in memory. The household count per subscriber is spread in the right direction in the same way, as long as the shape is written (1260, 1).

The fourth row is where the rule protects. If the household count per subscriber is written as a flat array, that is, in the shape (1260,), the right alignment compares the last dimensions: 3 and 1,260. They are not equal to each other and neither is 1, so the rule does not hold and the operation is not carried out. This is not a flaw but the one real guarantee of the rule: an attempt to align on the wrong axis is, more often than not, stopped. The eighth row is of the same kind; a coefficient array of two tiers cannot be broadcast onto a table of three periods.

The sixth and seventh rows, however, are not stopped. The rule holds, the operation is carried out, and the result has a shape other than the expected one. Two of the eight pairs are in this situation, and what the two have in common is this: the person writing has written one of the shapes transposed. The next section counts the cost of this.

The Silently Produced Result

# silent.py — MODEL: the cost of the silently produced result. The total consumption per
# subscriber is divided by the household count per subscriber. The household count is fictional.
import sys
from array import array

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)]
r3 = rng(31081977)
HOUSEHOLD = [1 + int(r3() * 5) for _ in range(SUBSCRIBERS)]
TOTALS = [sum(READINGS[i*PERIODS:(i+1)*PERIODS]) for i in range(SUBSCRIBERS)]

# Right form: (1260, 1) with (1260, 1)
RIGHT = array("d", [TOTALS[i] / HOUSEHOLD[i] for i in range(SUBSCRIBERS)])
# Silent form: (1260, 1) with (1, 1260) -> (1260, 1260)
WRONG = array("d", [TOTALS[i] / HOUSEHOLD[j] for i in range(SUBSCRIBERS) for j in range(SUBSCRIBERS)])
diagonal = [WRONG[i * SUBSCRIBERS + i] for i in range(SUBSCRIBERS)]
correct_count = sum(1 for i in range(SUBSCRIBERS) if abs(diagonal[i] - RIGHT[i]) < 1e-9)
z = row_fmt([26, 18, 20, 12])
print(z(["measure", "right form", "silent form", "ratio"]))
print(z(["result shape", "(1260, 1)", "(1260, 1260)", "-"]))
print(z(["element count", len(RIGHT), len(WRONG), f"{len(WRONG)//len(RIGHT)}x"]))
print(z(["bytes allocated", sys.getsizeof(RIGHT), sys.getsizeof(WRONG),
         f"{sys.getsizeof(WRONG)//sys.getsizeof(RIGHT)}x"]))
print(z(["per-person mean", f"{sum(RIGHT)/len(RIGHT):.6f}", f"{sum(WRONG)/len(WRONG):.6f}",
         "-"]))
print(f"\nof the {len(WRONG)} elements of the silent array, {correct_count} are correct (diagonal), "
      f"{len(WRONG)-correct_count} elements divide one subscriber's consumption by another subscriber's household count")
print(f"reported per-person mean {sum(WRONG)/len(WRONG):.6f} instead of "
      f"{sum(RIGHT)/len(RIGHT):.6f} m3, difference {abs(sum(WRONG)/len(WRONG)-sum(RIGHT)/len(RIGHT)):.6f} m3")

# Second silent case: period base (3,) with seasonal adjustment (3, 1).
PERIOD_MEAN = [sum(READINGS[i*PERIODS+d] for i in range(SUBSCRIBERS)) / SUBSCRIBERS for d in range(PERIODS)]
ADJUSTMENT = [0.98, 1.03, 1.10]   # fictional seasonal adjustment
adjusted = [PERIOD_MEAN[d] * ADJUSTMENT[d] for d in range(PERIODS)]
adjusted_wrong = [PERIOD_MEAN[d] * ADJUSTMENT[e] for e in range(PERIODS) for d in range(PERIODS)]
print(f"\nsecond silent case: in the product of period base (3,) and adjustment (3, 1) the result is "
      f"{len(adjusted_wrong)} values instead of {len(adjusted)}")
print(f"adjusted base total {sum(adjusted_wrong):.6f} m3 instead of {sum(adjusted):.6f} m3; "
      f"because the shape grew, the total rose to {sum(adjusted_wrong)/sum(adjusted):.4f} times")

measure                           right form         silent form       ratio
result shape                       (1260, 1)        (1260, 1260)           -
element count                           1260             1587600       1260x
bytes allocated                        10160            12700880       1250x
per-person mean                    36.938180           36.497564           -

of the 1587600 elements of the silent array, 1260 are correct (diagonal), 1586340 elements divide one subscriber's consumption by another subscriber's household count
reported per-person mean 36.497564 instead of 36.938180 m3, difference 0.440616 m3

second silent case: in the product of period base (3,) and adjustment (3, 1) the result is 9 values instead of 3
adjusted base total 251.846368 m3 instead of 83.901904 m3; because the shape grew, the total rose to 3.0017 times

The calculation is this: each subscriber’s three-period total consumption is to be divided by that subscriber’s household count. The correct result is 1,260 values. When the household count array is written as a row, that is, in the shape (1, 1260) instead of (1260, 1), the rule still holds — in the first dimension 1,260 meets 1, in the second dimension 1 meets 1,260, and in both of them one of the sides is 1. The result is an array of 1,587,600 elements, and in it every subscriber’s consumption stands divided by every subscriber’s household count.

The cost has three separate measures. The element count grows 1,260 times, the allocated bytes from 10,160 to 12,700,880, that is, 1,250 times; as long as it fits in memory this growth produces no warning, and precisely because it fits it goes unnoticed. The accuracy measure is sharper: of the 1,587,600 elements only 1,260 are correct, and these are the ones on the diagonal. The remaining 1,586,340 elements divide one subscriber’s consumption by another subscriber’s household count, that is, they have no real counterpart.

The reported number, however, reflects neither of these two measures. The per-person mean consumption comes out 36.938180 m³ in the right array and 36.497564 m³ in the silent array; the difference between them is 0.440616 m³, about 1.2%. This is the most dangerous side of broadcasting: the wrong result looks plausible. When a finding moves between 36.5 and 36.9, nobody thinks of checking the shape. The number itself does not give the error away; the only thing that gives the error away is the element count, and that number usually does not enter the report.

The second silent case is smaller and for exactly that reason more insidious. The average base of three periods is to be multiplied by three seasonal adjustments; when the adjustment is written as a column, the result is 9 values instead of 3. The adjusted base total comes out 251.846368 m³ instead of 83.901904 m³, that is, 3.0017 times. The multiplier here is not exactly 3, and this detail matters: had an exact multiple been seen, one could have said “I must have multiplied by three somewhere” and gone looking; 3.0017 looks like a number that comes from the calculation itself.

The working rule that comes out of these two cases is one sentence: in every operation that uses broadcasting, the shape of the result is written out and compared with the expected element count. It cannot be assumed that the result is correct because the rule holds; the rule holding says only that the calculation could be done. In two of the eight pairs above the calculation was not done, and that was the good news; in two it was done, and that was the bad news.

What the Rule Gains

# gain.py — MODEL: the gain of broadcasting. A base per period is to be subtracted.
# The by-hand way copies the base to 3,780 elements; broadcasting does not copy, it computes the index.
BASE = array("d", [9.80, 10.40, 11.20])   # fictional period base, m3

def by_copying(v, t):
    wide = array("d", [t[j % PERIODS] for j in range(len(v))])
    return array("d", [v[j] - wide[j] for j in range(len(v))]), sys.getsizeof(wide)

def by_broadcasting(v, t):
    return array("d", [v[j] - t[j % PERIODS] for j in range(len(v))]), 0

K1, extra1 = by_copying(READINGS, BASE)
K2, extra2 = by_broadcasting(READINGS, BASE)
deviating = sum(1 for a, b in zip(K1, K2) if abs(a - b) > 1e-9)
w = row_fmt([26, 16, 16, 12])
print(w(["measure", "by copying", "by broadcasting", "ratio"]))
print(w(["base array bytes", sys.getsizeof(BASE), sys.getsizeof(BASE), "1x"]))
print(w(["extra bytes allocated", extra1, extra2, "-"]))
print(w(["total bytes allocated", extra1 + sys.getsizeof(K1), extra2 + sys.getsizeof(K2),
         f"{(extra1+sys.getsizeof(K1))/(extra2+sys.getsizeof(K2)):.1f}x"]))
print(f"\nboth paths give the same result: {deviating} elements deviate above the 1e-9 m3 threshold")
print(f"broadcasting applies the 3-value base to {len(READINGS)} elements without copying; "
      f"in the copied version the base takes {extra1} bytes instead of {sys.getsizeof(BASE)}")

measure                         by copying by broadcasting       ratio
base array bytes                       104             104          1x
extra bytes allocated                30320               0           -
total bytes allocated                60640           30320        2.0x

both paths give the same result: 0 elements deviate above the 1e-9 m3 threshold
broadcasting applies the 3-value base to 3780 elements without copying; in the copied version the base takes 30320 bytes instead of 104

The gain of the rule stands in this table, and it has nothing to do with the computed result — both paths give the same 3,780 values, and there is not a single element deviating above the 1e-9 m³ threshold. The difference is in copying. The most direct way of applying the period base by hand is to expand the 3 values to 3,780 elements and subtract the two equal shapes element by element; this expansion takes 30,320 bytes, and the total allocated bytes rise to 60,640. Broadcasting never builds the expanded array; it computes the index of the small shape from the large shape, and the base stays at 104 bytes. The total allocated bytes fall by half, to 30,320.

The ratio here is 2 times because the expanded side is the same size as the result array. The smaller the ratio of the small array to the large one, the greater the gain: if the same operation were done with 3 tariff tiers instead of 3 periods and a table of 100,000 rows, the copied version would hold two full tables and broadcasting one. What broadcasting gains is not the calculation but the array that is not built, and it is worth noting that the cost in the previous section comes from the same place: because the array that is not built is cheap, the array of unexpected shape is also built cheaply and grows without drawing attention.

Summary

• The broadcasting rule aligns shapes from the right and looks in every dimension for equality or 1; the rule checks shape compatibility, not the meaning of the operation.
• Of the eight shape pairs, 4 give the correct result, 2 give an error, 2 silently produce an array of unexpected shape.
• When the household count is written as a row, the result is 1,587,600 elements instead of 1,260; the allocated bytes grow from 10,160 to 12,700,880, that is, 1,250 times.
• Of the 1,587,600 elements of the silent array only 1,260 are correct; the remaining 1,586,340 elements have no real counterpart, and the reported per-person mean comes out 36.497564 instead of 36.938180 m³.
• The period base adjustment produces 9 values instead of 3 and carries the adjusted base total from 83.901904 m³ to 251.846368 m³, to 3.0017 times.
• What the rule gains is the array that is not built: applying the period base by copying asks for 30,320 bytes of extra space, broadcasting asks for 0 bytes, and the total allocated bytes fall from 60,640 to 30,320.

Next Step

Broadcasting matches two arrays of the same shape element by element, and this is the simplest multi-array operation there is. The tariff calculation itself, however, is not an element-by-element match: when the unit prices of three tiers and the consumption of three periods are multiplied and summed, every result value is the joint calculation of a whole row and a whole column. The next lesson writes this operation — matrix multiplication and the solution problem derived from it — with a hand loop and with the abstraction, compares the pass counts and measures where numerical precision comes into play. Shape compatibility is decisive there too, but the rule changes: in broadcasting compatibility is sought by right alignment, while in multiplication a single equality is sought, and when that equality does not hold the operation cannot be done under any condition.