Difference in python vs rust code

[collord]

Sometimes I see vastly different results from your Python implementation vs your Rust code. Here is a minimal synthetic example--am I doing something wrong?

import numpy as np
from esat.model.sa import SA
from esat.metrics import q_loss

rng = np.random.default_rng(0)
m, n, p = 50, 12, 3

# Synthetic data with known factorization
F_true = rng.uniform(0.1, 1.0, size=(p, n))
G_true = rng.uniform(0.1, 1.0, size=(m, p))
V = G_true @ F_true          # W @ H in ESAT convention
U = 0.001 + (0.05 * np.abs(V))

for label, use_rust in [
    ("Rust   (optimized=True) ", True),
    ("Python (optimized=False)", False),
]:
    sa = SA(V=V, U=U, factors=p, method="ls-nmf", seed=42, parallel=False, verbose=False)
    sa.optimized = use_rust
    sa.train(max_iter=20000, converge_delta=0.01, converge_n=100)

    W = sa.W.astype(np.float64)
    H = sa.H.astype(np.float64)
    Q = q_loss(V, U, W, H)
    rel_err = np.linalg.norm(V - W @ H) / np.linalg.norm(V)
    print(f"{label}:  Q={Q:.2f}  rel_err={rel_err:.4f}  converged={sa.converged}")

[AliWaseem607]

Hello, I was taking a look through the git and noticed your issue. I think there is an issue in your data setup. I was able to recreate your issue, but when I tried it with my own dataset I had the same result between using python or rust code. Looking through the differences between our setups, I noticed I was using esat's dataloader rather than manually creating U and V and that appears to fix the issue.

So there may just be different broadcasting rules between python and rust. Below is the code that I ran that shows the same results between python and rust on my machine.

import numpy as np
import pandas as pd
from esat.data.datahandler import DataHandler
from esat.metrics import q_loss
from esat.model.sa import SA

rng = np.random.default_rng(0)
m, n, p = 50, 12, 3

# Synthetic data with known factorization
F_true = rng.uniform(0.1, 1.0, size=(p, n))
G_true = rng.uniform(0.1, 1.0, size=(m, p))
V = G_true @ F_true  # W @ H in ESAT convention
U = 0.001 + (0.05 * np.abs(V))

data_handler = DataHandler.load_dataframe(
    input_df=pd.DataFrame(V), uncertainty_df=pd.DataFrame(U)
)
V, U = data_handler.get_data()

for label, use_rust in [
    ("Rust   (optimized=True) ", True),
    ("Python (optimized=False)", False),
]:
    sa = SA(
        V=V, U=U, factors=p, method="ls-nmf", seed=42, parallel=False, verbose=False
    )
    sa.initialize(H=None, W=None, init_method="k_means", init_norm=True)
    sa.optimized = use_rust
    sa.train(max_iter=20000, converge_delta=0.01, converge_n=100)

    W = sa.W.astype(np.float64)
    H = sa.H.astype(np.float64)
    Q = q_loss(V, U, W, H)
    rel_err = np.linalg.norm(V - W @ H) / np.linalg.norm(V)
    print(f"{label}:  Q={Q:.2f}  rel_err={rel_err:.4f}  converged={sa.converged}")

[collord]

I apologize for being a pain with what is probably just nonstandard array initialization. Someone would try this perhaps corner case at some point so the project might want to handle it?

This is reported by Claude, which sounds believable:

LSNMF.update — there is no epsilon clipping after the multiplicative update:

H_num = np.matmul(W.T, WeV)
H_den = np.matmul(W.T, np.multiply(We, WH))
H = np.multiply(H, np.divide(H_num, H_den))   # no clip
W = np.multiply(W, np.divide(W_num, W_den))   # no clip
Now trace the column_mean initialization:


V_avg = np.sqrt(VH_mean / self.factors)   # ≈ 0.48
H = V_avg * rng.standard_normal(...)      # ≈ N(0, 0.48)
H[H <= 0.0] = 1e-12                       # ~50% of entries clipped to 1e-12
W[W <= 0.0] = 1e-12                       # same

With ~50% of W and H entries sitting at 1e-12, WH ≈ 0 for many elements. H_den = W.T @ (We * WH) ≈ 0. On the first update step, H_num / H_den → ∞ for any row of H whose corresponding W column is dominated by 1e-12 entries. In float32, 1e-12 * inf = inf without any NaN, so the update silently produces astronomically large H values. Subsequent W update makes W gigantic too. Q = Σ (V - WH)²/U² returns a finite but absurd number (~122k) rather than NaN, so the Rust convergence check never triggers the "NaN, stopping" guard and happily reports converged=True.

Why Python survives: numpy float64 arithmetic in the Python path is more numerically forgiving — the same update eventually self-corrects over many iterations because the multiplicative rule is self-regularizing once at least one column of W has a non-tiny entry to anchor the ratio. The Rust path in float32 hits a one-way overflow (1e-12 × (large ratio) → inf → stays inf).

Why k-means fixes it: k-means centroids are at data scale (~0.5–1.0). The cluster-assignment W matrix has entries ≈ 1/p = 0.33. No entry is 1e-12 at initialization, so H_den and W_den are always well-conditioned and the Rust arithmetic never overflows.

The Rust bug is therefore a float32 overflow in the first few update steps when initialized with many near-zero (1e-12) W/H entries. The Python multiplicative update eventually recovers; the Rust one does not. The issue to file against ESAT is in esat_rust.ls_nmf: it needs epsilon-clipping of H and W after each update step, matching what initialize() does before step 1.