End-to-End Joint NMF Worked Analysis

From Two Data Blocks to One Joint Analysis

A joint non-negative matrix factorization (joint NMF) seeks a compact set of sample profiles that can represent multiple omics blocks at once. Its defining choice is to use shared sample scores across blocks while allowing each block to have its own feature loadings. This makes the factors directly comparable across modalities, but it also means that preprocessing, scaling, rank, and fit diagnostics all matter to the result.

The analysis below uses two hypothetical omics blocks measured on the same four samples. The values are constructed for teaching and are not from a real study. They make it possible to inspect every reconstruction and residual while following the decisions that also arise in larger analyses.

Definition: Let \(X^{(m)}\) be a nonnegative matrix for omics block \(m\), with features in rows and the same \(n\) samples in columns. A rank-\(K\) joint NMF represents each block as \(X^{(m)} \approx W^{(m)}H^T\), where \(W^{(m)}\) is a nonnegative feature-loading matrix for that block and \(H\) is one nonnegative sample-score matrix shared across blocks. Each \(W^{(m)}\) has \(K\) columns, and \(H\) has \(n\) rows and \(K\) columns.

Before fitting, verify that columns refer to the same samples in the same order in every block. Decide how missing values, low-quality features, and transformations are handled; ordinary Frobenius-loss NMF requires finite, nonnegative entries. If a centering operation produces negative values, it is not suitable for this objective without a deliberate alternative formulation. Scaling the blocks is also substantive: a block with much larger numerical values can dominate a sum of squared errors. Use a stated scaling or block-weighting rule rather than assuming the two blocks contribute equally.

Specify the Objective Before Fitting

For squared-error fitting, a common weighted objective is

\min_{\{W^{(m)}\},H\geq 0} \frac{1}{2}\sum_{m=1}^{M}\lambda_m \left\|X^{(m)}-W^{(m)}H^T\right\|_F^2 .

Here \(\|\cdot\|_F^2\) is the sum of squared matrix entries, and \(\lambda_m\) is a prespecified nonnegative weight for block \(m\). The weights express the chosen relative contribution of each block after any preprocessing. They are not automatically estimates of biological importance. The notation \(\{W^{(m)}\},H\geq 0\) means that all factor entries are constrained to be nonnegative.

This objective assumes squared deviations are an appropriate measure of reconstruction error for the analyzed matrices. It does not, by itself, assert a probability model for the measurements. For count data or other data with a different error structure, a different loss may be more appropriate; the objective must match the data and scientific aim. The tutorial on NMF objectives for count matrices covers that alternative in more detail.

Justification for shared scores: If two blocks have the same sample order, the \(j\)th column of \(W^{(m)}H^T\) is \(W^{(m)}H_{j\cdot}^T\). Thus the same row \(H_{j\cdot}\) determines the component strengths for sample \(j\) in every block, while each block has different feature loadings \(W^{(m)}\). Sharing \(H\) therefore couples the reconstructions through samples, not through identical features.

The objective is not generally convex in all factors jointly. Alternating updates can reduce the loss, but a converged run need not be a global minimum. One common multiplicative update, shown to make the fitting mechanics explicit, is

W^{(m)} \leftarrow W^{(m)}\odot \frac{X^{(m)}H}{W^{(m)}H^TH}, \qquad H \leftarrow H\odot \frac{\sum_m\lambda_m(X^{(m)})^TW^{(m)}} {\sum_m\lambda_m H(W^{(m)})^TW^{(m)}} .

The symbol \(\odot\) denotes entrywise multiplication, and the division is entrywise. A small numerical safeguard is often used in denominators; zero entries require care because multiplicative updates can leave them at zero. These updates are an optimization procedure, not a guarantee of a unique or scientifically correct factorization. In practice, fit several initializations and retain the converged solution or solutions with the best objective under the same preprocessing and constraints.

Worked Example: Preparing and Aligning Two Blocks

Suppose four hypothetical samples, S1 through S4, have already been matched across two assays. After filtering and a stated nonnegative transformation, the matrices are

X^{(1)}= \begin{bmatrix} 4&0&4&8\\ 0&3&3&3\\ 1&1&2&3 \end{bmatrix}, \qquad X^{(2)}= \begin{bmatrix} 2&0&2&4\\ 0&5&5&5 \end{bmatrix}.

