NMF-SEM Main Estimation Algorithm

Fits the NMF-SEM model $$ Y_1 \approx X \bigl( \Theta_1 Y_1 + \Theta_2 Y_2 \bigr) $$ under non-negativity constraints with orthogonality and sparsity regularization. The function returns the estimated latent factors, structural coefficient matrices, and the implied equilibrium (input–output) mapping.

At equilibrium, the model can be written as $$ Y_1 \approx (I - X \Theta_1)^{-1} X \Theta_2 Y_2 \equiv M_{\mathrm{model}} Y_2, $$ where $M_{\mathrm{model}} = (I - X \Theta_1)^{-1} X \Theta_2$ is a Leontief-type cumulative-effect operator in latent space.

Internally, the latent feedback and exogenous loading matrices are stored as C1 and C2, corresponding to $\Theta_1$ and $\Theta_2$, respectively.

Usage

nmf.sem(
  Y1,
  Y2,
  rank = NULL,
  X.init = NULL,
  X.L2.ortho = 100,
  C1.L1 = 1,
  C2.L1 = 0.1,
  epsilon = 1e-06,
  maxit = 20000,
  seed = 123,
  ...
)

Arguments

Y1: A non-negative numeric matrix of endogenous variables with rows = variables (P1), columns = samples (N).
Y2: A non-negative numeric matrix of exogenous variables with rows = variables (P2), columns = samples (N). Must satisfy ncol(Y1) == ncol(Y2).
rank: Integer; number of latent factors $Q$. If NULL, $Q$ is taken from a hidden argument in ... or defaults to nrow(Y2).
X.init: Optional non-negative initialization for the basis matrix X ($P_1 \times Q$). If supplied, it is projected to be non-negative and column-normalized.
X.L2.ortho: L2 orthogonality penalty for X. This controls the penalty term $\lambda_X \lVert X^\top X - \mathrm{diag}(X^\top X) \rVert_F^2$. Default: 100.
C1.L1: L1 sparsity penalty for C1 (i.e., $\Theta_1$). Default: 1.0.
C2.L1: L1 sparsity penalty for C2 (i.e., $\Theta_2$). Default: 0.1.
epsilon: Relative convergence threshold for the objective function. Iterations stop when the relative change in reconstruction loss falls below this value. Default: 1e-6.
maxit: Maximum number of iterations for the multiplicative updates. Default: 20000.
seed: Random seed used to initialize X, C1, and C2. Default: 123.
...: Additional arguments (reserved for future use).

Value

A list with components:

X: Estimated basis matrix ($P_1 \times Q$).
C1: Estimated latent feedback matrix ($\Theta_1$, $Q \times P_1$).
C2: Estimated exogenous loading matrix ($\Theta_2$, $Q \times P_2$).
XC1: Feedback matrix $X \Theta_1$.
XC2: Direct-effect matrix $X \Theta_2$.
XC1.radius: Spectral radius $\rho(X \Theta_1)$.
XC1.norm1: Induced 1-norm $\lVert X \Theta_1 \rVert_{1,\mathrm{op}}$.
Leontief.inv: Leontief-type inverse $(I - X \Theta_1)^{-1}.$
M.model: Equilibrium mapping $M_{\mathrm{model}} = (I - X \Theta_1)^{-1} X \Theta_2$.
amplification: Latent amplification factor $\lVert M_{\mathrm{model}} \rVert_{1,\mathrm{op}} / \bigl\lVert X \Theta_2 \bigr\rVert_{1,\mathrm{op}}$.
amplification.bound: Geometric-series upper bound $1 / (1 - \lVert X \Theta_1 \rVert_{1,\mathrm{op}})$ if $\lVert X \Theta_1 \rVert_{1,\mathrm{op}} < 1$, otherwise Inf.
Q: Effective latent dimension used in the fit.
SC.cov: Correlation between sample and model-implied covariance (flattened) of $Y_1$.
MAE: Mean absolute error between $Y_1$ and its equilibrium prediction $\hat Y_1 = M_{\mathrm{model}} Y_2$.
objfunc: Vector of reconstruction losses per iteration.
objfunc.full: Vector of penalized objective values per iteration.
iter: Number of iterations actually performed.

References

Satoh, K. (2025). Applying non-negative matrix factorization with covariates to structural equation modeling for blind input-output analysis. arXiv:2512.18250. https://arxiv.org/abs/2512.18250

Examples

# Simple NMF-SEM with iris data (non-negative)
Y <- t(iris[, -5])
Y1 <- Y[1:2, ]  # Sepal
Y2 <- Y[3:4, ]  # Petal
result <- nmf.sem(Y1, Y2, rank = 2, maxit = 500)
result$MAE
#> [1] 1.692369