Matrix Factorization

The previous lesson built similarity directly from common observation: between two subscribers with no common item the similarity was not defined, and as sparsity rose, 105 of the 167 held-out subscribers were left with no neighbor. Yet two subscribers, one of whom applied for an installment plan and the other for debt restructuring, may carry the same billing behavior; the neighborhood routine cannot see this, because the only thing it looks at is the overlap of cells.

Matrix factorization removes this restriction: every subscriber and every item is represented by a vector on a small number of latent factors, and the prediction for a cell is the product of two vectors. Two subscribers with no common item can take high values in the same factor direction, so the similarity is read not from the cells but from the representation. The low-rank approximation of a full matrix was measured in the dimensionality reduction topic of the Unsupervised Learning course and is not repeated here. The only thing added is the missing entries, and it gives rise to a decision: of the 8,640 cells 8,230 are empty, and an empty cell either enters the computation as zero or does not enter at all. This lesson measures that decision.

• RS30. The setup, the split rule, the baseline and the measure carry on exactly from the previous two lessons: 360 subscribers, 24 items, 410 training interactions, 167 held-out subscribers, the measure hit rate@5, the baseline the order of the most chosen item.
• RS31. Factorization is the model and is written in the lesson with only the standard library: the subscriber matrix $P$ and the item matrix $Q$ are updated alternately, and at every step the penalized least squares is solved in closed form; the gradient descent of the Supervised Learning course is not repeated.
• RS32. The only difference between the two setups is which cells enter the sum in the update. In the zero-counting setup all 24 cells of the row enter and the target of the unobserved ones is 0. In the skipping setup only the observed cells enter and the unobserved ones are never written.
• RS33. The penalty coefficient is 1.0 in both setups and is not swept; what is measured is the difference made by the decision about the unobserved cell, not by the penalty coefficient.
• RS34. The number of factors $f$ is a hyperparameter; six candidates were tried (1, 2, 3, 4, 6, 8). The number of parameters is $384f$ and the observed interactions are 410; from $f=2$ on the number of parameters exceeds the number of observations.
• RS35. The alternating update runs for 10 rounds and the two setups start from the same initial point. The small square systems are solved by elimination; a threshold of $10^{-12}$ is written for the diagonal, because floating-point equality is not tested.
• RS36. Two errors are printed: the observed-cell error over the 410 filled cells and the all-cell error over all 8,640 cells, both root mean square deviation. Each setup is built to minimize one of them.
• RS37. The number of items shown and the share of the top-five lists falling on the five most shown items are the third number of the course; an item that is not shown does not enter the data of the next period.

One-Line Difference: Which Cells Enter the Computation

The first block builds the setup and the split in the form of the previous lessons; the numbers do not change.

# The setup catalog comes from the first lesson; here only the item's domain and base share are needed.
SEED, M32, N_SUBSCRIBERS = 20260218, 0xFFFFFFFF, 360
ITEMS = [("billing", 0.90), ("alert", 0.85), ("alert", 0.78), ("billing", 0.62),
         ("meter", 0.55), ("alert", 0.50), ("billing", 0.44), ("subscription", 0.40),
         ("meter", 0.36), ("alert", 0.32), ("billing", 0.28), ("subscription", 0.25),
         ("meter", 0.22), ("technical", 0.19), ("meter", 0.16), ("subscription", 0.14),
         ("technical", 0.12), ("subscription", 0.10), ("meter", 0.09), ("subscription", 0.08),
         ("technical", 0.07), ("billing", 0.06), ("alert", 0.05), ("technical", 0.04)]
SEGMENT = {
    "saver":    {"billing": 0.5, "meter": 1.7, "subscription": 0.6, "alert": 1.5, "technical": 0.6},
    "payer":    {"billing": 1.9, "meter": 0.5, "subscription": 0.6, "alert": 1.1, "technical": 0.3},
    "business": {"billing": 0.7, "meter": 0.9, "subscription": 1.6, "alert": 0.7, "technical": 1.9},
    "new":      {"billing": 0.8, "meter": 0.7, "subscription": 1.8, "alert": 1.3, "technical": 0.7}}
SEGMENT_NAMES, N_ITEMS = list(SEGMENT), len(ITEMS)


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


SUBSCRIBERS, INTERACTIONS = [], {}
for i in range(N_SUBSCRIBERS):
    r = rng(SEED + 7919 * i)
    a, segment = 20001 + i, SEGMENT_NAMES[int(r() * 4)]
    activity = 0.35 + 1.15 * r() ** 1.6
    SUBSCRIBERS.append(a)
    for j, (domain, share) in enumerate(ITEMS):
        propensity = min(0.95, SEGMENT[segment][domain] * share * activity * 0.55)
        if r() < min(0.95, share ** 0.55 * 0.72):
            if r() < propensity:
                INTERACTIONS[(a, j)] = r()
                if r() < 0.32:
                    r()
        else:
            r(), r(), r()