Rows are features and columns are S1, S2, S3, and S4 in both blocks. All entries are nonnegative, and neither matrix contains a missing value in this example. A valid joint analysis can therefore use the same four sample columns in both blocks. If the second matrix instead had columns ordered S1, S3, S2, S4, its columns would have to be reordered before fitting; otherwise the shared row of \(H\) would couple measurements from different samples.

For a transparent calculation, use the displayed values without additional rescaling and set \(\lambda_1=\lambda_2=1\). In a real analysis this is a modeling choice to justify, not a default that ensures equal influence: the blocks have different numbers of features and may have different scales. Any scaling and weights should be fixed before comparing candidate ranks.

Fit and Inspect a Candidate Factorization

For these matrices, consider rank \(K=2\), with

H= \begin{bmatrix} 1&0\\ 0&1\\ 1&1\\ 2&1 \end{bmatrix}, \quad W^{(1)}= \begin{bmatrix} 4&0\\ 0&3\\ 1&1 \end{bmatrix}, \quad W^{(2)}= \begin{bmatrix} 2&0\\ 0&5 \end{bmatrix}.

These factors are nonnegative. For sample S3, the score row is \((1,1)\). Its reconstruction in block 1 is \((4,0,1)(1)+(0,3,1)(1)=(4,3,2)\), and in block 2 it is \((2,0)(1)+(0,5)(1)=(2,5)\). This checks that both modality-specific reconstructions use the same sample scores.

Worked Example: Calculating the Reconstruction Loss

Using the proposed factors, block 1 is reconstructed as

W^{(1)}H^T= \begin{bmatrix} 4&0&4&8\\ 0&3&3&3\\ 1&1&2&3 \end{bmatrix}=X^{(1)}.

Likewise, \(W^{(2)}H^T=\begin{bmatrix}2&0&2&4\\0&5&5&5\end{bmatrix}=X^{(2)}\). Every residual entry is therefore zero. With \(\lambda_1=\lambda_2=1\), the objective value is

\frac{1}{2}\left( \sum_{i,j}(X^{(1)}_{ij}-(W^{(1)}H^T)_{ij})^2+ \sum_{i,j}(X^{(2)}_{ij}-(W^{(2)}H^T)_{ij})^2 \right) =\frac{1}{2}(0+0)=0.

A zero loss is possible here because the hypothetical matrices were constructed from the displayed rank-2 factors. In a real data set, exact reconstruction is unusual and not a target to pursue by increasing rank without limit. This calculation checks the model’s arithmetic; it does not establish that two components are biologically real or uniquely identified.

Assess Rank, Meaning, and Fit Together

Rank determines how many shared sample profiles are available. Training loss cannot increase when the feasible model is expanded, so a smaller loss at a larger rank is not, on its own, evidence that the extra component is useful. Rank assessment should combine reconstruction, stability across starts, interpretability, and—where feasible—performance on held-out entries or samples. Residual diagnostics should be considered at the selected rank, not replaced by a single scalar loss.

Worked Example: Comparing a Simple Rank-1 Candidate With Rank 2

To see why rank matters in the hypothetical data, consider a simple rank-1 candidate with every sample score equal to 1. For each feature, choose its across-sample mean as the loading; this gives the best constant reconstruction for that feature. The five loadings are 4, 2.25, 1.75, 2, and 3.75 for the three features of block 1 followed by the two features of block 2.

For the first feature of block 1, the residuals are \(4-4=0\), \(0-4=-4\), \(4-4=0\), and \(8-4=4\), giving squared error \(0+16+0+16=32\). For the remaining features, the squared errors are \(6.75\), \(2.75\), \(8\), and \(18.75\), respectively. Thus this particular rank-1 candidate has total squared error \(32+6.75+2.75+8+18.75=68.25\), or objective value \(68.25/2=34.125\) under equal block weights.

This is a concrete rank-1 reference, not a claim that it is the best rank-1 joint NMF: an optimized rank-1 solution can have an objective no larger than 34.125. The rank-2 factors above achieve objective zero for these constructed data. Rank 3 can also achieve zero, so training loss alone cannot distinguish ranks 2 and 3. Given the known construction in this example, rank 2 is the parsimonious choice. With real data, the choice would require repeat-fit stability and validation or other evidence, because the generating rank is not known.

