feat: data-driven low-rank factor initialization for FactorEM

[jamgochiana]

Summary

Replaces the arbitrary 0.1*I initialization of the FactorEM low-rank factor F with a data-driven scheme built from each component's residual covariance.

Closes #17

The old scheme

initialize_gmm set F = zeros(n_features, rank), and m_step! detected the all-zero factor on the first inner iteration and replaced it with 0.1 * I. Problems: the 0.1 scale is arbitrary and data-scale-independent, it only populated the first rank features (identity columns), and it was not data-driven.

The new scheme

For each cluster/component, F is initialized to the leading rank principal directions of that cluster's residual covariance, scaled by the square root of the corresponding eigenvalues:

lambda, Q = eigen(Symmetric(residual_covariance))   # ascending
F = Q_top * Diagonal(sqrt.(max.(lambda_top, 0)))     # top  eigenpairs

This is exactly the low-rank factorization PCAEM already uses (F = P * Q * Diagonal(sqrt(lambda))). The diagonal D is set to the per-feature residual variance not explained by F (total_var - diag(F F')), floored at epsilon * mean(global_var) so D stays strictly positive (required by the LRDMvNormal constructor).

Edge cases handled:

• Empty clusters -> zero factor, global mean/variance.
• Clusters smaller than rank, or rank == 0 -> remaining F columns zero-padded.
• :rand method -> analogous eigendecomposition-based init from a random data subset (no longer the arbitrary identity).

The all(iszero, F) branch in m_step! is kept as a safety net for degenerate all-zero factors, but now scaled by the residual variance rather than a fixed 0.1.

Fitted results change by design

This changes fitted results (approved). The public API (FactorEM struct, fit signatures and return types) is unchanged.

Non-regression evidence

On fixed-seed synthetic data with genuine low-rank-plus-diagonal cluster structure (X = F_true*z + cluster_mean + diagonal noise, d=15, rank 3, 3 clusters):

At convergence (15 EM iters), new init matches old (delta ~ 0):

seed	old 0.1*I	new data-driven	delta
1	-17.566	-17.566	0.000
2	-18.739	-18.739	0.000
3	-18.353	-18.353	0.000
4	-17.029	-17.029	0.000
5	-18.313	-18.313	0.000

With few EM iterations the new init is strictly better (starts near the true structure):

nIter	mean delta (new - old) across 5 seeds
1	+6.2 nats
2	+0.8 nats
3	+0.2 nats

Tests

Added to test/test_fitting.jl (all deterministic, Random.seed!):

• Structure: fitted model has n_components LRDMvNormal comps, correct rank, finite logpdf, weights sum to 1.
• F-init sanity: after initialize_gmm, each component's F is finite and correctly shaped and D is strictly positive, across both methods and ranks {0, 2, 5}; factor is non-zero on structured data.
• Non-regression: converged mean log-likelihood on structured data reaches a good absolute value and is >= the old 0.1*I scheme (within noise); strictly greater with few iterations.

Full suite passes (191 fitting tests, all suites green). Unrelated GaussianMixtures "Too low occupancy" warnings are expected.

🤖 Generated with Claude Code