Activation Functions

The previous lesson measured that nonlinear separation comes from the activation function, but the two functions it used were not chosen with any justification. What an activation adds to a network is not only its curve; that function also determines the gradient that flows backward during training. The gradient is multiplied layer by layer through the chain rule to reach the earlier layers, and one of the factors is the derivative of the activation.

The measure of this lesson lies here. If the derivative of a function falls to nearly zero in one region of the input, a unit that stays in that region receives no gradient and its learning stops. This has two separate forms, and they are confused: the derivative shrinks but does not become zero — saturation; or the derivative is exactly zero — a dead unit. The two are measured separately.

• NN8. The set, the split, the seed and the twelve columns are as in 01; the training set is again reduced to 300 rows. The baseline on the test set is 0.7579.
• NN9. The network is the network of 01: 12 inputs, a hidden layer of eight units, a single output, 113 parameters. The only thing that changes is the activation of the hidden layer.
• NN10. The output layer is a sigmoid in all three setups, because the output has to be a probability. The comparison is made only through the hidden layer.
• NN11. The derivative is written in terms of the output of the activation: $a(1-a)$ for the sigmoid, $1-a^2$ for the hyperbolic tangent, and for the rectified unit 1 when the output is positive and 0 otherwise.
• NN12. A unit is not counted as active on a row if its derivative in the hidden layer stays below 5% of its ceiling. A unit that is not active on any of the 300 rows of the training set is dead.
• NN13. The gradient reaching the first layer is the mean absolute magnitude of the weight gradients of the first epoch. Since the three setups start from the same initial values, the numbers are comparable.
• NN14. The learning rate sweep is not a tuning search; it is made to show at which rate a dead unit appears, and it makes no choice.

Saturation Is a Number

When saturation is described in words, it remains a description of a curve. Once it is turned into a number, how sharp a boundary it is becomes visible: by how many times the derivative of the same function differs at two points of the input.

# activations.py — MODEL: the same FICTIONAL subscriber table, the same seed and the same split;
# the training set is reduced to 300 rows as in 01.
import math
import statistics

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
REGION = [("north", 0.28, 21, 0.00), ("south", 0.22, 17, -0.10),
          ("east", 0.18, 26, 0.30), ("west", 0.14, 14, -0.05),
          ("central", 0.18, 23, 0.05)]


def rng(t):
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

    def next_value():
        s[0] = (s[0] * 1664525 + 1013904223) & M32
        return s[0] / 4294967296
    return next_value


def pick(u, w):
    t = 0.0
    for i, x in enumerate(w):
        t += x
        if u < t:
            return i
    return len(w) - 1


SUBSCRIBERS, DATA = [], []
for i in range(RAW):
    r = rng(SEED + i)
    b = REGION[pick(r(), [x[1] for x in REGION])]
    household = 1 + pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05])
    rating = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17])
    if r() >= 0.046:
        SUBSCRIBERS.append({"id": 10001 + i, "region": b, "household": household,
                            "satisfaction": rating})
for k in SUBSCRIBERS:
    r, v = rng(SEED + 7000 + k["id"]), []
    for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        v.append(0.0 if r() < 0.038 else math.floor(
            k["region"][2] * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - n * 0.05) * 100 + 0.5) / 100)
        r()
    if not v:
        continue
    mean, vol, r = sum(v) / len(v), max(v) - min(v), rng(SEED + 51000 + k["id"])
    z = (0.052 * (mean - 20) + 0.85 * (min(v) == 0.0) + 0.026 * vol
         + 0.24 * (k["household"] >= 5) - 0.20 * (k["satisfaction"] >= 4) + k["region"][3]
         + 0.9 * (k["id"] - 10001) / 1399 + (r() + r() + r() - 1.5) * 1.30)
    x = {"household": k["household"], "satisfaction": k["satisfaction"], "period": len(v),
         "mean_consumption": round(mean, 2), "volatility": round(vol, 2),
         "zero_readings": int(min(v) == 0.0), "suspicious": int(z > 1.35)}
    x["per_person"] = round(x["mean_consumption"] / x["household"], 3)
    for b in REGION:
        x["area_" + b[0]] = int(k["region"][0] == b[0])
    DATA.append(x)


def shuffle(data, seed):
    r, s = rng(seed), list(range(len(data)))
    for i in range(len(s) - 1, 0, -1):
        j = int(r() * (i + 1))
        s[i], s[j] = s[j], s[i]
    return [data[i] for i in s]


K = shuffle(DATA, SEED + 90000)
TRAIN, VAL, TEST = K[:756], K[756:1008], K[1008:]
FIELDS = ["mean_consumption", "volatility", "household", "satisfaction", "period",
          "zero_readings", "per_person", "area_north", "area_south", "area_east",
          "area_west", "area_central"]