Interpret the two columns of \(W^{(m)}\) within each modality and the corresponding columns of \(H\) together. Here component 1 has positive loading on feature 1 in both blocks, and component 2 has positive loading on feature 2 in both. S1 has scores \((1,0)\), S2 has \((0,1)\), and S3 has \((1,1)\). These descriptions identify patterns in the fitted representation; they do not establish that a component causes an outcome or corresponds to a particular biological mechanism. Feature names and independent subject-matter evidence are needed for that interpretation.

NMF factors have scale indeterminacy: multiplying a column of \(W^{(m)}\) by a positive constant and dividing the corresponding column of \(H\) by that constant leaves the reconstruction unchanged. Consequently, compare runs using a consistent normalization or a scale-invariant comparison. The earlier material on aligning components across NMF runs addresses the additional issue that component order can differ across runs.

Worked Example: Checking an Individual Residual

Suppose a fitted model predicts 2.7 for the third feature of block 1 in sample S4, while the observed value is 3. The residual, defined as observed minus fitted, is \(3-2.7=0.3\) in the transformed measurement units. Its squared contribution to the block’s unweighted Frobenius loss is \(0.3^2=0.09\). With \(\lambda_1=2\), its contribution to the full objective, including the factor \(1/2\), is \((1/2)(2)(0.09)=0.09\).

A small residual for this cell does not show that the whole fit is adequate. Inspect residuals across features and samples, and look for systematic structure: for example, a feature with consistently positive residuals, a sample with unusually large residuals in both blocks, or a modality with much larger weighted residual loss. Such patterns can indicate an inadequate rank, an unmodeled sample effect, an outlier, or a scaling and preprocessing problem. The earlier material on residual diagnostics for NMF gives further diagnostic tools.

Practical Checks and Limitations

A useful workflow is to record the sample alignment, transformations, feature filtering, scaling, weights, rank, initialization strategy, stopping rule, and final objective. Then check that the objective is finite and that all fitted factors are nonnegative; compare several starts; inspect blockwise and cellwise residuals; and assess whether the selected components recur and have coherent feature patterns. If an unusually good fit depends on one initialization, or residuals retain strong structure, report that uncertainty rather than presenting a single factorization as definitive.

Joint NMF is not always the best alternative. Separate NMF permits each block to have its own sample scores, which can be useful when cross-block patterns are weak, but it does not impose a shared representation. A count-appropriate loss may be preferable to squared error when measurement variability is not well described by squared deviations. A shared-score model is attractive when samples are genuinely matched and common profiles are scientifically plausible; it can be misleading when those conditions are not met. No factorization alone resolves batch effects, identifies causal mechanisms, or guarantees unique components.

Do not infer sample alignment from matching column positions. Verify sample identifiers and reorder columns explicitly.
Do not treat block weights as harmless constants. Changing them changes the fitted compromise between modalities.
Do not select rank only by training loss. Larger ranks can fit at least as well and may add unstable or redundant components.
Do not interpret scores without loadings. A high score has meaning only relative to the features contributing to that component and the scaling used.

Key takeaway: An end-to-end joint NMF analysis is a chain of linked decisions: align samples, prepare nonnegative matrices, define the weighted objective, fit shared scores and block-specific loadings, assess rank, interpret factors in context, and inspect residuals and stability. A low objective is evidence of reconstruction, not by itself evidence of validity.

Check Your Understanding

Answer each question, then check the solutions below.

In the joint model, which factor is shared across omics blocks, and which factors are block-specific?
Why must the sample columns be aligned before fitting?
For a residual of \(-0.4\) and block weight \(\lambda_m=1.5\), what is its contribution to the objective that includes the factor \(1/2\)?
Why does a zero training loss at rank 2 not establish that rank 2 is scientifically correct?
Name two checks, besides training loss, that can inform rank choice or fit adequacy.

Solutions. (1) \(H\), the sample-score matrix, is shared; each \(W^{(m)}\) is specific to one block. (2) A row of \(H\) is intended to describe the same sample in every block, so mismatched columns would couple different individuals. (3) The contribution is \((1/2)(1.5)(-0.4)^2=(1/2)(1.5)(0.16)=0.12\). (4) A larger rank can fit at least as well, and even an exact reconstruction does not establish stability, unique factors, or biological meaning. (5) Examples include repeated-start stability, held-out prediction, interpretability of loadings, and inspection of residual patterns.