Regularized Estimation in High-Dimensional Factor Models

Why Regularize Factor Loadings?

A factor model represents each sample’s measurements using a smaller number of latent factors. In genomic data, the number of features $p$ may greatly exceed the number of samples $n$. An unrestricted fit can then give unstable estimates of the loading patterns: small changes in the samples may produce large changes in estimated loadings, especially when factors are weak or hard to distinguish. Regularization adds a constraint or penalty to discourage overly large or overly complex loading estimates.

The goal is not simply to make every loading small. Shrinkage and sparsity express different preferences: ridge shrinkage reduces loading magnitudes continuously, while an elementwise sparsity penalty can set some loadings exactly to zero. Those choices can alter the reconstructed profiles and, when factors are re-estimated, the fitted factors themselves.

Model, Scaling, and Estimation Objective

Definition: Let $X$ be an $n \times p$ matrix of centered measurements, with one row per sample and one column per feature. A $q$-factor model writes $$ X = F\Lambda^\mathsf{T} + E, $$ where $F$ is an $n \times q$ matrix of factor scores, $\Lambda$ is a $p \times q$ matrix of loadings, and $E$ is the residual matrix. The fitted matrix is $\hat X = \hat F\hat\Lambda^\mathsf{T}$. A row of $\Lambda$ describes how one feature relates to the $q$ factors.

A conventional probabilistic version assumes that the rows of $F$ have mean zero and covariance $I_q$, and that the residuals have mean zero, are independent of the factors, and have covariance matrix $\Psi$, often taken to be diagonal. A squared-error fit is especially natural under independent, mean-zero residuals with a common variance; unequal residual variances call for a weighted loss or a likelihood that accounts for them. The penalty does not repair an unsuitable residual model.

The factorization alone does not determine unique scores and loadings: rotating the factors and loadings together leaves $F\Lambda^\mathsf{T}$ unchanged. There is also a scale ambiguity: multiplying $F$ by a constant and dividing $\Lambda$ by that constant leaves the product unchanged. To make a loading penalty meaningful, fix the score scale. In the conditional fitting steps below, take $F^\mathsf{T}F/n=I_q$; when scores are updated, keep a stated normalization or convention. Factor signs and permutations can still vary without changing the fitted matrix.

For fixed, normalized scores, estimate each feature’s loading vector by minimizing a residual loss plus a penalty. With $x_j$ the $n$-vector of measurements for feature $j$, and $\lambda_j$ its $q$-vector of loadings, a ridge objective is $$ \frac{1}{2n}\|x_j-F\lambda_j\|_2^2 + \frac{\gamma}{2}\|\lambda_j\|_2^2, \qquad \gamma \geq 0. $$ The corresponding sparsity objective replaces the squared penalty by $\tau\|\lambda_j\|_1$, where $\tau \geq 0$ and $\|\lambda_j\|_1$ is the sum of the absolute loading values. Applying either penalty across all features sums the feature-specific objectives. These definitions describe conditional loading updates; a joint factor fit must also update $F$ under its scale constraint.

Two Penalties, Two Loading Effects

With $F^\mathsf{T}F/n=I_q$, define the unpenalized coefficient vector $b_j=F^\mathsf{T}x_j/n$. Differentiating the ridge objective with respect to $\lambda_j$ gives $$ -\frac{F^\mathsf{T}x_j}{n}+\frac{F^\mathsf{T}F}{n}\lambda_j+\gamma\lambda_j=0. $$ Thus $\hat\lambda_j=b_j/(1+\gamma)$. Ridge reduces every coefficient by the same positive factor in this orthonormal-score update, so it generally preserves which coefficients are zero and their relative signs, but it does not create exact zeros from nonzero coefficients.

For the sparsity penalty and orthonormal scores, the minimization separates by coefficient. Each estimate is the soft-thresholded value $$ \hat\lambda_{jk}=\operatorname{sign}(b_{jk})\max(|b_{jk}|-\tau,0). $$ A coefficient with magnitude at most $\tau$ becomes zero. This can give a simpler loading pattern, but the retained coefficients are not just uniformly rescaled: their relative magnitudes and the directions of the loading vectors can change.

Worked Example: Ridge Shrinkage for Two Features

Suppose there are $n=2$ samples and one fixed score vector $F=(1,-1)^\mathsf{T}$, so $F^\mathsf{T}F/n=1$. Two hypothetical centered feature vectors are $x_1=(3,-3)^\mathsf{T}$ and $x_2=(1,-1)^\mathsf{T}$. Their unpenalized coefficients are $$ b_1=\frac{F^\mathsf{T}x_1}{n}=\frac{6}{2}=3, \qquad b_2=\frac{F^\mathsf{T}x_2}{n}=\frac{2}{2}=1. $$

Use ridge penalty $\gamma=0.5$. For the first feature, $\hat\lambda_1=3/(1+0.5)=2$; for the second, $\hat\lambda_2=1/(1+0.5)=0.667$. With the scores held fixed, the fitted values are $F\hat\lambda_1=(2,-2)^\mathsf{T}$ and $F\hat\lambda_2=(0.667,-0.667)^\mathsf{T}$. Both loading magnitudes have decreased, and neither is zero. This conditional calculation isolates the loading update; it does not claim that the scores would remain fixed in a joint refit.

Worked Example: Thresholding a Weak Loading to Zero

Consider one fixed normalized score direction and two hypothetical features with unpenalized coefficients $b_1=0.8$ and $b_2=0.3$. Use the elementwise sparsity penalty $\tau=0.5$. Soft thresholding gives $$ \hat\lambda_1=\max(0.8-0.5,0)=0.3, \qquad \hat\lambda_2=\max(0.3-0.5,0)=0. $$