LOG, TRAIN, HELD_OUT = {}, {}, {}
for (a, j), stamp in INTERACTIONS.items():
    LOG.setdefault(a, []).append((stamp, j))
SEEN = {a: set() for a in SUBSCRIBERS}
for a, entries in LOG.items():
    entries.sort()
    if len(entries) >= 2:
        HELD_OUT[a] = entries[-1][1]
        entries = entries[:-1]
    for stamp, j in entries:
        TRAIN[(a, j)] = 1
        SEEN[a].add(j)

EVAL_SUBS = sorted(HELD_OUT)
N_CHOSEN = [sum(1 for (a, j) in TRAIN if j == jj) for jj in range(N_ITEMS)]
ITEM_SEEN = {j: set(a for a in SUBSCRIBERS if j in SEEN[a]) for j in range(N_ITEMS)}


def measure(score, k=5):
    hit, rank, shown = 0, 0, set()
    for a in EVAL_SUBS:
        s = score(a)
        candidates = sorted((j for j in range(N_ITEMS) if j not in SEEN[a]),
                            key=lambda j: (-s[j], j))
        shown.update(candidates[:k])
        p = candidates.index(HELD_OUT[a])
        hit, rank = hit + (p < k), rank + p + 1
    return hit / len(EVAL_SUBS), rank / len(EVAL_SUBS), len(shown)


print("training", len(TRAIN), "interactions -- empty cell", N_SUBSCRIBERS * N_ITEMS - len(TRAIN),
      "-- held-out subscribers", len(EVAL_SUBS))
print("baseline -- most chosen item: hit rate@5 %.4f, mean rank %.2f, shown %d/24"
      % measure(lambda a: N_CHOSEN))

training 410 interactions -- empty cell 8230 -- held-out subscribers 167
baseline -- most chosen item: hit rate@5 0.6527, mean rank 4.74, shown 9/24

The difference between the two setups is the single sum opened by the count_zero flag. When the flag is on, the square system is built over all items and is the same for every subscriber; when it is off, only over the observed items of the subscriber. The target vector is the same in both, because cells with a zero target add nothing to the sum.

def solve(A, b):                   # elimination on a small square system; f is at most 8
    n = len(b)
    M = [A[i][:] + [b[i]] for i in range(n)]
    for i in range(n):
        p = max(range(i, n), key=lambda t: abs(M[t][i]))
        M[i], M[p] = M[p], M[i]
        if abs(M[i][i]) < 1e-12:   # threshold: the penalty term keeps the diagonal away from zero
            continue
        d = M[i][i]
        for k in range(i, n + 1):
            M[i][k] /= d
        for t in range(n):
            if t != i and M[t][i]:
                c = M[t][i]
                for k in range(i, n + 1):
                    M[t][k] -= c * M[i][k]
    return [M[i][n] for i in range(n)]


def factorize(f, count_zero, rounds=10, penalty=1.0):
    # THE MODEL: matrix factorization by alternating updates, standard library only.
    r = rng(SEED + 5100 + f)
    P = {a: [0.1 + 0.2 * r() for _ in range(f)] for a in SUBSCRIBERS}
    Q = [[0.1 + 0.2 * r() for _ in range(f)] for _ in range(N_ITEMS)]
    blank = [[0.0] * f for _ in range(f)]
    for _ in range(rounds):
        G = ([[sum(q[x] * q[y] for q in Q) for y in range(f)] for x in range(f)]
             if count_zero else blank)
        for a in SUBSCRIBERS:
            idx = SEEN[a]
            A = [[G[x][y] + (0 if count_zero else sum(Q[j][x] * Q[j][y] for j in idx))
                  + (penalty if x == y else 0) for y in range(f)] for x in range(f)]
            P[a] = solve(A, [sum(Q[j][x] for j in idx) for x in range(f)])
        H = ([[sum(P[a][x] * P[a][y] for a in SUBSCRIBERS) for y in range(f)] for x in range(f)]
             if count_zero else blank)
        for j in range(N_ITEMS):
            users = ITEM_SEEN[j]
            A = [[H[x][y] + (0 if count_zero else sum(P[a][x] * P[a][y] for a in users))
                  + (penalty if x == y else 0) for y in range(f)] for x in range(f)]
            Q[j] = solve(A, [sum(P[a][x] for a in users) for x in range(f)])
    return P, Q


def predict(P, Q):
    return lambda a: [sum(P[a][x] * Q[j][x] for x in range(len(Q[0])))
                      for j in range(N_ITEMS)]


