NMF-FFB Main Estimation Algorithm (formerly NMF-SEM)

Fits the NMF-FFB 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.ffb(
  Y1,
  Y2,
  rank = NULL,
  X.init = "nndsvd",
  X.L2.ortho = 100,
  C1.L1 = 1,
  C2.L1 = 0.1,
  epsilon = 1e-06,
  maxit = 5000,
  seed = 123,
  ...
)

nmf.sem(
  Y1,
  Y2,
  rank = NULL,
  X.init = "nndsvd",
  X.L2.ortho = 100,
  C1.L1 = 1,
  C2.L1 = 0.1,
  epsilon = 1e-06,
  maxit = 5000,
  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

Initialization strategy for the basis matrix X ($P_1 \times Q$). One of:

"nndsvd" (default): Non-negative Double SVD with additive randomness (NNDSVDar; Boutsidis & Gallopoulos 2008), computed internally via .nndsvdar(Y1, Q). Requires $Q \le \min(P_1, N)$ (over-rank case falls back to "runif"). Uses a full SVD of $Y_1$, so for very large $Y_1$ consider switching to "kmeans" to avoid SVD memory / compute cost.
"kmeans": k-means on the columns of $Y_1$ (samples clustered into $Q$ groups); the transposed cluster centers become $X$. Scales well for large $Y_1$; this is the default of nmfkc.
"kmeansar": "kmeans" followed by filling zero entries of $X$ with $\mathrm{Uniform}(0, \bar Y_1 / 100)$ (NNDSVDar-style additive randomness to escape trivial stationary points).
"runif": Uniform random entries in $[0, 1]$.
A numeric $P_1 \times Q$ matrix supplied by the user; negative entries are projected to 0.
NULL: backward-compatible alias for "nndsvd".

In all cases the result is column-normalized to colSums(X) = 1 before iteration. The menu mirrors nmfkc's X.init option for consistency across the package.

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: 5000 (matches nmfkc and other MU functions in the package).

seed

Random seed used to initialize X, C1, and C2. Default: 123.

...

Additional hidden arguments controlling the optional feedforward baseline (used both as an $X$ warm-start and as the reference for SC.map, the input-output structural fidelity defined in Satoh (2025) §4.SC.map):

nmfkc.baseline

Controls whether a feedforward nmfkc(Y1, A = Y2) fit is used as baseline. Possible values:

Default (not given) — nmf.sem runs nmfkc internally when X.init is a string method ("nndsvd", "kmeans", ...) or NULL, forwarding X.init, X.L2.ortho, epsilon, maxit, seed. The fitted $X$ of the baseline is then used as warm-start for the nmf.sem MU iterations, and SC.map is computed. This means nmf.sem(Y1, Y2, rank = Q) runs end-to-end without a prior nmfkc call.
TRUE — same as above, but force the internal nmfkc call even when X.init is a user-supplied matrix (the matrix is overridden).
FALSE — opt out; no internal call, SC.map = NA (pre-v0.6.8 behavior).
An nmfkc result (list with $X and $C) — use as the baseline directly (no internal call); also adopted as X.init when the latter is a string / NULL.

M.simple

Optional $P_1 \times P_2$ pre-computed baseline mapping. Takes precedence over nmfkc.baseline for the SC.map calculation but does not affect warm-start.

Q

Backward-compat alias for rank.

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$. See second-moment fidelity in Satoh (2025).
SC.map: Correlation between the equilibrium operator $M_{\mathrm{model}}$ and a feedforward baseline mapping $M_{\mathrm{simple}} = X_0 \Theta_0$, computed only when the baseline is supplied via M.simple or nmfkc.baseline in ...; otherwise NA. See input-output structural fidelity in Satoh (2025).
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-FFB 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)
#> Warning: maximum iterations (500) reached...
result$MAE
#> [1] 1.691285