Why Joint NMF Results Can Mislead
Joint non-negative matrix factorization (joint NMF) can produce a compact representation of several omics blocks, but a converged algorithm does not guarantee a trustworthy interpretation. One block may dominate the objective because of its units or variability; a loss may not reflect the data-generating scale; an optimization run may settle in a poor local minimum; and factor magnitudes may look biologically absolute even when the model does not identify them that way. These errors can change which samples appear similar, which features seem important, and which components receive attention.
Consider \(M\) non-negative data matrices \(X^{(m)}\), one for each block \(m\), with the same \(n\) aligned samples and \(p_m\) features. A common joint NMF approximation uses a shared sample-score matrix \(W\) and a block-specific loading matrix \(H^{(m)}\):
Here \(K\) is the chosen rank, and \(\mathbb{R}_{+}\) denotes non-negative entries. One general joint objective is
where \(D_m\) measures discrepancy for block \(m\), and \(\lambda_m\) sets its contribution to the total objective. The matrices must represent the same samples in the same row order, entries must meet the chosen NMF method's non-negativity requirements, and the selected loss and weights must be specified. A small objective value only says that this particular weighted objective is small; it does not by itself establish that the model is suitable or that its components are scientifically meaningful.
Block-Scale Imbalance
With squared Frobenius loss, the block discrepancy is \(\|X^{(m)}-WH^{(m)}\|_F^2\), the sum of squared residuals across entries. If one block has much larger numerical values, units, or residual variance, it can dominate the objective even when the other block is equally important to the scientific question. A joint fit may then devote most of its flexibility to the numerically dominant block.
Worked Example: Finding a Dominant Block
Suppose a hypothetical analysis has two blocks, each with four entries, and the residual matrices at a candidate fit are
For an unweighted squared-error objective, block 1 contributes \(\|R^{(1)}\|_F^2=10^2+(-10)^2+10^2+(-10)^2=400\). Block 2 contributes \(\|R^{(2)}\|_F^2=1^2+(-1)^2+1^2+(-1)^2=4\). Thus block 1 contributes \(400/(400+4)=0.9901\), or about 99.0%, of the objective. A modeler who reports only that the total objective is 404 could miss how little this fit is being judged by block 2.
If block 1 was recorded in units ten times larger, dividing its values and corresponding fitted values by 10 would divide its squared residual contribution by 100, from 400 to 4. This shows why units and preprocessing matter; it does not establish that rescaling is always appropriate. The analysis should use transformations and block weights justified by measurement scale and scientific goals, then report how conclusions respond to reasonable alternatives.
Useful diagnostics include tabulating each weighted block contribution \(\lambda_m D_m\), comparing residual summaries on interpretable scales, and checking whether conclusions change when defensible normalization or weight choices are used. A weight selected only to make results look balanced can conceal rather than resolve a modeling problem. Likewise, equal weights do not imply equal influence when the discrepancy scales differ.
Loss Function Must Match the Measurement Scale
Squared Frobenius loss is natural when residuals are plausibly treated as continuous, approximately equal-variance deviations and large deviations should be penalized quadratically. It is not automatically suitable for every non-negative matrix. Count matrices often have mean-dependent variability; a count-oriented objective, such as generalized Kullback–Leibler divergence, can instead measure relative discrepancy in a way related to a Poisson likelihood. Such a likelihood interpretation requires an appropriate count sampling model, including its conditional mean structure; choosing the divergence alone does not prove that the data are Poisson.
Worked Example: Comparing Losses for Counts
Suppose one hypothetical count feature has observed values \((0,1,100)\), and a candidate fit gives \((2,3,98)\). The squared-error loss is
The generalized Kullback–Leibler loss is
The two totals are on different scales, so 12 cannot be declared worse than 2.922 by direct comparison. More importantly, they penalize discrepancies differently: squared error assigns a penalty of 4 to each of the three residuals, while the divergence assigns about 2 to the zero-versus-two discrepancy and about 0.0203 to the 100-versus-98 discrepancy. A count analysis should justify its loss using the count scale and its aims, and assess whether a Poisson-like variance assumption is credible. If counts are overdispersed or contain structural zeros, a simple Poisson interpretation may be inadequate.
A mismatch can be detected by plotting or summarizing residuals against fitted values and by checking whether the assumed variability is plausible across the fitted range. Persistent patterns, disproportionate influence from a few large observations, or systematic failures around zeros are warnings. Alternative losses are not interchangeable repairs: their objectives have different assumptions and meanings, and results should be compared using suitable diagnostics rather than raw objective values across losses.
Local Minima and Restart Sensitivity
Joint NMF objectives are generally non-convex in the factor matrices jointly. An iterative algorithm can converge to a local minimum or another stationary solution, so one initialization is not a reliable check of optimization quality. Before interpreting a solution, use multiple initializations under the same model, rank, loss, weights, and stopping rules. Compare final objective values and component patterns, accounting for possible component permutations and scaling when comparing runs.
Worked Example: A Single Run Hides a Better Fit
In a hypothetical two-block fit, three initializations use the same rank, objective, and convergence tolerance. Their final block contributions are:
| Initialization | Block 1 contribution | Block 2 contribution | Total objective |
|---|---|---|---|
| A | 82 | 28 | 110 |
| B | 49 | 31 | 80 |
| C | 52 | 30 | 82 |
Initialization A is worse under the specified objective because \(110>80\). Its termination is not evidence that it found the best solution available. Runs B and C are closer in objective, differing by \(82-80=2\), but that alone does not establish that their factors are equivalent or that 80 is the global minimum. The next checks are to run additional starts, verify convergence criteria, inspect block residuals, and compare the fitted component patterns after aligning components. If scientifically important patterns vary across near-optimal solutions, report that instability rather than choosing the most convenient run.
The comparison is meaningful only when the runs optimize the same objective on the same preprocessed data. A lower loss under different weights or a different divergence is not evidence that one initialization did better. Multi-start optimization is a diagnostic, not a proof of global optimality. It also does not address rank misspecification, poor preprocessing, or whether the loss represents the data adequately.
Factor Magnitudes Are Not Absolute Quantities
Even with a fixed fitted product, the factor magnitudes are not uniquely determined by the reconstruction objective. For any positive diagonal matrix \(D\),
Both transformed matrices remain non-negative, while the product is unchanged. Thus, without additional scale constraints or conventions, the size of a column of \(W\) cannot be interpreted as an absolute sample-level quantity independently of the corresponding row of every \(H^{(m)}\). Permuting component indices also leaves the reconstruction unchanged. Penalization or other constraints may change which transformations preserve the full objective, so the ambiguity should be assessed for the model actually fitted.
Worked Example: Same Reconstruction, Different Factor Sizes
For a single component, let the hypothetical score vector be \(W=(2,1)^\mathsf{T}\) and the loading vector be \(H=(3,4)\). Their product is
Choose \(c=10\). Replacing \(W\) by \(10W=(20,10)^\mathsf{T}\) and \(H\) by \(H/10=(0.3,0.4)\) gives \((10W)(H/10)=WH\), exactly the same reconstructed matrix. The first score vector has values ten times larger, despite no change in the fitted data. Reporting “a score of 20” as an absolute amount is therefore unjustified unless the model specifies and applies a normalization convention. A defensible report describes the convention and interprets scores and loadings relative to it.
Normalization can make numerical comparisons easier, but it is a reporting convention, not evidence that factor scale is intrinsically identified. When comparing runs, normalizing and aligning components can separate meaningful pattern differences from arbitrary scale and ordering differences. Do not compare raw loading magnitudes across analyses that used different scaling conventions as though they were directly commensurate.
A Practical Diagnostic Sequence
A disciplined review separates data checks, objective checks, optimization checks, and interpretation checks. Each addresses a different failure mode; passing one does not compensate for failing another.
Confirm sample alignment, non-negativity, preprocessing, rank, losses, and block weights. Record these choices so that fitted objectives can be compared meaningfully.
Calculate each \(\lambda_mD_m\) and inspect residual patterns within blocks. Investigate a block that dominates or is consistently poorly represented, rather than relying on the total alone.
Assess whether the chosen loss matches the measurement scale and plausible variability. Look for systematic residual structure, sensitivity to extreme values, and problems around zeros.
Use multiple starts with identical settings; compare objectives and aligned patterns. State any scaling convention before interpreting factor magnitudes.
The main alternatives address distinct needs. Block normalization or explicit weights can reduce undue scale dominance, but require a reasoned choice and sensitivity checks. A different loss may better suit counts or continuous measurements, but changes the target and assumptions. Multiple starts reduce dependence on a single initialization but cannot certify a global optimum. Explicit normalization resolves an arbitrary scale for presentation, but does not make an unconstrained factor magnitude an absolute biological measurement.
Common Mistakes and Limits
- “The total residual is small, so the joint fit is good.” A large block can conceal poor fit elsewhere. Report blockwise contributions and residual diagnostics.
- “The algorithm converged, so the solution is optimal.” Convergence can occur at a local solution. Compare multiple starts and state that they do not prove global optimality.
- “Counts are non-negative, so squared error is appropriate.” Non-negativity alone does not justify a loss. Match the objective to the scale and variability, and check its assumptions.
- “A larger score means a larger absolute biological quantity.” Rescaling the score and inversely rescaling the loadings leaves the reconstruction unchanged. Interpret magnitudes only with a stated convention.
These diagnostics cannot guarantee that the selected rank is correct, that preprocessing has removed unwanted technical effects, or that a component has a causal meaning. A joint NMF result is a model-based decomposition, not by itself proof of biological mechanism. Conclusions should reflect the fit, the assumptions, the sensitivity analyses, and the limits of the chosen representation.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Two blocks contribute 90 and 10 to an unweighted squared-error objective. What fraction of the objective comes from the first block, and what should be checked before deciding to reweight it?
- For a count observation \(x=0\) and fitted value \(y=4\), what is the generalized Kullback–Leibler contribution?
- Three starts under the same objective finish at 71, 73, and 96. What can and cannot be concluded from these values?
- If a component score column is multiplied by 5, how must its loading row change to preserve the reconstruction?
Solutions. (1) The fraction is \(90/(90+10)=0.90\), or 90%. Check units, preprocessing, residual scales, and the scientific role of both blocks; a dominant contribution is a warning, not by itself a prescription for weights. (2) By the zero convention, \(d(0\|4)=4\). (3) The run with objective 71 has the smallest observed objective among these starts, and 96 is substantially worse under the same objective. The three runs do not prove that 71 is a global minimum or that the two near-minimal runs have the same component patterns. (4) Divide the corresponding loading row by 5; then the score increase and loading decrease cancel in their product.