def errors(P, Q):
    f, so, sa = len(Q[0]), 0.0, 0.0
    for a in SUBSCRIBERS:
        for j in range(N_ITEMS):
            d = (1.0 if j in SEEN[a] else 0.0) - sum(P[a][x] * Q[j][x]
                                                     for x in range(f))
            sa += d * d
            if j in SEEN[a]:
                so += d * d
    return (so / len(TRAIN)) ** 0.5, (sa / (N_SUBSCRIBERS * N_ITEMS)) ** 0.5


def top5(score, a):
    s = score(a)
    return set(sorted((j for j in range(N_ITEMS) if j not in SEEN[a]),
                      key=lambda j: (-s[j], j))[:5])


for flag, label in ((True, "zero counted"), (False, "skipped     ")):
    P, Q = factorize(3, flag)
    oe, ae = errors(P, Q)
    print(f"f=3, unobserved {label}: hit rate@5 {measure(predict(P, Q))[0]:.4f}  "
          f"observed-cell error {oe:.4f}  all-cell error {ae:.4f}")

f=3, unobserved zero counted: hit rate@5 0.4731  observed-cell error 0.7168  all-cell error 0.1614
f=3, unobserved skipped     : hit rate@5 0.6826  observed-cell error 0.2125  all-cell error 0.5055

The observed-cell error of the skipping setup is 0.2125, that of the zero-counting one 0.7168, more than three times as much. In the all-cell error the order reverses: 0.1614 for the zero-counting setup, 0.5055 for the skipping one. Each setup reduces the sum it tries to reduce and loses under the other’s measure: the skipping setup sees no constraint on the unobserved cells and leaves a value around 0.5 there, while the zero-counting setup, pulling 8,230 empty cells to 0, finds no capacity to pull the 410 filled cells to 1.

The real number is the third. The setup with the smaller reconstruction error gives the worse ranking: 0.6826 for the skipping, 0.4731 for the zero-counting.

The Number of Factors Is a Hyperparameter

The number of factors is swept over six candidates; next to it are printed the number of parameters, the two errors, the number of items shown, the overlap with the baseline list and the overlap of the two setups with each other.

BASE5 = {a: top5(lambda x: N_CHOSEN, a) for a in EVAL_SUBS}
print(f"{'f':>2} {'par':>5} | {'A hit':>6} {'A obs':>7} {'A all':>6} {'A shw':>5} {'A ovl':>5}"
      f" | {'B hit':>6} {'B obs':>7} {'B all':>6} {'B shw':>5} {'B ovl':>5} | {'A~B':>5}")
for f in (1, 2, 3, 4, 6, 8):
    row, lists = [f"{f:>2} {(N_SUBSCRIBERS + N_ITEMS) * f:>5}"], []
    for flag in (True, False):
        P, Q = factorize(f, flag)
        sc = predict(P, Q)
        oe, ae = errors(P, Q)
        hit, rank, shw = measure(sc)
        L = {a: top5(sc, a) for a in EVAL_SUBS}
        lists.append(L)
        row.append(f"| {hit:>6.4f} {oe:>7.4f} {ae:>6.4f} {shw:>5} "
                   f"{sum(len(L[a] & BASE5[a]) for a in EVAL_SUBS) / len(EVAL_SUBS):>5.2f}")
    row.append(f"| {sum(len(lists[0][a] & lists[1][a]) for a in EVAL_SUBS) / len(EVAL_SUBS):>5.2f}")
    print(" ".join(row))

 f   par |  A hit   A obs  A all A shw A ovl |  B hit   B obs  B all B shw B ovl |   A~B
 1   384 | 0.6707  0.8377 0.1930    10  4.23 | 0.6707  0.2161 0.5207     9  4.34 |  4.02
 2   768 | 0.6048  0.7667 0.1760    11  4.16 | 0.6766  0.2124 0.5072     9  4.37 |  3.65
 3  1152 | 0.4731  0.7168 0.1614    12  3.43 | 0.6826  0.2125 0.5055     9  4.35 |  3.10
 4  1536 | 0.4251  0.6591 0.1479    15  2.68 | 0.7006  0.2115 0.5000     9  4.31 |  2.60
 6  2304 | 0.3593  0.5466 0.1263    20  1.45 | 0.6886  0.2113 0.4940     9  4.46 |  1.44
 8  3072 | 0.2395  0.4758 0.1066    20  0.82 | 0.6826  0.2107 0.4937     9  4.54 |  0.86

The last column is the headline of the lesson. With one factor the top-five lists of the two setups share 4.02 items, with eight factors 0.86. Same data, routine, penalty and starting point, and a one-line difference changes nearly all of the lists. The decision about the unobserved cell is not an implementation detail; it is the definition of what the model learns.

