model_selection() reports a cumulative running average across K, not each K's own held-out error

[mkb-phoenix]

Where: cp_utils.py, model_selection(), lines 108-110.

for k in n_topics:
    ...
    error = mean_squared_error(A.X.toarray(), A_pred)
    errors.append(error)
    error_mean = np.mean(errors)      # mean over EVERY k tried so far in the sweep, not just this k
    mean.append((k, error_mean))
...
print('The optimal number of factors for each rep is:', all_errors.idxmin())

Symptom: this is CellPie's own documented model-selection routine — the recommended way to choose the number of topics K via bi-cross-validation. The value recorded and plotted for a given K is not that K's own held-out reconstruction error; it's the running mean of every error computed so far across the whole K-sweep (errors is accumulated outside the loop and never reset). Because later terms in a running mean move less than earlier ones, the reported curve is smoothed toward monotonic improvement almost regardless of the true per-K trend, and idxmin() — the package's own suggested selection rule — is biased toward picking a K near the top of whatever range was tested.

Concretely, on our data: the buggy curve suggested K=16 for two different tissue regions. Hand-inverting the reported cumulative values back into each K's true isolated error (error_K = mean_K·i − mean_{K-1}·(i-1) for the i-th K in the sequence) recovered a real elbow at K=5-6 and K=7-9 for those same two regions — both fully masked by the cumulative-averaging artifact.

Suggested fix: record each K's own error only.

for k in n_topics:
    ...
    error = mean_squared_error(A.X.toarray(), A_pred)
    mean.append((k, error))           # no running average across k

Also worth flagging separately (not a code bug, a methodology gap): even fixed, this routine is a single 50/50 row+column split with no repeats and no seed-variance estimate. On our data a single split's noise step between adjacent K values was up to 165x larger than the margin that decided the reported "best K" — i.e. a single-fold BCV curve here isn't distinguishable from noise without a repeated/shuffled-split variant and some standard-error-based selection rule (we implemented multiple shuffled splits + a 1-SE rule on top of intNMF to work around this; happy to share that code if useful, but it's a separate script, not a patch to the package itself).

---

Found while running CellPie's own cp_utils.model_selection() bi-cross-validation at production scale (700K-1.5M spatial units per section). Happy to share our patched fork if useful (commit pinned at 8880d673c, main branch).