The first feature retains a positive loading, reduced from $0.8$ to $0.3$; the second contributes no fitted signal along this factor in this loading update. The zero is a consequence of the chosen penalty and data, not proof that the feature has no biological relationship to the underlying factor.

How Penalties Change Fitted Factors

For fixed scores, fitted values for feature $j$ are $F\hat\lambda_j$, so changing the loading changes that feature’s fitted profile. In a joint model, scores are also estimated from the loading patterns. Consequently, a penalty can change both the representation of features and the sample coordinates on the factors. A ridge penalty often retains a dense loading pattern but reduces its magnitude. A sparsity penalty can make the factors depend on fewer features, while changing the relative contributions of the features that remain.

The orthonormal-score convention matters when interpreting these changes. If all loadings for a factor are multiplied by a constant and scores are free to rescale inversely, the product $F\Lambda^\mathsf{T}$ can remain unchanged. A loading penalty by itself does not identify a useful scale in an unconstrained joint fit. The fitting procedure must enforce scale and rotation conventions, or include other identification restrictions. Otherwise, apparent shrinkage may reflect a change of scale rather than a change in the fitted signal.

Worked Example: Comparing Ridge and Sparse Fitted Profiles

Use the first hypothetical feature in the ridge example, whose unpenalized coefficient along the fixed score vector is $b=3$. Under ridge with $\gamma=0.5$, the loading is $2$, producing fitted values $(2,-2)^\mathsf{T}$. Now use an elementwise sparsity penalty $\tau=1.2$. Soft thresholding gives $\hat\lambda=3-1.2=1.8$, producing fitted values $(1.8,-1.8)^\mathsf{T}$. The sparse estimate is nonzero, but is smaller here than the ridge estimate.

For a second hypothetical feature with $b=1$, ridge still gives $1/(1+0.5)=0.667$, while the sparse estimate is $\max(1-1.2,0)=0$. Thus the ridge fit retains a smaller contribution for this feature, whereas the sparse fit omits it from the fitted profile. The comparison shows why the methods can yield different feature patterns; it does not establish that either pattern is the true biological mechanism.

Why the High-Dimensional Setting Matters

When $p>n$, the centered sample covariance matrix has rank at most $n-1$, so it cannot be inverted as a full-rank $p \times p$ matrix. Low-rank methods such as principal components can still extract a limited number of directions, but estimated directions may be sensitive to noise, especially when eigenvalues are close or the signal is weak. A penalty can stabilize the loading estimate by limiting its magnitude or complexity. It does not create information absent from the samples, and it cannot guarantee that the estimated factors correspond to distinct biological processes.

Ridge and sparsity are not interchangeable. Ridge is rotation-invariant for the loading matrix under an orthogonal rotation: the sum of squared loadings is unchanged. An elementwise absolute-value penalty is generally not rotation-invariant. Rotating factors can spread a sparse loading pattern across several factors, so a sparse solution depends on the chosen orientation and identification convention. In contrast, ridge tends to retain many features and may be preferable when weak contributions are plausible across a broad set of features. Sparse penalties are attractive when a compact loading pattern is scientifically useful, but their zeros can be unstable when correlated features can substitute for one another.

Fitting, Diagnostics, and Limitations

Joint estimation is generally not a single convex optimization: updating scores and loadings alternately can converge to different local solutions. A practical analysis should record the normalization and rotation convention, the penalty strength, and the initialization or fitting procedure. Check whether the algorithm converged, whether the fitted residuals show systematic structure, and whether loadings or fitted subspaces change substantially under reasonable perturbations of the data or initialization. For sparse fits, examine how often selected features change; for ridge fits, inspect whether shrinkage is so strong that meaningful factor signal is nearly removed.

A poor residual pattern can indicate an inadequate number of factors, unmodeled dependence, outliers, or a mismatch between the loss and residual variance structure. Penalty choice also involves a trade-off: stronger regularization may reduce variance but increase bias. No penalty ensures stable inference if the factors are not identifiable, the sample does not contain enough information, or the model assumptions fail. A sparse loading is a conditional feature-selection result under a particular model and tuning choice, not a causal claim.

Common mistake: “The lasso set this loading to zero, so the feature is unrelated to the factor.” The estimate is zero because its fitted contribution did not overcome the chosen threshold, conditional on the model, score orientation, and other included features. Correlated features may share or exchange signal, and another penalty or sample could retain this feature.

Another common mistake is to compare loading magnitudes without checking scale and orientation. Factor signs may flip, factors may be permuted, and rotations may change individual loadings while preserving the fitted matrix. Compare factors only after applying a consistent identification convention, and distinguish changes in loadings from changes in the reconstructed profiles.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

With normalized fixed scores, an unpenalized loading coefficient is $-2$. What is its ridge estimate when $\gamma=1$?
For an unpenalized coefficient $b=0.4$, what does soft thresholding with $\tau=0.6$ return?
Why does a loading penalty require a score-scale convention in a joint factor fit?
State one difference between ridge and elementwise sparsity penalties and one limitation of interpreting a zero sparse loading.

Solutions. (1) Ridge gives $-2/(1+1)=-1$, retaining the sign and halving the magnitude. (2) Soft thresholding gives $\max(0.4-0.6,0)=0$. (3) Without a scale constraint, multiplying scores by a constant and dividing loadings by it leaves the fitted matrix unchanged, so a loading penalty can be reduced by rescaling rather than by changing the fit. (4) Ridge generally shrinks coefficients without setting nonzero values exactly to zero, whereas sparsity can set coefficients to zero. A zero is conditional on the model and penalty; it does not establish that the feature has no biological relationship to the factor.