MEAN = {a: statistics.fmean(x[a] for x in TRAIN) for a in FIELDS}
STD = {a: max(1e-9, statistics.pstdev([x[a] for x in TRAIN])) for a in FIELDS}
BASELINE = sum(x["suspicious"] == 0 for x in TEST) / len(TEST)


def make_set(data):
    return ([[(x[a] - MEAN[a]) / STD[a] for a in FIELDS] for x in data],
            [float(x["suspicious"]) for x in data])


X_TR, Y_TR = make_set(TRAIN[:300])          # REDUCED training set
X_VA, Y_VA = make_set(VAL)
X_TE, Y_TE = make_set(TEST)


def squash(z):                              # S-shaped function, overflow-protected
    if z >= 0:
        return 1.0 / (1.0 + math.exp(-z)) if z < 700 else 1.0
    e = math.exp(z) if z > -700 else 0.0
    return e / (1.0 + e)


ACTIVATIONS = {                             # (function, derivative — its input is the activation output)
    "sigmoid": (squash, lambda a: a * (1.0 - a)),
    "hyperbolic tangent": (math.tanh, lambda a: 1.0 - a * a),
    "rectified": (lambda z: z if z > 0.0 else 0.0, lambda a: 1.0 if a > 0.0 else 0.0),
}
Z = (0.0, 1.0, 2.0, 4.0, 6.0, 8.0)
print(f"train {len(X_TR)}, validation {len(X_VA)}, test {len(X_TE)}, "
      f"baseline {BASELINE:.4f}")
print("\n" + f"{'derivative':<20}" + "".join(f"{'z=' + str(int(z)):>10}" for z in Z))
for name, (f, deriv) in ACTIVATIONS.items():
    print(f"{name:<20}" + "".join(f"{deriv(f(z)):>10.4f}" for z in Z))
for name, (f, deriv) in ACTIVATIONS.items():
    print(f"{name:<20} at z=1 {deriv(f(1.0)):.6f}, at z=6 {deriv(f(6.0)):.6f}, "
          f"ratio {deriv(f(1.0)) / max(deriv(f(6.0)), 1e-12):.1f}")

train 300, validation 252, test 252, baseline 0.7579

derivative                 z=0       z=1       z=2       z=4       z=6       z=8
sigmoid                 0.2500    0.1966    0.1050    0.0177    0.0025    0.0003
hyperbolic tangent      1.0000    0.4200    0.0707    0.0013    0.0000    0.0000
rectified               0.0000    1.0000    1.0000    1.0000    1.0000    1.0000
sigmoid              at z=1 0.196612, at z=6 0.002467, ratio 79.7
hyperbolic tangent   at z=1 0.419974, at z=6 0.000025, ratio 17088.4
rectified            at z=1 1.000000, at z=6 1.000000, ratio 1.0

The derivative of the sigmoid does not exceed 0.25 anywhere. This alone is a limit: the gradient coming from an input of twelve columns falls to one quarter at best as it passes the hidden layer. When the input rises to 6, the derivative becomes 0.0025 and the ratio reaches 79.7. The hyperbolic tangent is more generous at the ceiling — its derivative starts from 1 — but it saturates faster: the ratio between the same two points is 17,088. The derivative of the rectified unit is a constant 1 in the positive region and exactly 0 in the negative region.

Three Activations, Two Runs

The table alone does not say which is better, because in which region of the input the network actually works has not been measured. The same network is trained with three activations, from the same initial values, with two runs each.

def init_weights(seed, sizes):              # gain-scaled initialization, bias zero
    r, net = rng(seed), []
    for a, b in zip(sizes, sizes[1:]):
        scale = math.sqrt(2.0 / a)
        net.append(([[(r() * 2 - 1) * scale for _ in range(a)] for _ in range(b)], [0.0] * b))
    return net


def forward(net, x, f):                     # FORWARD PASS, output layer is always a sigmoid
    a, acts = x, [x]
    for k, (W, bias) in enumerate(net):
        z = [sum(w * v for w, v in zip(s, a)) + b for s, b in zip(W, bias)]
        a = [squash(z[0])] if k == len(net) - 1 else [f(t) for t in z]
        acts.append(a)
    return acts