The zero-counting setup degrades with the number of factors: from 0.6707 to 0.2395. The reason is in the second column. The number of parameters is $384f$; with one factor it is 384, while the observed interactions are 410. From two factors on, the number of parameters exceeds the number of observations, and the model starts to carve out per-subscriber exceptions among the 8,230 zeros. The overlap with the baseline list falls from 4.23 to 0.82: the list breaks away from the baseline, but the place it breaks away to is not the right place.

The skipping setup is almost insensitive to the number of factors: its best is 0.7006 with four factors, 0.0479 above the baseline and 0.0060 above the best neighborhood setup of the previous lesson. The reason is structural: the unobserved cell is never looked at, so the extra factors have nothing to learn wrongly, and a subscriber’s ranking remains the items that lie close to their own items in the latent factor space. This is the item-based neighborhood of the previous lesson written over factors.

The Catalog Shown by the Two Rankings

Hit rate ranks the two setups in one direction; the items shown, the third number of the course, rank them in the other. The top-five lists of both setups are counted next to the baseline.

def distribution(score):
    counts = [0] * N_ITEMS
    for a in EVAL_SUBS:
        for j in top5(score, a):
            counts[j] += 1
    d = sorted(counts, reverse=True)
    return sum(1 for v in counts if v == 0), d[0], sum(d[:5]) / sum(d)


print(f"{'setup':>13} {'f':>2} {'never shown':>18} {'most shown':>18} "
      f"{'top-5 item share':>19}")
for label, f, flag in (("zero counted", 1, True), ("zero counted", 8, True),
                       ("skipped", 1, False), ("skipped", 8, False)):
    print("%13s %2d %18d %18d %19.4f"
          % ((label, f) + distribution(predict(*factorize(f, flag)))))
print("%13s %2s %18d %18d %19.4f" % (("baseline", "--") + distribution(lambda a: N_CHOSEN)))

        setup  f        never shown         most shown    top-5 item share
 zero counted  1                 14                150              0.7904
 zero counted  8                  4                 84              0.4467
      skipped  1                 15                146              0.7940
      skipped  8                 15                146              0.7856
     baseline --                 15                146              0.7868

The skipping setup shows the same catalog as the baseline at both numbers of factors: 15 items never appear, the most shown item reaches 146 of the 167 subscribers, and 0.7856 of the impressions fall on five items. The setup with the highest hit rate does not widen the visible part of the catalog at all; it amounts to reordering the baseline list a little.

With eight factors the zero-counting setup does the opposite: the items never shown go from 15 to 4, the number of subscribers reached by the most shown item from 146 to 84, and the share of the top five items from 0.7904 to 0.4467. Only counting the unobserved cell as zero pushes the tail of the catalog up, because the zeros of popular items also enter the sum and the direction of popularity is suppressed. The price is paid in hit rate: 0.2395, which is 0.4132 below the baseline.

No choice is made between the two numbers; choosing the setup determines both. An item that is not shown is not chosen, and an item that is not chosen does not enter the training set of the next period: a setup that never shows 15 items today has also decided which items tomorrow’s 410 interactions will come from.

Summary

• Matrix factorization places every subscriber and every item on a small number of latent factors; two subscribers with no common item can meet in the same factor direction. In a matrix with missing entries the single decision is what the unobserved cell is counted as, and the rest of the routine is exactly the same in the two setups.
• With three factors the observed-cell error of the skipping setup is 0.2125 and that of the zero-counting one 0.7168; in the all-cell error the order reverses (0.5055 and 0.1614). A small error does not mean a better ranking.
• The number of factors is a hyperparameter and six candidates were tried: the zero-counting setup degrades from 0.6707 to 0.2395, and the skipping setup gives its best with four factors (0.7006).
• The number of items shared by the top-five lists of the two setups is 4.02 with one factor and 0.86 with eight: two separate rankings from the same data.
• Coverage moves in the opposite direction: the skipping setup never shows the same 15 items as the baseline; the zero-counting setup lowers this to 4 and pulls the share of the top five items from 0.7868 to 0.4467.

Next Step

In four lessons the recommendation was always produced from the same place: from the existing interactions. The baseline counted the items chosen in training, the content-based setup extracted a profile from what the subscriber chose, the neighborhood looked for the common choices of two subscribers, and factorization placed the 410 filled cells on factors. The silent assumption of all four is that at least one interaction has come from the subscriber. Yet in the same setup 80 of the 360 subscribers have no interactions at all and the 15 items in the last table are shown to no subscriber; for these two groups the four methods produce the same thing. The next lesson takes up these two gaps by name, the subscriber with no interactions and the item never shown, and counts how many subscribers and how many items each solution brings into coverage and which measure it gives up in return.