def train_model(net, X, Y, rate, epochs, name):   # MODEL: full-batch gradient descent
    f, deriv = ACTIVATIONS[name]
    n, first = len(X), [0.0]
    for t in range(epochs):
        gW = [[[0.0] * len(W[0]) for _ in W] for W, _ in net]
        gb = [[0.0] * len(bias) for _, bias in net]
        for i in range(n):
            acts = forward(net, X[i], f)
            d = [acts[-1][0] - Y[i]]
            for k in range(len(net) - 1, -1, -1):
                W, a = net[k][0], acts[k]
                for j, dj in enumerate(d):
                    if dj:
                        s = gW[k][j]
                        for p in range(len(a)):
                            s[p] += dj * a[p]
                        gb[k][j] += dj
                if k:
                    prev = [0.0] * len(a)
                    for j, dj in enumerate(d):
                        if dj:
                            s = W[j]
                            for p in range(len(a)):
                                prev[p] += dj * s[p]
                    d = [prev[p] * deriv(a[p]) for p in range(len(a))]
        if t == 0:                          # mean gradient magnitude reaching the first layer
            first[0] = sum(abs(v) for s in gW[0] for v in s) / (len(gW[0]) * len(gW[0][0]))
        for k, (W, bias) in enumerate(net):
            for j in range(len(W)):
                s, g = W[j], gW[k][j]
                for p in range(len(s)):
                    s[p] -= rate * g[p] / n
                bias[j] -= rate * gb[k][j] / n
    return net, first[0]


def accuracy(net, X, Y, name):
    f = ACTIVATIONS[name][0]
    return sum((forward(net, x, f)[-1][0] > 0.5) == (y > 0.5) for x, y in zip(X, Y)) / len(X)


def unit_stats(net, X, name):               # active rows per unit and mean derivative
    f, deriv = ACTIVATIONS[name]
    W, bias = net[0]
    ceiling = 0.25 if name == "sigmoid" else 1.0
    counts, total = [0] * len(W), 0.0
    for x in X:
        for j in range(len(W)):
            t = deriv(f(sum(w * v for w, v in zip(W[j], x)) + bias[j]))
            total += t
            if t > 0.05 * ceiling:
                counts[j] += 1
    return counts, total / (len(X) * len(W))


RESULTS = {}
print(f"{'activation':<20} {'run':>6} {'validation':>10} {'test':>8} "
      f"{'mean_deriv':>10} {'first_layer_grad':>17}")
for name in ACTIVATIONS:
    RESULTS[name] = []
    for run in (1, 2):
        net, first = train_model(init_weights(SEED + run, [12, 8, 1]), X_TR, Y_TR, 0.3, 100, name)
        counts, mean_deriv = unit_stats(net, X_TR, name)
        RESULTS[name].append((accuracy(net, X_VA, Y_VA, name), accuracy(net, X_TE, Y_TE, name),
                              mean_deriv, first))
        print(f"{name:<20} {run:>6} {RESULTS[name][-1][0]:>10.4f} {RESULTS[name][-1][1]:>8.4f} "
              f"{mean_deriv:>10.4f} {first:>17.4f}")
CHOSEN = max(ACTIVATIONS, key=lambda a: RESULTS[a][0][0])
print(f"\nactivation chosen on validation: {CHOSEN}")
print(f"baseline {BASELINE:.4f}, network number {RESULTS[CHOSEN][0][1]:.4f}, "
      f"second run {RESULTS[CHOSEN][1][1]:.4f}")
print(f"first-layer gradient: {CHOSEN} {RESULTS[CHOSEN][0][3]:.4f}, rectified "
      f"{RESULTS['rectified'][0][3]:.4f}, ratio "
      f"{RESULTS['rectified'][0][3] / RESULTS[CHOSEN][0][3]:.2f} times")

activation              run validation     test mean_deriv  first_layer_grad
sigmoid                   1     0.8294   0.7937     0.1971            1.0456
sigmoid                   2     0.8333   0.7976     0.2100            1.5007
hyperbolic tangent        1     0.7976   0.7659     0.5913            3.1693
hyperbolic tangent        2     0.8016   0.7698     0.6587            4.8762
rectified                 1     0.8214   0.7857     0.5183            4.4498
rectified                 2     0.8056   0.7698     0.5012            5.2971

activation chosen on validation: sigmoid
baseline 0.7579, network number 0.7937, second run 0.7976
first-layer gradient: sigmoid 1.0456, rectified 4.4498, ratio 4.26 times

The gradient measure confirms the table: the mean gradient magnitude reaching the first layer is 1.0456 with the sigmoid and 4.4498 with the rectified unit — 4.26 times. This difference comes directly from the derivative ceiling and is the counterpart in the network of the saturation described in words.

But the test column says the opposite of what is expected. The saturating function gives the highest numbers: 0.7937 and 0.7976. The reason is in the mean derivative column: the mean derivative of the sigmoid units is 0.1971, 79% of its ceiling of 0.2500. The network is not working in the saturated region, because with scaled inputs the weighted sum stays small. Saturation is not a property of a function but a consequence of the region in which that function is run. At this scale a small gradient is not an obstacle; it will be one when the layers deepen and the same factor stacks up.

The triple is again in place: baseline 0.7579, network number 0.7937, second run 0.7976. The difference between the two runs is one subscriber; but the largest difference between the three activations is 0.7937 against 0.7659, that is, seven subscribers. In this lesson the choice produces a difference above the run-to-run noise.

Dead Unit and Saturated Unit

The way to tell the two forms of degradation apart is to increase the learning rate. As the step grows, the weights are thrown further and the units enter the region where the derivative is small.

print(f"{'activation':<20} {'rate':>6} {'test':>8} {'dead units':>10} "
      f"{'mean_deriv':>10}  active rows per unit")
for name in ACTIVATIONS:
    for rate in (0.3, 20.0, 30.0):
        net, _ = train_model(init_weights(SEED + 1, [12, 8, 1]), X_TR, Y_TR, rate, 100, name)
        counts, mean_deriv = unit_stats(net, X_TR, name)
        print(f"{name:<20} {rate:>6} {accuracy(net, X_TE, Y_TE, name):>8.4f} "
              f"{sum(1 for c in counts if c == 0):>10} {mean_deriv:>10.4f}  {counts}")

activation             rate     test dead units mean_deriv  active rows per unit
sigmoid                 0.3   0.7937          0     0.1971  [300, 300, 300, 300, 300, 300, 300, 300]
sigmoid                20.0   0.7579          0     0.0536  [176, 127, 118, 119, 132, 169, 141, 127]
sigmoid                30.0   0.6071          0     0.0295  [89, 68, 66, 70, 53, 95, 85, 71]
hyperbolic tangent      0.3   0.7659          0     0.5913  [296, 299, 291, 291, 283, 289, 288, 282]
hyperbolic tangent     20.0   0.7381          0     0.0528  [28, 37, 62, 27, 35, 36, 11, 47]
hyperbolic tangent     30.0   0.6667          0     0.0259  [5, 26, 17, 22, 15, 18, 15, 24]
rectified               0.3   0.7857          0     0.5183  [151, 185, 143, 142, 154, 134, 132, 203]
rectified              20.0   0.7579          4     0.0825  [0, 30, 91, 33, 44, 0, 0, 0]
rectified              30.0   0.7540          8     0.0000  [0, 0, 0, 0, 0, 0, 0, 0]

The distinction is in the last column. With the sigmoid and the hyperbolic tangent no unit dies: even at rate 30 every unit stays active on a part of the 300 rows, between 53 and 95. The mean derivative falls from 0.1971 to 0.0295 — a six-fold shrinkage — and the number drops to 0.6071, far below the baseline. But the derivative is not zero; the network can keep learning, even if slowly.

With the rectified unit the situation is different. At rate 20, four of the eight units die; at rate 30, all eight do, and the mean derivative becomes exactly 0.0000. The gradient flowing back through a dead unit is multiplied by zero; its weights are never updated again. This is an irreversible loss and it shows in the table as a number: when eight of eight units are dead, the output of the network is constant and the test number is 0.7540, that is, below even the baseline. With four units dead at rate 20, the number is exactly 0.7579, that is, the network tells everyone “not suspicious”.

Summary

• Saturation is a number: the derivative of the sigmoid is 0.1966 at z=1 and 0.0025 at z=6, a ratio of 79.7; for the hyperbolic tangent the same ratio is 17,088. The derivative of the rectified unit is a constant 1 in the positive region.
• In the same network the mean gradient reaching the first layer is 1.0456 with the sigmoid and 4.4498 with the rectified unit — 4.26 times. The difference comes directly from the derivative ceiling.
• The triple: baseline 0.7579, network number 0.7937 (sigmoid), second run 0.7976. The largest difference between the activations is seven subscribers, the difference between the two runs is one subscriber.
• Saturation is not a property of a function but a consequence of the region in which it is run: the mean derivative of the sigmoid units is 0.1971, that is 79% of its ceiling, because with scaled inputs the sum stays small.
• A dead unit and a saturated unit are different things. At rate 30 no unit dies with the sigmoid and the hyperbolic tangent, while with the rectified unit eight of eight die and the mean derivative becomes 0.0000; the number falls to 0.7540, below the baseline.

Next Step

In the two lessons the network always stayed the same: twelve inputs, eight units, a single output. How many units the hidden layer would carry was not chosen, and what would change if it were chosen was not measured. The next lesson opens the forward pass layer by layer and ties two separate numbers to this choice: how the parameter count and the number of multiplications made for a single example grow as the unit count increases, and how much the held-out set count responds to that growth. The matrix multiplication is written by hand, because the cost of a layer can be counted only when it is written out explicitly.