| Title: | Non-Negative Matrix Factorization with Kernel Covariates |
|---|---|
| Description: | Performs Non-negative Matrix Factorization (NMF) with Kernel Covariates. Given an observation matrix and kernel covariates, it optimizes both a basis matrix and a parameter matrix. Notably, if the kernel matrix is an identity matrix, the method simplifies to standard NMF. Also provides NMF with Random Effects (NMF-RE) via nmfre(), which estimates a mixed-effects model combining covariate-driven scores with unit-specific random effects together with wild bootstrap inference, and NMF with feed-forward and feedback structure (NMF-FFB) via nmf.ffb(), which fits a two-block input-output model for blind source separation and path analysis. References: Satoh (2025) <doi:10.48550/arXiv.2403.05359>; Satoh (2026) <doi:10.1007/s42081-026-00349-x>; Satoh (2025) <doi:10.48550/arXiv.2512.18250>; Satoh (2026) <doi:10.48550/arXiv.2603.01468>; Satoh and Tokuda (2026) <doi:10.48550/arXiv.2607.27474>; Satoh (2026) <doi:10.1007/s42081-025-00314-0>. |
| Authors: | Kenichi Satoh [aut, cre] (ORCID: <https://orcid.org/0000-0003-4436-9347>) |
| Maintainer: | Kenichi Satoh <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.0.2 |
| Built: | 2026-10-29 22:27:55 UTC |
| Source: | https://github.com/ksatohds/nmfkc |
Returns the coefficients data frame from a fitted NMF model
that has been passed through an inference function
(nmfkc.inference, nmf.rrr.inference,
nmfre.inference).
If inference has not been run, returns the parameter matrix
() directly.
For nmf.ffb objects, returns (exogenous block) as fallback.
## S3 method for class 'nmf' coef(object, ...) ## S3 method for class 'nmf.ffb' coef(object, ...)## S3 method for class 'nmf' coef(object, ...) ## S3 method for class 'nmf.ffb' coef(object, ...)
object |
A fitted model object of class |
... |
Not used. |
A data frame of coefficients (if inference was performed),
or the parameter matrix .
nmfkc.inference, nmf.rrr.inference,
nmfre.inference, nmf.ffb.inference
Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) coef(result) # returns C matrix result2 <- nmfkc.inference(result, Y, A) coef(result2) # returns coefficients data frameY <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) coef(result) # returns C matrix result2 <- nmfkc.inference(result, Y, A) coef(result2) # returns coefficients data frame
Extract the covariate-coefficient matrix from an NMF-GMM fit
## S3 method for class 'nmf.gmm' coef(object, ...)## S3 method for class 'nmf.gmm' coef(object, ...)
object |
An object of class |
... |
Ignored. |
The coefficient matrix ().
Part of the experimental NMF-GMM family: see nmf.gmm.
The interface may change in future versions — argument names, defaults
and the contents of the returned object are not yet stable.
nmf.gmm, nmf.gmm.inference, nmf.gmm.select
Returns the reconstructed matrix from a fitted
NMF model.
For nmf.ffb objects, returns the equilibrium prediction
if available. Supply Y1 and
Y2 to get the direct reconstruction
instead.
## S3 method for class 'nmf' fitted(object, ...) ## S3 method for class 'nmfae' fitted(object, ...) ## S3 method for class 'nmfre' fitted(object, type = c("blup", "fixed"), ...) ## S3 method for class 'nmf.ffb' fitted(object, ...)## S3 method for class 'nmf' fitted(object, ...) ## S3 method for class 'nmfae' fitted(object, ...) ## S3 method for class 'nmfre' fitted(object, type = c("blup", "fixed"), ...) ## S3 method for class 'nmf.ffb' fitted(object, ...)
object |
A fitted model object of class |
... |
For |
type |
For |
The fitted matrix .
nmfkc, nmf.rrr, nmfre,
nmf.ffb, residuals.nmf
result <- nmfkc(matrix(runif(50), 5, 10), rank = 2) fitted(result)result <- nmfkc(matrix(runif(50), 5, 10), rank = 2) fitted(result)
Fitted (responsibility-averaged) reconstruction of an NMF-GMM fit
## S3 method for class 'nmf.gmm' fitted(object, ...)## S3 method for class 'nmf.gmm' fitted(object, ...)
object |
An object of class |
... |
Ignored. |
The fitted matrix .
Part of the experimental NMF-GMM family: see nmf.gmm.
The interface may change in future versions — argument names, defaults
and the contents of the returned object are not yet stable.
nmf.gmm, nmf.gmm.inference, nmf.gmm.select
Computes the clustering-quality criteria silhouette,
CPCC, and dist.cor for a list of models fitted at
different ranks (or a single fit), returning one row per rank. These
are clustering-stability diagnostics (how decisively and
faithfully the samples cluster), conceptually separate from the
rank-selection *.rank functions (which use r.squared, effective
rank, and ECV) and complementary to nmf.cluster.flow
(which shows how the hard clustering itself changes across ranks).
Hard sample clustering requires a non-negative coefficient/score
matrix (so the columns form a membership simplex); when a model's
coefficient is signed (e.g.\ nmfkc.signed, nmfae.signed,
nmfre fits whose coefficient has negative entries) the
hard-label silhouette is NA while the distance-based
CPCC and dist.cor are still computed.
nmf.cluster.criteria(fits, Y, Y2 = NULL, names = NULL, plot = TRUE, ...)nmf.cluster.criteria(fits, Y, Y2 = NULL, names = NULL, plot = TRUE, ...)
fits |
A list of fitted models, one per rank, all over the same
|
Y |
The original data matrix used to fit the models ( |
Y2 |
Exogenous block, required only for |
names |
Optional character vector (length |
plot |
Logical; draw the diagnostics plot immediately
(default |
... |
When |
An object of class "nmf.cluster.criteria" (returned
invisibly): a list with criteria (a data frame with one row
per result and columns rank, silhouette, CPCC,
dist.cor, and hard) and labels (the x-axis
labels). Results are kept in the given order (not sorted).
Y <- t(as.matrix(iris[, 1:4])) fits <- lapply(2:6, function(q) nmfkc(Y, Q = q, print.dims = FALSE)) cc <- nmf.cluster.criteria(fits, Y, plot = FALSE) cc$criteria plot(cc)Y <- t(as.matrix(iris[, 1:4])) fits <- lapply(2:6, function(q) nmfkc(Y, Q = q, print.dims = FALSE)) cc <- nmf.cluster.criteria(fits, Y, plot = FALSE) cc$criteria plot(cc)
Visualizes how the hard sample clustering changes across a sequence of
fitted models – typically the same model at increasing ranks, but
also different models at the same rank. Each individual is a line
flowing left-to-right across the results (x-axis); its vertical
position at each result is determined by its cluster, and clusters are
reordered (barycenter method) to reduce crossings. Lines are coloured
by the individual's cluster in the reference result, so one can
see how the reference clusters split or merge. The adjusted Rand
index (ARI) between each pair of adjacent results is printed along the
top of the figure. X-axis ticks default to each result's $rank
and can be overridden with names.
Works for any non-negative multiplicative-update family
(nmfkc, nmfae, nmfkc.net, nmfre,
and the signed variants); the hard label is the argmax of the
coefficient/score matrix.
nmf.cluster.flow(fits, reference = NULL, names = NULL, plot = TRUE, ...)nmf.cluster.flow(fits, reference = NULL, names = NULL, plot = TRUE, ...)
fits |
A list (length |
reference |
The index (1-based position in |
names |
Optional character vector (length |
plot |
Logical; draw the diagram immediately by calling
|
... |
When |
An object of class "nmf.cluster.flow" (returned
invisibly): a list with clusters (the table:
rows = individuals, columns = results, entries = cluster number = the
dominant-factor index of each fit, so it matches the factor/basis
numbering of fits; a factor that never dominates leaves an
empty, unused cluster number),
ypos (the layout positions), ranks (each result's
rank), labels (the x-axis labels), reference (the
reference index), ref.cluster (the reference hard labels),
ARI (adjusted Rand index between each pair of adjacent
results, length ), and colors
(the default per-individual reference colour). Call
plot on it to (re)draw the diagram.
plot.nmf.cluster.flow,
nmf.cluster.criteria, nmfkc.rank
Y <- t(as.matrix(iris[, 1:4])) fits <- lapply(2:6, function(q) nmfkc(Y, Q = q, print.dims = FALSE)) fl <- nmf.cluster.flow(fits, reference = 2, plot = FALSE) # 2nd result head(fl$clusters) plot(fl, lwd = 2, main = "iris cluster flow")Y <- t(as.matrix(iris[, 1:4])) fits <- lapply(2:6, function(q) nmfkc(Y, Q = q, print.dims = FALSE)) fl <- nmf.cluster.flow(fits, reference = 2, plot = FALSE) # 2nd result head(fl$clusters) plot(fl, lwd = 2, main = "iris cluster flow")
Fits the NMF-FFB model
in which a non-negative basis () generates the
endogenous block from latent scores that are driven by
the exogenous block (feed-forward, ) and fed back
by itself (latent feedback, ). The function
returns the estimated basis, the structural coefficient matrices and the
implied equilibrium (input-output) mapping.
At equilibrium, the model can be written as
where 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 and ,
respectively.
Two estimators are available through method:
"fiml" (default)A two-stage likelihood-based estimator.
Stage 1 estimates the basis with the feed-forward fit
nmfkc(Y1, A = Y2) (X.init, X.L2.ortho,
epsilon, maxit, seed are forwarded), unless a
basis is supplied through X. Columns are normalized to unit
sum. Stage 2 holds fixed and fits the Gaussian
working model , by full-information maximum likelihood
(L-BFGS-B with an analytic gradient) under the non-negativity of
and an exclusion restriction on
. It fits (a) the feed-forward null
(a non-negative MIMIC factor model with correlated factors), (b) the
unpenalized feedback model on the admitted entries, and (c) an L1
path over C1.L1.path. The penalized problem is non-convex, so
every point of the path is fitted from several starting points
(starts); every distinct support proposed by any (penalty,
start) pair is re-estimated without penalty (itself from three
starts, keeping the best log-likelihood) and BIC is minimized over
all distinct candidates together with the null and the unpenalized
model. The support with the smallest BIC is the selected
model, reported in C1, C2, Phi, psi.
C1.L1 and C2.L1 are not used by this estimator. The likelihood-ratio statistics against the null are
returned without p-values: puts the null on the
boundary of the parameter space and the BIC refit is a post-selection
statistic, so a chi-square reference is not valid. Their calibration
by parametric bootstrap is done by nmf.ffb.inference.
"mu"The legacy estimator: joint multiplicative updates of
for the structural-form squared error
with an
orthogonality penalty on and L1 penalties C1.L1,
C2.L1. Kept, bit-identical, so that analyses published with it
reproduce. Because the structural and reduced forms have the same
fit once is free, this estimator cannot separate
from ; prefer "fiml".
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, ..., method = c("fiml", "mu"), X = NULL, C1.restriction = c("union", "none"), C1.restriction.threshold = 0.05, Phi.restriction = c("full", "diag"), C1.L1.path = NULL, select = c("BIC", "none") )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, ..., method = c("fiml", "mu"), X = NULL, C1.restriction = c("union", "none"), C1.restriction.threshold = 0.05, Phi.restriction = c("full", "diag"), C1.L1.path = NULL, select = c("BIC", "none") )
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 |
rank |
Integer; number of latent factors |
X.init |
Initialization strategy for the basis matrix
To impose anchor rows on the stage-1 basis instead, pass
|
X.L2.ortho |
L2 orthogonality penalty for |
C1.L1 |
L1 sparsity penalty for |
C2.L1 |
L1 sparsity penalty for |
epsilon |
Relative convergence threshold for the objective function.
Iterations stop when the relative change in reconstruction loss
falls below this value. Default: |
maxit |
Maximum number of iterations for the multiplicative updates.
Default: |
seed |
Random seed used to initialize |
... |
Additional hidden arguments. For
|
method |
|
X |
Optional basis for |
C1.restriction |
Exclusion restriction on |
C1.restriction.threshold |
Loading threshold for
|
Phi.restriction |
Covariance of the latent disturbance |
C1.L1.path |
Numeric vector of L1 penalties on |
select |
|
An object of class c("nmf.ffb", "nmf"), a list
with components:
X |
Estimated basis matrix ( |
C1 |
Estimated latent feedback matrix ( |
C2 |
Estimated exogenous loading matrix ( |
XC1 |
Feedback matrix |
XC2 |
Direct-effect matrix |
XC1.radius |
Spectral radius |
XC1.norm1 |
Induced 1-norm |
Leontief.inv |
Leontief-type inverse |
M.model |
Equilibrium mapping
|
amplification |
Latent amplification factor
|
amplification.bound |
Geometric-series upper bound
|
rank |
Effective latent dimension used in the fit. |
SC.cov |
Correlation between sample and model-implied covariance
(flattened) of |
SC.map |
Correlation between the equilibrium operator
|
mae |
Mean absolute error between |
objfunc |
Vector of reconstruction losses per iteration
( |
objfunc.penalized |
Vector of penalized objective values per iteration
( |
iter, maxit, epsilon, converged
|
Convergence bookkeeping. For
|
method |
|
The following are present for method = "fiml" only
(SC.cov and SC.map are then NULL):
Phi, psi, loglik, npar
|
Latent-disturbance covariance
( |
null |
The feed-forward null: list |
full |
The unpenalized feedback fit: list |
path |
Data frame with one row per ( |
candidates, supports, support.selected
|
One row per distinct
support ( |
C1.free, C1.L1.path, C1.L1.selected, support
|
The free-entry matrix used
( |
LR, LR.df
|
Likelihood-ratio statistics |
BIC, AIC
|
Named vectors |
call |
The matched call, from which
|
Feedback is identified only through exclusion restrictions: an outcome
may not feed back into a factor on which it has more than a negligible
loading, because such an entry is nearly equivalent to a change of the
outcome's loading and uniqueness. With
C1.restriction = "union" (default) entry of is
excluded if (the dominant factor) or
C1.restriction.threshold. Both halves are needed:
blocking the dominant factor alone would leave an outcome with a
substantial second loading free to feed that factor, and blocking only the
factors above the threshold would leave an outcome whose largest loading is
below the threshold free to feed its own.
"none" frees every entry (not recommended: the model is then
identified only through the non-negativity and the covariance structure).
A user-supplied 0/1 matrix is used as given.
The restriction is derived from the estimated basis and therefore from the same
that is subsequently tested; see the Calibration section
of nmf.ffb.inference for what this implies.
method = "fiml" became the default in version 0.9.8, as did
C1.restriction = "union" (earlier fiml fits blocked the dominant
factor only), and
nmf.ffb.inference gained the calibration argument.
method = "mu" is the legacy estimator, kept for the reproducibility
of published analyses; it will be deprecated in a later release.
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
nmf.ffb.inference, nmf.ffb.cv,
nmf.ffb.split, nmf.ffb.DOT,
summary.nmf.ffb
# Simple NMF-FFB with iris data (non-negative) Y <- t(iris[, -5]) Y1 <- Y[1:2, ] # Sepal Y2 <- Y[3:4, ] # Petal result <- nmf.ffb(Y1, Y2, rank = 2) result$LR # feedback vs feed-forward null (no p-value here) result$BIC result$path # Legacy multiplicative-update estimator result.mu <- nmf.ffb(Y1, Y2, rank = 2, maxit = 500, method = "mu") result.mu$mae# Simple NMF-FFB with iris data (non-negative) Y <- t(iris[, -5]) Y1 <- Y[1:2, ] # Sepal Y2 <- Y[3:4, ] # Petal result <- nmf.ffb(Y1, Y2, rank = 2) result$LR # feedback vs feed-forward null (no p-value here) result$BIC result$path # Legacy multiplicative-update estimator result.mu <- nmf.ffb(Y1, Y2, rank = 2, maxit = 500, method = "mu") result.mu$mae
Performs K-fold cross-validation to evaluate the equilibrium mapping of the NMF-FFB model.
For each fold, nmf.ffb is fitted on the training samples,
yielding an equilibrium mapping .
The held-out endogenous variables are then predicted from
using this mapping, and the mean absolute error (MAE) over all entries in the
test block is computed. The returned value is the average MAE across folds.
This implements the hyperparameter selection strategy described in the paper: hyperparameters are chosen by predictive cross-validation rather than direct inspection of the internal structural matrices.
method = "fiml". When method = "fiml" is passed
(through ...), the column-wise scheme above is not run: the
equilibrium prediction is a reduced-form
quantity, and the reduced form is the same with and without feedback, so
predicting held-out columns cannot select (that is what the
BIC path and the bootstrap-calibrated LR test in nmf.ffb /
nmf.ffb.inference are for). The only tunable quantity is the
stage-1 rank, and the call is delegated to the element-wise CV of the
feed-forward fit, nmfkc.ecv(Y = Y1, A = Y2, rank = rank,
...), whose object is returned ($objfunc holds the CV error per
rank; rank may be a vector, default 1:min(3, P1)).
C1.L1 and C2.L1 are ignored in that case.
nmf.ffb.cv( Y1, Y2, rank = NULL, X.init = "nndsvd", X.L2.ortho = 100, C1.L1 = 1, C2.L1 = 0.1, epsilon = 1e-06, maxit = 5000, ... )nmf.ffb.cv( Y1, Y2, rank = NULL, X.init = "nndsvd", X.L2.ortho = 100, C1.L1 = 1, C2.L1 = 0.1, epsilon = 1e-06, maxit = 5000, ... )
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 |
rank |
Integer; rank (number of latent factors) passed to |
X.init |
Initialization strategy for |
X.L2.ortho |
L2 orthogonality penalty for |
C1.L1 |
L1 sparsity penalty for |
C2.L1 |
L1 sparsity penalty for |
epsilon |
Convergence threshold for |
maxit |
Maximum number of iterations for |
... |
Additional arguments passed to |
A numeric scalar: mean MAE across CV folds (the fits use
method = "mu" unless method is given). With
method = "fiml", the nmfkc.ecv object of the
stage-1 rank sweep.
Y <- t(iris[, -5]) Y1 <- Y[1:2, ] Y2 <- Y[3:4, ] mae <- nmf.ffb.cv(Y1, Y2, rank = 2, maxit = 500, nfolds = 3) mae # rank selection for the likelihood-based estimator ecv <- nmf.ffb.cv(Y1, Y2, rank = 1:2, method = "fiml", nfolds = 3) ecv$objfuncY <- t(iris[, -5]) Y1 <- Y[1:2, ] Y2 <- Y[3:4, ] mae <- nmf.ffb.cv(Y1, Y2, rank = 2, maxit = 500, nfolds = 3) mae # rank selection for the likelihood-based estimator ecv <- nmf.ffb.cv(Y1, Y2, rank = 1:2, method = "fiml", nfolds = 3) ecv$objfunc
Two quantities that describe what the selected feedback is, as opposed
to whether it is there (which is the job of the calibrated likelihood-ratio
test in nmf.ffb.inference).
Cycle structure. The equilibrium operator is on
the indicators and on the factors, with the same
spectral radius . The strongly connected components of the
directed graph of are the cycles. with a non-empty
support means is nilpotent: the selected feedback is a cascade of
directed cross-factor paths without a closed loop. Since ,
Perron–Frobenius gives a non-negative leading eigenvector, and its
normalized entries say how much of the loop each factor carries.
Factor-level alignment. Feedback of the form
enters each factor through the factor
scores of the factors and is observationally equivalent, under the Gaussian
working model, to a change of and of the factor covariance
; only the departure from that family is identified, through the
uniquenesses . On the free set of the exclusion restriction the
family is a subspace of dimension at most , and
omega is the share of the squared norm of that lies in
it. omega has no natural zero: for an isotropic direction its
expectation is omega0 = dim / n.free, so it should be read against
omega0 and against a simulation in which the true feedback is
indicator-level (there it averages about 1.3 * omega0).
nmf.ffb.diagnostics(object, which = c("selected", "full"))nmf.ffb.diagnostics(object, which = c("selected", "full"))
object |
A fitted |
which |
Which feedback matrix to describe: |
A list with cycles (A.factor, rho,
cycle, n.cycles, perron) and omega
(omega, omega0, dim, n.free).
set.seed(1) Y2 <- matrix(runif(2 * 200), 2, 200) X <- matrix(0.02, 6, 2); X[1:3, 1] <- 1; X[4:6, 2] <- 1 X <- sweep(X, 2, colSums(X), "/") T1 <- matrix(0, 2, 6); T1[1, 5] <- 0.6; T1[2, 2] <- 0.5 L <- solve(diag(6) - X %*% T1) Y1 <- pmax(L %*% (X %*% (matrix(c(1, 0.5, 0.3, 1), 2, 2) %*% Y2) + matrix(rnorm(6 * 200, 0, 0.05), 6, 200)), 0) fit <- nmf.ffb(Y1, Y2, Q = 2) d <- nmf.ffb.diagnostics(fit) d$cycles$rho # spectral radius d$cycles$cycle # which factors lie on a cycle c(d$omega$omega, d$omega$omega0)set.seed(1) Y2 <- matrix(runif(2 * 200), 2, 200) X <- matrix(0.02, 6, 2); X[1:3, 1] <- 1; X[4:6, 2] <- 1 X <- sweep(X, 2, colSums(X), "/") T1 <- matrix(0, 2, 6); T1[1, 5] <- 0.6; T1[2, 2] <- 0.5 L <- solve(diag(6) - X %*% T1) Y1 <- pmax(L %*% (X %*% (matrix(c(1, 0.5, 0.3, 1), 2, 2) %*% Y2) + matrix(rnorm(6 * 200, 0, 0.05), 6, 200)), 0) fit <- nmf.ffb(Y1, Y2, Q = 2) d <- nmf.ffb.diagnostics(fit) d$cycles$rho # spectral radius d$cycles$cycle # which factors lie on a cycle c(d$omega$omega, d$omega$omega0)
Creates a Graphviz DOT script that visualizes the structural network
estimated by nmf.ffb.
The resulting diagram displays:
endogenous observed variables (),
exogenous observed variables (),
latent factors (, ..., ),
together with the non-negative path coefficients whose magnitudes exceed a user-specified threshold.
Directed edges represent estimated relationships:
: entries of C2 (exogenous loadings),
: rows of X (factor-to-endogenous mappings),
: entries of C1 (feedback paths).
Edge widths are scaled by coefficient magnitude, and nodes are placed in optional visual clusters. Only variables participating in edges above the threshold are displayed, while latent factors are always shown.
nmf.ffb.DOT( result, weight_scale = 5, weight_scale_c2 = weight_scale, weight_scale_x1 = weight_scale, weight_scale_feedback = weight_scale, threshold = 0.01, sig.level = 0.1, rankdir = "LR", fill = TRUE, cluster.box = c("normal", "faint", "invisible", "none"), cluster.labels = NULL, hide.isolated = TRUE, model = c("selected", "null", "full"), ... )nmf.ffb.DOT( result, weight_scale = 5, weight_scale_c2 = weight_scale, weight_scale_x1 = weight_scale, weight_scale_feedback = weight_scale, threshold = 0.01, sig.level = 0.1, rankdir = "LR", fill = TRUE, cluster.box = c("normal", "faint", "invisible", "none"), cluster.labels = NULL, hide.isolated = TRUE, model = c("selected", "null", "full"), ... )
result |
A list returned by |
weight_scale |
Base scaling factor for edge widths. |
weight_scale_c2 |
Scaling factor for edges
|
weight_scale_x1 |
Scaling factor for edges
|
weight_scale_feedback |
Scaling factor for feedback edges
|
threshold |
Minimum coefficient value needed for an edge to be drawn. |
sig.level |
Significance level for filtering structural edges
( |
rankdir |
Graphviz rank direction (e.g., |
fill |
Logical; whether to use filled node shapes. |
cluster.box |
Character string controlling the visibility and style
of cluster frames around Y2, factors, and Y1 blocks.
One of |
cluster.labels |
Optional character vector of length 3 giving custom labels for the Y2, factor, and Y1 clusters. |
hide.isolated |
Logical. If |
model |
Which fit of a |
... |
For backward compatibility: accepts deprecated names
|
A character string representing a valid Graphviz DOT script.
nmf.ffb, nmf.ffb.inference,
plot.nmfkc.DOT
Y <- t(iris[, -5]) Y1 <- Y[1:2, ] Y2 <- Y[3:4, ] result <- nmf.ffb(Y1, Y2, rank = 2, maxit = 500) dot <- nmf.ffb.DOT(result) cat(dot)Y <- t(iris[, -5]) Y1 <- Y[1:2, ] Y2 <- Y[3:4, ] result <- nmf.ffb(Y1, Y2, rank = 2, maxit = 500) dot <- nmf.ffb.DOT(result) cat(dot)
Step 1 of the NMF-FFB workflow: select by element-wise
cross-validation of the feed-forward fit. Held-out entries of
are predicted from the remaining ones, so the criterion measures how
well a -factor non-negative basis describes the outcomes, which is what
is for.
nmf.ffb.ecv( Y1, Y2, rank = NULL, X.init = "nndsvd", X.L2.ortho = 100, epsilon = 1e-06, maxit = 5000, ... )nmf.ffb.ecv( Y1, Y2, rank = NULL, X.init = "nndsvd", X.L2.ortho = 100, epsilon = 1e-06, maxit = 5000, ... )
Y1, Y2
|
Endogenous and exogenous matrices (variables in rows, units in columns). |
rank |
Candidate values of |
X.init |
Stage-1 initialization, passed to |
X.L2.ortho |
Orthogonality penalty on the basis in stage 1. Keep the
value that will be used in |
epsilon, maxit
|
Convergence tolerance and iteration cap for stage 1. |
... |
Passed to |
Feedback cannot be cross-validated, and that is not a limitation of the
implementation. By the reduced-form equivalence, the equilibrium mapping
is reproduced exactly by a
feed-forward model with the exogenous matrix : the two models
predict from identically, so no prediction criterion –
column-wise CV, element-wise CV, or a test-sample error – can prefer one over
the other. Feedback is identified by the conditional covariance, not by
prediction, which is why it is addressed by nmf.ffb.test and not
here.
The object returned by nmfkc.ecv: the held-out error by
rank, with the selected rank.
ecv <- nmf.ffb.ecv(Y1, Y2, rank = 1:5) # 1. choose Q <- this function fit <- nmf.ffb(Y1, Y2, rank = ecv$rank) # 2. estimate; BIC selects the support tst <- nmf.ffb.test(fit, Y1, Y2) # 3. test the feed-forward null dgn <- nmf.ffb.diagnostics(fit) # 4. cycles, spectral radius, identifiability inf <- nmf.ffb.inference(fit, Y1, Y2) # 5. intervals for the retained entries
nmf.ffb, nmf.ffb.test,
nmfkc.ecv
set.seed(1) Y <- t(as.matrix(iris[, 1:4])) Y1 <- Y[1:2, ]; Y2 <- Y[3:4, ] ecv <- nmf.ffb.ecv(Y1, Y2, rank = 1:2, nfolds = 3) ecv$rankset.seed(1) Y <- t(as.matrix(iris[, 1:4])) Y1 <- Y[1:2, ]; Y2 <- Y[3:4, ] ecv <- nmf.ffb.ecv(Y1, Y2, rank = 1:2, nfolds = 3) ecv$rank
nmf.ffb.inference performs statistical inference on the structural
coefficient matrices (latent feedback, ) and
(exogenous loading, ) from a fitted
nmf.ffb model, after the feed-forward null has been
rejected by nmf.ffb.test. The
basis is held fixed throughout, which avoids label switching
and gives a clean conditional interpretation: uncertainty of the structural
coefficients given the measurement model. Which resampling scheme is run
depends on object$method.
method = "fiml": a parametric bootstrap from the selected
model. data sets are drawn from the selected model (with its
) and are re-estimated
with and the selected support held fixed. Confidence intervals are
the centred (basic) percentile intervals ; support_rate is
the share of replicates with threshold (0 for
entries outside the selected support), and prob.unsupported is
.
These are post-selection quantities, conditional on the support that
BIC picked, and they can overstate the evidence. They answer "how large is
this entry, given that the procedure kept it", not "is there feedback at
all"; the calibrated statement about the presence of feedback is
nmf.ffb.test, which is the step that belongs before this one.
Each replicate seeds itself (seed + B + b), so the result does not
depend on cores.
method = "mu" (legacy fits, and objects without a
method field). The procedure is a full pair bootstrap that holds the basis
matrix from the original fit fixed across all replicates
(which avoids label switching and gives a clean conditional
interpretation: “uncertainty of the structural coefficients given the
measurement model”):
For each replicate , resample column indices
with replacement from
and form , .
Re-estimate by running the
nmf.ffb multiplicative updates with
held fixed (no update; no centroid sort), using the same
C1.L1, C2.L1 as the original fit.
Discard replicates that violate stationarity
() or have an amplification ratio
exceeding the geometric-series bound by more than 1\
Because are non-negative by construction, exact
zeros are essentially never observed in the bootstrap distribution.
Significance is assessed via a support rate at a small display
threshold (default 0.01):
This is a one-sided counterpart of the classical -value:
large support_rate indicates strong evidence that the entry is
meaningfully positive. Significance markers follow the lavaan
convention with the natural correspondence :
* (sup > 0.95), ** (sup > 0.99), ***
(sup > 0.999). Cutoffs use strict greater-than so the rule
mirrors the standard R convention for p-values (p < 0.05 / 0.01 /
0.001 → //), translated to support_rate via
.
nmf.ffb.inference( object, Y1, Y2, B = 1000L, threshold = 0.01, boot.level = 0.95, C1.L1 = 1, C2.L1 = 0.1, seed = 123L, ... )nmf.ffb.inference( object, Y1, Y2, B = 1000L, threshold = 0.01, boot.level = 0.95, C1.L1 = 1, C2.L1 = 0.1, seed = 123L, ... )
object |
A fitted object returned by |
Y1 |
Endogenous variable matrix (P1 x N). Must match the data
used in |
Y2 |
Exogenous variable matrix (P2 x N). Same. |
B |
Number of bootstrap replicates. Default |
threshold |
Display threshold |
boot.level |
Confidence level for the bootstrap CI.
Default |
C1.L1, C2.L1
|
L1 sparsity penalties used by the original
|
seed |
Base RNG seed for the bootstrap. For |
... |
Hidden options. Shared:
|
The input object with additional bootstrap inference
components:
coefficients |
Data frame with rows for every entry of |
C1.support.rate, C2.support.rate
|
Per-element support rates (Q x P1 and Q x P2 matrices). |
C1.ci.lower, C1.ci.upper, C2.ci.lower, C2.ci.upper
|
Per-element
CI bounds (percentile for |
C1.boot.draws, C2.boot.draws
|
Bootstrap distributions: 3D arrays of shape
B x Q x P1 (and B x Q x P2). Invalid replicates contain |
rho.boot.draws |
Per-replicate spectral radius |
AR.boot, iter.boot
|
( |
#'
boot.B, boot.threshold, boot.level
|
Inputs recorded for reproducibility. |
boot.n.valid, boot.n.invalid
|
Validity counts (for |
boot.method |
How the replicates were drawn. |
identified.approx, identification
|
Whether |
The calibrated test of the feed-forward null – LR.p.boot,
prob.select.null, C1.restriction.change.rate and the null
bootstrap behind them – is not returned here. It moved to
nmf.ffb.test in 0.9.8, so that an object carrying coefficient
intervals cannot be mistaken for a test of feedback (see NEWS).
This function's interface changed at v0.6.8: the legacy 1-step Newton
wild bootstrap (with sandwich SE) has been replaced by the full pair
bootstrap described above, following the paper revision. The fields
sigma2.used, C2.se, C2.se.boot, C2.p.side
that the previous implementation produced are no longer present.
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
Y <- t(iris[, -5]) Y1 <- Y[1:2, ]; Y2 <- Y[3:4, ] res <- nmf.ffb(Y1, Y2, rank = 2) ## intervals for the entries the test has already licensed reading. inf <- nmf.ffb.inference(res, Y1, Y2, B = 20) # use B = 1000 in practice head(inf$coefficients) # estimate, interval and support rate per entry res.mu <- nmf.ffb(Y1, Y2, rank = 2, method = "mu") res.mu2 <- nmf.ffb.inference(res.mu, Y1, Y2, B = 200) head(res.mu2$coefficients)Y <- t(iris[, -5]) Y1 <- Y[1:2, ]; Y2 <- Y[3:4, ] res <- nmf.ffb(Y1, Y2, rank = 2) ## intervals for the entries the test has already licensed reading. inf <- nmf.ffb.inference(res, Y1, Y2, B = 20) # use B = 1000 in practice head(inf$coefficients) # estimate, interval and support rate per entry res.mu <- nmf.ffb(Y1, Y2, rank = 2, method = "mu") res.mu2 <- nmf.ffb.inference(res.mu, Y1, Y2, B = 200) head(res.mu2$coefficients)
Infers a heuristic partition of observed variables into exogenous ()
and endogenous () blocks for use in NMF-FFB.
The method is based on positive-SEM logic, causal ordering, and optional
sign alignment using the first principal component (PC1).
The procedure:
internally standardizes variables (mean 0, sd 1),
optionally flips signs so that most variables align positively with PC1,
infers a causal ordering by repeatedly regressing each variable on the remaining ones and selecting the variable with the largest minimum standardized coefficient,
determines an exogenous block by scanning the ordering from upstream
and stopping at the first variable whose strongest parent coefficient
exceeds threshold.
If n.exogenous is supplied, it overrides the automatic threshold rule.
nmf.ffb.split( x, n.exogenous = NULL, threshold = 0.1, auto.flipped = TRUE, verbose = FALSE, ... )nmf.ffb.split( x, n.exogenous = NULL, threshold = 0.1, auto.flipped = TRUE, verbose = FALSE, ... )
x |
A numeric matrix or data frame with rows = samples and columns = observed variables. |
n.exogenous |
Optional integer specifying the number of exogenous variables
( |
threshold |
Standardized regression-coefficient threshold used in the
automatic exogenous–endogenous split. A variable is treated as endogenous
once its maximum standardized parent coefficient exceeds this value.
(Default: |
auto.flipped |
Logical; if |
verbose |
Logical; if |
... |
Reserved for future use; currently unused (also accepted
by the |
A list with:
endogenous.variables |
Character vector of variables selected as endogenous ( |
exogenous.variables |
Character vector of variables selected as exogenous ( |
ordered.variables |
Variables in inferred causal order (from exogenous to endogenous). |
is.flipped |
Logical vector indicating which variables were sign-flipped during processing. |
n.exogenous |
Integer giving the number of exogenous variables. |
# Infer exogenous/endogenous split from iris sp <- nmf.ffb.split(iris[, -5], n.exogenous = 2) sp$endogenous.variables sp$exogenous.variables# Infer exogenous/endogenous split from iris sp <- nmf.ffb.split(iris[, -5], n.exogenous = 2) sp$endogenous.variables sp$exogenous.variables
Calibrated test of the hypothesis that the outcomes need no indicator-level
feedback: the feed-forward model
with a full random-effect covariance against the feedback model with
. This is the only inferential step of the procedure
(see the workflow below); nmf.ffb.inference is for the
uncertainty of the coefficients after the null has been rejected, and
nmf.ffb.diagnostics for what the selected feedback looks like.
nmf.ffb.test( object, Y1, Y2, B = 1000L, calibration = c("procedure", "conditional"), seed = 123L, ... )nmf.ffb.test( object, Y1, Y2, B = 1000L, calibration = c("procedure", "conditional"), seed = 123L, ... )
object |
A fit from |
Y1, Y2
|
The endogenous and exogenous matrices the fit was computed from (variables in rows, units in columns). |
B |
Number of bootstrap replicates (default 1000). The smallest
reportable |
calibration |
What is re-applied to each null replicate; see above.
|
seed |
Base seed; replicate |
... |
Passed to the bootstrap: |
The likelihood ratio does not have a null distribution: the
constraint puts the null on the boundary of the
parameter space, and the support of is chosen from the same
data. It is therefore calibrated by a parametric bootstrap from the fitted
feed-forward null,
with negative draws clipped at zero, and the reported -value is
.
calibration = "procedure" (the default) re-applies the whole
procedure: stage 1 is re-run on each , the basis is
re-estimated, the exclusion restriction is re-derived from that basis, and
stage 2 follows. The -value is then the operating characteristic of
the complete exploratory procedure, which is what is applied to a fresh data
set. This is the calibration to use.
"conditional" holds and the exclusion restriction at
their fitted values and re-runs stage 2 only. It is valid only if basis and
restriction are independent of the being tested, which they are not when
all three come from the same sample; it is provided for comparison with that
practice, and it is anti-conservative.
Sample splitting – basis and restriction from one half of the units, test on
other – also removes the dependence, and was offered here in 0.9.8. It was
withdrawn in 0.9.8: measured against "procedure" it has the same size and
lower power, because the test uses units, and it cannot be run at
all when the halves are too small for stage 1.
Report LR.p.boot together with prob.select.null and, for
calibration = "procedure", C1.restriction.change.rate. The last two say
whether the exclusion restriction is determined well enough for the test to
mean anything: if the dominant factor of an outcome moves from one null
replicate to the next, a previously blocked entry becomes free and a feedback
coefficient can absorb loading structure that the feed-forward model
attributes to . A large prob.select.null with a large null
quantile is the signature of a data set in which the procedure finds feedback
whether or not there is any.
An object of class "nmf.ffb.test":
Observed c(full, selected) likelihood ratios and the
number of free feedback entries left by the exclusion restriction.
Calibrated -values of the two statistics.
The 0.95 quantiles of the null distributions.
matrix of the null replicates.
Null false-selection rate: the probability that BIC retains at least one feedback entry when the feed-forward model is true.
("procedure" only) share of null replicates in which
the exclusion restriction moved.
Free entries per null replicate (varies under "procedure").
Usable replicates, and how many did not meet the optimizer tolerance.
ecv <- nmf.ffb.ecv(Y1, Y2, rank = 1:5) # 1. choose Q fit <- nmf.ffb(Y1, Y2, rank = ecv$rank) # 2. estimate; BIC selects the support tst <- nmf.ffb.test(fit, Y1, Y2) # 3. test the feed-forward null <- only inference dgn <- nmf.ffb.diagnostics(fit) # 4. cycles, spectral radius, identifiability inf <- nmf.ffb.inference(fit, Y1, Y2) # 5. intervals for the retained entries
nmf.ffb, nmf.ffb.ecv,
nmf.ffb.inference, nmf.ffb.diagnostics
This function is experimental and still under development. The interface may change in future versions: argument names, defaults and the contents of the returned object are not yet stable. Code written against it today may need adjusting after an update. The rest of the package does not carry this caveat.
nmf.gmm fits the model
by a generalized EM algorithm: a -component Gaussian mixture is placed
on the latent NMF scores , whose class-conditional mean combines
a covariate regression ( in the paper, the
coefficient matrix) with a class shift . The
shared basis stays non-negative and column-normalized. Clustering is
model based, through the posterior responsibilities. (tied
covariance) reduces exactly to NMF with random effects (nmfre).
Like the rest of the package, nmf.gmm performs optimization
only; use nmf.gmm.inference for standard errors and tests of
, and nmf.gmm.select to choose .
nmf.gmm(Y, A = NULL, rank, K = 1, ...)nmf.gmm(Y, A = NULL, rank, K = 1, ...)
Y |
Data matrix |
A |
Covariate matrix |
rank |
Integer rank |
K |
Integer number of mixture components. Default 1 (= NMF-RE). |
... |
Additional arguments:
|
An object of class "nmf.gmm": a list with X
(basis), X0 (the column-normalized initial basis the EM started
from), C (, Q x R), mu (class means, Q x K),
tau2, sigma2, xi (mixing proportions), gamma
(responsibilities, N x K), cluster (hard labels), loglik,
BIC, ICL, n.params, entropy, Yhat,
objfunc(.iter), iter, converged, stop_by,
K,
rank, cov, dims and runtime. With
class.effects, also class.effects (the freed rows of
A, by name), class.rows (their indices) and C.class
(a Q x |S| x K array of class-specific coefficients; the class-specific
effect in observation space is X %*% C.class[, , k]).
nmf.gmm.inference, nmf.gmm.select,
nmfre (the special case).
set.seed(1) Y <- matrix(abs(rnorm(20 * 60)) + 1, 20, 60) A <- rbind(1, rnorm(60)) # intercept + one covariate fit <- nmf.gmm(Y, A, rank = 2, K = 2) fit$BIC table(fit$cluster)set.seed(1) Y <- matrix(abs(rnorm(20 * 60)) + 1, 20, 60) A <- rbind(1, rnorm(60)) # intercept + one covariate fit <- nmf.gmm(Y, A, rank = 2, K = 2) fit$BIC table(fit$cluster)
This function is experimental and still under development. The interface may change in future versions: argument names, defaults and the contents of the returned object are not yet stable. Code written against it today may need adjusting after an update. The rest of the package does not carry this caveat.
nmf.gmm.inference adds Wald inference on the covariate-coefficient
matrix () of a fitted nmf.gmm object,
conditional on the estimated basis and mixture. It reports the
analytic outer-product (mixture-information) standard error — the
coverage-calibrated choice at — together with a
wild-bootstrap standard error and percentile confidence interval,
mirroring nmfre.inference / nmfkc.inference.
Supported for the tied-covariance variant.
nmf.gmm.inference(object, Y, A = object$A, ...)nmf.gmm.inference(object, Y, A = object$A, ...)
object |
A fitted |
Y |
The data matrix used in the fit. |
A |
The covariate matrix used in the fit. Default |
... |
Additional arguments, named as in the other wild-bootstrap
inference functions of the package: |
object with class c("nmf.gmm.inference", "nmf.gmm")
and added components: coefficients (data frame with Basis,
Covariate, Estimate, SE, BSE, z_value,
p_value, CI_low, CI_high), C.se.outer,
C.se.sandwich, C.se.boot (Q x R matrices) and Xi
().
set.seed(1) Y <- matrix(abs(rnorm(20 * 60)) + 1, 20, 60); A <- rbind(1, rnorm(60)) fit <- nmf.gmm(Y, A, rank = 2, K = 2) fit <- nmf.gmm.inference(fit, Y, A, wild.B = 300) head(fit$coefficients)set.seed(1) Y <- matrix(abs(rnorm(20 * 60)) + 1, 20, 60); A <- rbind(1, rnorm(60)) fit <- nmf.gmm(Y, A, rank = 2, K = 2) fit <- nmf.gmm.inference(fit, Y, A, wild.B = 300) head(fit$coefficients)
This function is experimental and still under development. The interface may change in future versions: argument names, defaults and the contents of the returned object are not yet stable. Code written against it today may need adjusting after an update. The rest of the package does not carry this caveat.
nmf.gmm.select fits nmf.gmm over a vector of
values and reports the log-likelihood, BIC and ICL for each, selecting
K.best by BIC (K.best.icl by ICL). If a vector of known labels
is supplied via truth, the adjusted Rand index of the hard clustering
is added.
nmf.gmm.select(Y, A = NULL, rank, K = 1:5, ...)nmf.gmm.select(Y, A = NULL, rank, K = 1:5, ...)
Y, A, rank
|
As in |
K |
Integer vector of candidate component counts. Default |
... |
Additional arguments passed to |
An object of class "nmf.gmm.select": a list with the
table (data frame over K), K.best (min BIC),
K.best.icl (min ICL) and fits (the fitted objects).
set.seed(1) Y <- matrix(abs(rnorm(20 * 80)) + 1, 20, 80); A <- rbind(1, rnorm(80)) sel <- nmf.gmm.select(Y, A, rank = 2, K = 1:4, verbose = FALSE) sel$K.bestset.seed(1) Y <- matrix(abs(rnorm(20 * 80)) + 1, 20, 80); A <- rbind(1, rnorm(80)) sel <- nmf.gmm.select(Y, A, rank = 2, K = 1:4, verbose = FALSE) sel$K.best
This function is experimental and still under development. The interface may change in future versions: argument names, defaults and the contents of the returned object are not yet stable. Code written against it today may need adjusting after an update. The rest of the package does not carry this caveat.
nmf.gmm.twostage runs the two-stage route that
nmf.gmm is designed to improve on, as a matched baseline:
(1) estimate least-squares scores on the initial basis, (2) regress the
scores on the covariates blind to the class and keep the residuals,
(3) reconstitute the residuals in observation space, shift them to
non-negativity, and (4) refit an intercept-only nmf.gmm from the
same basis initialization. Both routes start from the same basis and
use the same mixture family and settings. They differ in the order of
adjustment and clustering and, unless fixX = TRUE, in the basis each
refits to its own data: the two-stage route refits it to the
reconstituted, shifted residuals. When the covariate is associated with the
class, the class-blind regression of stage (2) displaces the class means
(the oracle proposition of Satoh 2026); when the covariate is
(near-)mean-independent of the class, that displacement vanishes in the
population, so both routes target the same class means, although the
fitted partitions need not coincide in a finite sample.
nmf.gmm.twostage(Y, A = NULL, rank, K = 1, ...)nmf.gmm.twostage(Y, A = NULL, rank, K = 1, ...)
Y |
Data matrix |
A |
Covariate matrix |
rank |
Integer rank |
K |
Integer number of mixture components. |
... |
Additional arguments as in |
An object of class c("nmf.gmm.twostage", "nmf.gmm"): the
stage-2 fit (all nmf.gmm fields and S3 methods apply),
plus a twostage list with the non-negativity shift, the
covariate matrix A that was removed, and the shared basis
initialization X0.
nmf.gmm (the joint route this baselines).
set.seed(1) Y <- matrix(abs(rnorm(12 * 40)) + 1, 12, 40) A <- rbind(1, rnorm(40)) ts <- nmf.gmm.twostage(Y, A, rank = 2, K = 2, nstart = 2, maxit = 100) table(ts$cluster)set.seed(1) Y <- matrix(abs(rnorm(12 * 40)) + 1, 12, 40) A <- rbind(1, rnorm(40)) ts <- nmf.gmm.twostage(Y, A, rank = 2, K = 2, nstart = 2, maxit = 100) table(ts$cluster)
nmfae fits a three-layer nonnegative matrix factorization model
, where is a decoder basis
(column sum 1), is a bottleneck parameter matrix,
is an encoder basis (row sum 1), and is the input matrix.
When Y2 = Y1, the model acts as a non-negative autoencoder.
When Y1 != Y2, it acts as a heteroencoder.
Initialization uses a three-step NMF procedure via nmfkc:
(1) nmfkc(Y1, rank=Q) to obtain ,
(2) nmfkc(Y1, A=Y2, rank=Q) with fixed to obtain ,
(3) nmfkc(Y2, rank=R) to factor into and .
nmf.rrr( Y1, Y2 = Y1, rank1 = 2, rank2 = NULL, epsilon = 1e-04, maxit = 5000, verbose = FALSE, ... )nmf.rrr( Y1, Y2 = Y1, rank1 = 2, rank2 = NULL, epsilon = 1e-04, maxit = 5000, verbose = FALSE, ... )
Y1 |
Output matrix |
Y2 |
Input matrix |
rank1 |
Integer. Rank of the response basis |
rank2 |
Integer. Rank of the covariate basis |
epsilon |
Positive convergence tolerance. Default is |
maxit |
Maximum number of multiplicative update iterations. Default is 5000. |
verbose |
Logical. If |
... |
Additional arguments:
Rank aliases accepted here for backward compatibility:
|
An object of class "nmfae", a list with components:
X1 |
Decoder basis matrix (P1 x Q), column sum 1. |
C |
Parameter matrix (Q x R). |
X2 |
Encoder basis matrix (R x P2), row sum 1. |
Y1hat |
Fitted values |
B1 |
Decoder-side scores |
H |
Deprecated; identical to |
B2 |
Encoder-side scores |
B1.prob, B2.prob
|
Column-normalized |
B1.cluster, B2.cluster
|
Hard label per sample (argmax over the columns
of the corresponding |
X1.prob |
Row-normalized |
X1.cluster |
Hard response-group label per response variable
(argmax over |
X2.prob |
Column-normalized |
X2.cluster |
Hard covariate-group label per covariate variable
(argmax over |
rank |
Named integer vector |
method |
Objective used ( |
objfunc |
Final objective value. |
objfunc.iter |
Objective values by iteration. |
r.squared |
|
r.squared.uncentered |
Uncentered |
r.squared.centered |
Row-mean centered |
niter |
Number of iterations performed. |
runtime |
Elapsed time as a |
n.missing |
Number of missing (or zero-weighted) elements in |
n.total |
Total number of elements in |
This function is experimental. The interface may change in future versions.
Satoh, K. (2025). Applying Non-negative Matrix Factorization with Covariates to Multivariate Time Series. Japanese Journal of Statistics and Data Science.
Satoh, K. and Tokuda, Y. (2026). Co-clustering of Response and Covariate Variables by Tri-Factorizing Their Non-negative Regression Coefficient Matrix. arXiv preprint arXiv:2607.27474. doi:10.48550/arXiv.2607.27474
Lee, D. D. and Seung, H. S. (2001). Algorithms for Non-negative Matrix Factorization. Advances in Neural Information Processing Systems, 13.
Saha, S. et al. (2022). Hierarchical Deep Learning Neural Network (HiDeNN): An Artificial Intelligence (AI) Framework for Computational Science and Engineering. Computer Methods in Applied Mechanics and Engineering, 399.
nmf.rrr.inference, predict.nmfae, nmf.rrr.ecv, nmf.rrr.DOT, nmfkc
# Autoencoder example Y <- matrix(c(1,0,1,0, 0,1,0,1, 1,1,0,0), nrow=3, byrow=TRUE) res <- nmf.rrr(Y, rank1=2, rank2=2) res$r.squared # Heteroencoder example Y1 <- matrix(c(1,0,0,1), nrow=2) Y2 <- matrix(runif(8), nrow=4) res2 <- nmf.rrr(Y1, Y2, rank1=2, rank2=2)# Autoencoder example Y <- matrix(c(1,0,1,0, 0,1,0,1, 1,1,0,0), nrow=3, byrow=TRUE) res <- nmf.rrr(Y, rank1=2, rank2=2) res$r.squared # Heteroencoder example Y1 <- matrix(c(1,0,0,1), nrow=2) Y2 <- matrix(runif(8), nrow=4) res2 <- nmf.rrr(Y1, Y2, rank1=2, rank2=2)
nmfae.cv performs k-fold cross-validation by splitting columns (samples)
of and into div folds. For each fold, the model
is fitted on the training samples and
predictive performance is evaluated on the held-out samples.
When Y2 is a kernel matrix created by nmfkc.kernel
(detected via attributes), the symmetric kernel splitting convention is used:
Y2[train, train] for training and Y2[train, test] for prediction.
nmf.rrr.cv(Y1, Y2 = Y1, rank1 = 2, rank2 = NULL, ...)nmf.rrr.cv(Y1, Y2 = Y1, rank1 = 2, rank2 = NULL, ...)
Y1 |
Output matrix |
Y2 |
Input matrix |
rank1 |
Integer. Rank of the response basis. Default is 2. |
rank2 |
Integer. Rank of the covariate basis. Default ( |
... |
Additional arguments passed to Rank aliases accepted here for backward compatibility:
|
A list with components:
objfunc |
Mean squared error per valid element over all folds. |
sigma |
Residual standard error (RMSE), same scale as |
objfunc.block |
Per-fold squared error totals. |
block |
Integer vector of fold assignments (1, ..., |
This function is experimental. The interface may change in future versions; details are to be described in an upcoming paper.
nmf.rrr, nmf.rrr.ecv, nmf.rrr.kernel.beta.cv,
nmfkc.cv
Y <- t(iris[1:30, 1:4]) res <- nmf.rrr.cv(Y, rank1 = 2, rank2 = 2, nfolds = 5, maxit = 500) res$sigmaY <- t(iris[1:30, 1:4]) res <- nmf.rrr.cv(Y, rank1 = 2, rank2 = 2, nfolds = 5, maxit = 500) res$sigma
nmfae.DOT generates a DOT language string for visualizing the structure
of a three-layer NMF model. Two graph types are supported:
"XCX" shows encoder factors, , and decoder factors;
"YXCXY" shows the full structure from through , ,
to .
Edge widths are proportional to matrix element values, and edges below
threshold are omitted for clarity.
nmf.rrr.DOT( result, type = c("XCX", "YXCXY"), threshold = 0.01, sig.level = 0.1, rankdir = "LR", fill = TRUE, weight_scale = 5, weight_scale_x1 = weight_scale, weight_scale_theta = weight_scale, weight_scale_x2 = weight_scale, Y1.label = NULL, X1.label = NULL, X2.label = NULL, Y2.label = NULL, Y1.title = "Output (Y1)", X1.title = "Response (X1)", X2.title = "Covariate (X2)", Y2.title = "Input (Y2)", hide.isolated = TRUE )nmf.rrr.DOT( result, type = c("XCX", "YXCXY"), threshold = 0.01, sig.level = 0.1, rankdir = "LR", fill = TRUE, weight_scale = 5, weight_scale_x1 = weight_scale, weight_scale_theta = weight_scale, weight_scale_x2 = weight_scale, Y1.label = NULL, X1.label = NULL, X2.label = NULL, Y2.label = NULL, Y1.title = "Output (Y1)", X1.title = "Response (X1)", X2.title = "Covariate (X2)", Y2.title = "Input (Y2)", hide.isolated = TRUE )
result |
An object of class |
type |
Character. Graph type: |
threshold |
Numeric. Edges with values below this are omitted. Default is 0.01. |
sig.level |
Numeric or |
rankdir |
Character. Graph direction for DOT layout. Default is |
fill |
Logical. If |
weight_scale |
Numeric. Base scale factor for edge widths. Default is 5. |
weight_scale_x1 |
Numeric. Scale factor for |
weight_scale_theta |
Numeric. Scale factor for |
weight_scale_x2 |
Numeric. Scale factor for |
Y1.label |
Character vector of output variable labels. |
X1.label |
Character vector of decoder basis labels. |
X2.label |
Character vector of encoder basis labels. |
Y2.label |
Character vector of input variable labels. |
Y1.title |
Character. Title for output node group. Default is |
X1.title |
Character. Title for the response-basis node group. Default is |
X2.title |
Character. Title for the covariate-basis node group. Default is |
Y2.title |
Character. Title for input node group. Default is |
hide.isolated |
Logical. If |
A character string containing the DOT graph specification.
This function is experimental. The interface may change in future versions; details are to be described in an upcoming paper.
set.seed(1) Y <- matrix(runif(20), nrow = 4) res <- nmf.rrr(Y, rank1 = 2) dot <- nmf.rrr.DOT(res)set.seed(1) Y <- matrix(runif(20), nrow = 4) res <- nmf.rrr(Y, rank1 = 2) dot <- nmf.rrr.DOT(res)
nmfae.ecv performs k-fold element-wise cross-validation by randomly
holding out individual elements of , assigning them a weight of 0
via Y1.weights, and evaluating the reconstruction error on those
held-out elements.
This method (also known as Wold's CV) is suitable for determining the optimal
rank pair in three-layer NMF. Both rank1 and rank2 accept
vector inputs. When rank2 = NULL (default), rank2 is set equal to rank1
and pairs are evaluated element-wise (i.e., ).
When rank.encoder is explicitly specified, all combinations of rank and rank.encoder
are evaluated via expand.grid.
nmf.rrr.ecv(Y1, Y2 = Y1, rank1 = 1:2, rank2 = NULL, ...)nmf.rrr.ecv(Y1, Y2 = Y1, rank1 = 1:2, rank2 = NULL, ...)
Y1 |
Output matrix |
Y2 |
Input matrix |
rank1 |
Integer vector of response-basis ranks to evaluate. Default is |
rank2 |
Integer vector of covariate-basis ranks to evaluate. Default is |
... |
Additional arguments passed to Rank aliases accepted here for backward compatibility:
|
A list with components:
objfunc |
Named numeric vector of mean MSE for each (Q, R) pair. |
sigma |
Named numeric vector of RMSE (square root of MSE) for each pair. |
objfunc.fold |
Named list of per-fold MSE vectors for each pair. |
folds |
List of length |
QR |
Data frame with columns |
This function is experimental. The interface may change in future versions; details are to be described in an upcoming paper.
Y <- t(iris[1:30, 1:4]) # Default: rank2=NULL -> paired rank1=rank2 res <- nmf.rrr.ecv(Y, rank1 = 1:2, nfolds = 3, maxit = 200) res$sigma # Explicit rank2: full grid (kept small so the example stays under CRAN's # 10s budget -- a real sweep would use wider ranges and a larger maxit) res2 <- nmf.rrr.ecv(Y, rank1 = 1:2, rank2 = 1:2, nfolds = 3, maxit = 200) res2$sigmaY <- t(iris[1:30, 1:4]) # Default: rank2=NULL -> paired rank1=rank2 res <- nmf.rrr.ecv(Y, rank1 = 1:2, nfolds = 3, maxit = 200) res$sigma # Explicit rank2: full grid (kept small so the example stays under CRAN's # 10s budget -- a real sweep would use wider ranges and a larger maxit) res2 <- nmf.rrr.ecv(Y, rank1 = 1:2, rank2 = 1:2, nfolds = 3, maxit = 200) res2$sigma
nmfae.heatmap displays the three factor matrices , ,
and as side-by-side heatmaps. This provides an alternative to DOT graph
visualization, especially when has many variables (e.g., kernel matrix).
nmf.rrr.heatmap( x, Y1.label = NULL, X1.label = NULL, X2.label = NULL, Y2.label = NULL, palette = NULL, ... )nmf.rrr.heatmap( x, Y1.label = NULL, X1.label = NULL, X2.label = NULL, Y2.label = NULL, palette = NULL, ... )
x |
An object of class |
Y1.label |
Character vector of output variable names (rows of |
X1.label |
Character vector of decoder basis labels (columns of |
X2.label |
Character vector of encoder basis labels (rows of |
Y2.label |
Character vector of input variable names (columns of |
palette |
Color palette vector. Default is white-orange-red (64 colors). |
... |
Not used. |
Invisible NULL. Called for its side effect (plot).
This function is experimental. The interface may change in future versions; details are to be described in an upcoming paper.
nmf.rrr, plot.nmfae, nmf.rrr.DOT
set.seed(1) Y <- matrix(runif(20), nrow = 4) res <- nmf.rrr(Y, rank1 = 2) nmf.rrr.heatmap(res)set.seed(1) Y <- matrix(runif(20), nrow = 4) res <- nmf.rrr(Y, rank1 = 2) nmf.rrr.heatmap(res)
Performs post-estimation inference for in the three-layer NMF model
, conditional on .
Uses sandwich covariance estimation and one-step wild bootstrap with
non-negative projection.
nmf.rrr.inference(object, Y1, Y2 = Y1, wild.bootstrap = TRUE, ...)nmf.rrr.inference(object, Y1, Y2 = Y1, wild.bootstrap = TRUE, ...)
object |
An object of class |
Y1 |
Output matrix |
Y2 |
Input matrix |
wild.bootstrap |
Logical. If |
... |
Additional arguments:
|
An object of class c("nmfae.inference", "nmfae"), inheriting all
components from the input object, with additional inference components:
sigma2.used |
Estimated |
C.se |
Sandwich standard errors for |
C.se.boot |
Bootstrap standard errors for |
C.ci.lower |
Lower CI bounds for |
C.ci.upper |
Upper CI bounds for |
coefficients |
Data frame with Estimate, SE, BSE, z, p-value for each element of |
C.p.side |
P-value type used. |
identified.approx |
Whether the factorization is unique up to
permutations of both latent spaces, i.e. under
|
identification |
The two checks, |
This function is experimental. The interface may change in future versions; details are to be described in an upcoming paper.
nmf.rrr, summary.nmfae.inference
Y <- matrix(c(1,0,1,0, 0,1,0,1, 1,1,0,0), nrow=3, byrow=TRUE) res <- nmf.rrr(Y, rank1=2, rank2=2) res2 <- nmf.rrr.inference(res, Y) summary(res2)Y <- matrix(c(1,0,1,0, 0,1,0,1, 1,1,0,0), nrow=3, byrow=TRUE) res <- nmf.rrr(Y, rank1=2, rank2=2) res2 <- nmf.rrr.inference(res, Y) summary(res2)
nmfae.kernel.beta.cv selects the optimal beta parameter of the
kernel function by evaluating nmf.rrr.cv for each candidate value.
The kernel matrix replaces in the three-layer
NMF model.
When beta = NULL, candidate values are automatically generated via
nmfkc.kernel.beta.nearest.med.
nmf.rrr.kernel.beta.cv( Y1, rank1 = 2, rank2 = NULL, U, V = NULL, beta = NULL, plot = TRUE, ... )nmf.rrr.kernel.beta.cv( Y1, rank1 = 2, rank2 = NULL, U, V = NULL, beta = NULL, plot = TRUE, ... )
Y1 |
Output matrix |
rank1 |
Integer. Rank of the response basis. Default is 2. |
rank2 |
Integer. Rank of the covariate basis. Default ( |
U |
Covariate matrix |
V |
Covariate matrix |
beta |
Numeric vector of candidate beta values. If |
plot |
Logical. If |
... |
Additional arguments. Kernel-specific args ( Rank aliases accepted here for backward compatibility:
|
A list with components:
beta |
The beta value that minimizes the cross-validation objective. |
objfunc |
Named numeric vector of objective function values for each candidate beta. |
nmf.rrr.cv, nmfkc.kernel,
nmfkc.kernel.beta.cv
Y <- matrix(cars$dist, nrow = 1) U <- matrix(cars$speed, nrow = 1) res <- nmf.rrr.kernel.beta.cv(Y, rank1 = 1, rank2 = 1, U = U, beta = c(0.01, 0.02, 0.05), nfolds = 5) res$betaY <- matrix(cars$dist, nrow = 1) U <- matrix(cars$speed, nrow = 1) res <- nmf.rrr.kernel.beta.cv(Y, rank1 = 1, rank2 = 1, U = U, beta = c(0.01, 0.02, 0.05), nfolds = 5) res$beta
Fits nmf.rrr with a paired decoder/encoder rank
() across a range of ranks and reports r.squared,
the effective rank (of the latent encoding ), and the
element-wise CV error sigma.ecv, with the same concise plot as
nmfkc.rank. For a full grid use
nmf.rrr.ecv with rank.encoder and its heatmap.
nmf.rrr.rank( Y1, Y2 = Y1, rank1 = 1:5, detail = c("full", "fast"), plot = TRUE, ... )nmf.rrr.rank( Y1, Y2 = Y1, rank1 = 1:5, detail = c("full", "fast"), plot = TRUE, ... )
Y1 |
Endogenous matrix ( |
Y2 |
Exogenous matrix; defaults to |
rank1 |
Integer vector of (paired) ranks to evaluate (both bases use
the same value). Legacy |
detail |
|
plot |
Logical; draw the diagnostics plot (default |
... |
Passed on to |
A list with rank.best and criteria
(rank, effective.rank, effective.rank.ratio,
r.squared, sigma.ecv).
Roy, O., & Vetterli, M. (2007). The effective rank: A measure of
effective dimensionality. Proc. EUSIPCO, 606–610.
(effective.rank)
Wold, S. (1978). Cross-validatory estimation of the number of
components in factor and principal components models.
Technometrics, 20(4), 397–405. (sigma.ecv)
nmf.rrr, nmf.rrr.ecv,
nmfkc.rank
Assigns user-specified names to the decoder (X1 columns) and encoder
(X2 rows) bases of an nmfae object. The names propagate to
, the coefficients table, and all downstream displays
such as summary, nmfae.DOT, and nmfae.heatmap.
nmf.rrr.rename(x, X1.colnames = NULL, X2.rownames = NULL)nmf.rrr.rename(x, X1.colnames = NULL, X2.rownames = NULL)
x |
An object of class |
X1.colnames |
Character vector of length |
X2.rownames |
Character vector of length |
A modified copy of x with updated names.
set.seed(1) Y <- matrix(runif(15), nrow = 3) res <- nmf.rrr(Y, rank1 = 2, rank2 = 2) res <- nmf.rrr.rename(res, X1.colnames = c("Resp1", "Resp2"), X2.rownames = c("Cov1", "Cov2")) summary(res)set.seed(1) Y <- matrix(runif(15), nrow = 3) res <- nmf.rrr(Y, rank1 = 2, rank2 = 2) res <- nmf.rrr.rename(res, X1.colnames = c("Resp1", "Resp2"), X2.rownames = c("Cov1", "Cov2")) summary(res)
nmfae.signed fits a three-layer non-negative matrix factorization
autoencoder with a signed bottleneck, solving
where is the signed bottleneck. The basis
matrices (columns sum to 1) and (rows sum to 1) retain
their non-negative "parts-based" interpretability, while can
express anti-correlations (e.g., refractive index up vs. Abbe number down).
The algorithm uses Direct Multiplicative Updates derived from
Ding et al. (2010) sign-splitting technique, applied block-wise to the
four non-negative blocks . Each block update
monotonically decreases the true objective
(Lee-Seung auxiliary
function method).
Relation to nmf.rrr: When suffices
(the nmfae case), nmfae.signed reduces to nmfae up to
the parameterization. Use nmfae.signed when the
data exhibit negative cross-property correlations that tri-NMF-AE cannot
express (e.g., high refractive index <-> low Abbe number trade-off).
nmf.rrr.signed( Y1, Y2 = Y1, rank1 = 2, rank2 = NULL, epsilon = 1e-04, maxit = 5000, verbose = FALSE, ... )nmf.rrr.signed( Y1, Y2 = Y1, rank1 = 2, rank2 = NULL, epsilon = 1e-04, maxit = 5000, verbose = FALSE, ... )
Y1 |
Output matrix |
Y2 |
Input matrix |
rank1 |
Integer. Response-basis rank Q. Default 2. |
rank2 |
Integer. Covariate-basis rank R. Default ( |
epsilon |
Relative convergence tolerance on the objective.
Default |
maxit |
Maximum iterations. Default 5000. |
verbose |
Logical. Print progress. Default |
... |
Additional arguments:
Rank aliases accepted here for backward compatibility:
|
An object of class c("nmfae.signed", "nmfae", "nmf") with:
X1 |
Decoder basis (P1 x Q), column sum 1. |
Cp, Cn
|
Non-negative parts of |
C |
Signed bottleneck |
X2 |
Encoder basis (R x P2), row sum 1. |
Y1hat |
Fitted values |
B1 |
Decoder-side scores |
B2 |
Encoder-side scores |
H |
Deprecated; identical to |
B.prob, B.cluster
|
Sample-level soft/hard clustering from a
non-negative clip of the (possibly signed) |
X1.prob |
Row-normalized |
X1.cluster |
Hard response-group label per response variable. |
X2.prob |
Column-normalized |
X2.cluster |
Hard covariate-group label per covariate variable. |
rank |
|
dims |
|
objfunc, objfunc.iter
|
Final and per-iteration objective values. |
r.squared |
|
r.squared.uncentered |
Uncentered |
r.squared.centered |
Row-mean centered |
sigma, mae
|
Residual SE and mean absolute error. |
niter, runtime
|
Iterations and elapsed seconds. |
Y.signed |
Logical; whether |
call |
Matched call. |
This function is experimental; interface may change.
Satoh, K. and Tokuda, Y. (2026). Co-clustering of Response and Covariate Variables by Tri-Factorizing Their Non-negative Regression Coefficient Matrix. arXiv preprint arXiv:2607.27474. doi:10.48550/arXiv.2607.27474
Ding, C.H.Q., Li, T., and Jordan, M.I. (2010). Convex and Semi-Nonnegative Matrix Factorizations. IEEE TPAMI, 32(1), 45-55.
Satoh, K. (2026). Signed-Bottleneck NMF-AE: Signed-Bottleneck 3-Layer NMF (research memo, 2026-04-18).
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1), 45–55.
nmf.rrr, predict.nmfae.signed,
summary.nmfae.signed, nmf.rrr.signed.rename
set.seed(1) Y1 <- matrix(abs(rnorm(12)), 3, 4) Y2 <- matrix(abs(rnorm(20)), 5, 4) res <- nmf.rrr.signed(Y1, Y2, rank1 = 2, rank2 = 2, maxit = 500) summary(res)set.seed(1) Y1 <- matrix(abs(rnorm(12)), 3, 4) Y2 <- matrix(abs(rnorm(20)), 5, 4) res <- nmf.rrr.signed(Y1, Y2, rank1 = 2, rank2 = 2, maxit = 500) summary(res)
Element-wise k-fold cross-validation for nmf.rrr.signed to
select the decoder / encoder ranks . Mirrors
nmf.rrr.ecv but uses the weighted Signed-Bottleneck NMF-AE fit
path (Y1.weights): test-fold elements are zero-weighted during
fitting, and held-out MSE is computed on those elements.
nmf.rrr.signed.ecv(Y1, Y2 = Y1, rank1 = 1:2, rank2 = NULL, ...)nmf.rrr.signed.ecv(Y1, Y2 = Y1, rank1 = 1:2, rank2 = NULL, ...)
Y1 |
Output matrix (P1 x N). |
Y2 |
Input matrix (P2 x N). Default |
rank1 |
Integer vector of candidate response-basis ranks. Default |
rank2 |
Integer vector of candidate covariate-basis ranks, or |
... |
Additional arguments:
Rank aliases accepted here for backward compatibility:
|
An object of class c("nmfae.signed.ecv", "nmfae.ecv") with
objfunc (MSE per pair), sigma (RMSE), objfunc.fold
(per-fold MSE), folds, QR, paired.
This function is experimental. The interface may change in future versions.
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1), 45–55.
Displays the factor blocks of a nmf.rrr.signed fit as
side-by-side heatmaps. Non-negative blocks () use the white-orange-red palette; the signed combined
bottleneck is rendered with a diverging
blue-white-red palette so positive and negative weights are visually
distinguishable.
nmf.rrr.signed.heatmap( x, Y1.label = NULL, X1.label = NULL, X2.label = NULL, Y2.label = NULL, palette.pos = NULL, palette.signed = NULL, show.C = TRUE, ... )nmf.rrr.signed.heatmap( x, Y1.label = NULL, X1.label = NULL, X2.label = NULL, Y2.label = NULL, palette.pos = NULL, palette.signed = NULL, show.C = TRUE, ... )
x |
An object of class |
Y1.label |
Character vector for rows of |
X1.label |
Decoder basis labels. |
X2.label |
Encoder basis labels. |
Y2.label |
Input variable labels. |
palette.pos |
Palette for non-negative blocks. Default white-orange-red. |
palette.signed |
Palette for signed |
show.C |
Logical. If |
... |
Not used. |
Invisible NULL. Called for its side effect (plot).
This function is experimental. The interface may change in future versions.
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1), 45–55.
nmf.rrr.signed, nmf.rrr.heatmap
set.seed(1) Y1 <- matrix(abs(rnorm(12)), 3, 4) Y2 <- matrix(abs(rnorm(20)), 5, 4) res <- nmf.rrr.signed(Y1, Y2, rank1 = 2, rank2 = 2, maxit = 200) nmf.rrr.signed.heatmap(res)set.seed(1) Y1 <- matrix(abs(rnorm(12)), 3, 4) Y2 <- matrix(abs(rnorm(20)), 5, 4) res <- nmf.rrr.signed(Y1, Y2, rank1 = 2, rank2 = 2, maxit = 200) nmf.rrr.signed.heatmap(res)
Post-estimation inference for the signed bottleneck
in the Signed-Bottleneck NMF-AE model
, conditional on
. Uses sandwich covariance and wild bootstrap
without the non-negativity projection that nmf.rrr.inference
applies (because is unconstrained in sign here).
nmf.rrr.signed.inference(object, Y1, Y2 = Y1, wild.bootstrap = TRUE, ...)nmf.rrr.signed.inference(object, Y1, Y2 = Y1, wild.bootstrap = TRUE, ...)
object |
A fitted |
Y1 |
Output matrix used during fitting. |
Y2 |
Input matrix used during fitting. Default |
wild.bootstrap |
Logical. Default |
... |
Additional arguments:
|
An object of class c("nmfae.signed.inference",
"nmfae.inference", "nmfae.signed", "nmfae", "nmf") with added fields:
sigma2.used |
Estimated |
C.se, C.se.boot
|
Sandwich / bootstrap SEs for |
C.ci.lower, C.ci.upper
|
Bootstrap CIs. |
coefficients |
Data frame with Estimate, SE, BSE, z, p-value, CI. |
C.p.side |
P-value side used. |
This function is experimental. The interface may change in future versions.
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1), 45–55.
nmf.rrr.signed, nmf.rrr.inference
Fits nmf.rrr.signed with a paired decoder/encoder
rank () across a range of ranks and reports
r.squared, the effective rank (of the latent encoding ),
and the element-wise CV error sigma.ecv, with the same concise
plot as nmfkc.rank. For a full grid use
nmf.rrr.signed.ecv.
nmf.rrr.signed.rank( Y1, Y2 = Y1, rank1 = 1:5, detail = c("full", "fast"), plot = TRUE, ... )nmf.rrr.signed.rank( Y1, Y2 = Y1, rank1 = 1:5, detail = c("full", "fast"), plot = TRUE, ... )
Y1 |
Endogenous matrix ( |
Y2 |
Exogenous matrix; defaults to |
rank1 |
Integer vector of (paired) ranks to evaluate (both bases use
the same value). Legacy |
detail |
|
plot |
Logical; draw the diagnostics plot (default |
... |
Passed on to |
A list with rank.best and criteria
(rank, effective.rank, effective.rank.ratio,
r.squared, sigma.ecv).
Roy, O., & Vetterli, M. (2007). The effective rank: A measure of
effective dimensionality. Proc. EUSIPCO, 606–610.
(effective.rank)
Wold, S. (1978). Cross-validatory estimation of the number of
components in factor and principal components models.
Technometrics, 20(4), 397–405. (sigma.ecv)
nmf.rrr.signed, nmf.rrr.signed.ecv,
nmfkc.rank
Replaces the default "Resp1", "Resp2", ... (response basis / X1 columns
and Cp/Cn/C rows) and "Cov1", "Cov2", ... (covariate basis / X2 rows and
Cp/Cn/C columns) with user-supplied labels. Propagates to coefficients
tables if present (e.g., from nmfae.inference).
nmf.rrr.signed.rename(x, X1.colnames = NULL, X2.rownames = NULL)nmf.rrr.signed.rename(x, X1.colnames = NULL, X2.rownames = NULL)
x |
An |
X1.colnames |
Character vector of length Q for decoder labels. |
X2.rownames |
Character vector of length R for encoder labels. |
The renamed object (same class).
This function is experimental. The interface may change in future versions.
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1), 45–55.
nmfkc fits a nonnegative matrix factorization with kernel covariates
under the tri-factorization model .
This function supports two major input modes:
Matrix Mode (Existing): nmfkc(Y=matrix, A=matrix, ...)
Formula Mode (New): nmfkc(formula=Y_vars ~ A_vars, data=df, rank=Q, ...)
The rank of the basis matrix can be specified using either the rank argument
(preferred for formula mode) or the hidden Q argument (for backward compatibility).
nmfkc( Y, A = NULL, rank = NULL, data, epsilon = 1e-04, maxit = 5000, verbose = TRUE, ... )nmfkc( Y, A = NULL, rank = NULL, data, epsilon = 1e-04, maxit = 5000, verbose = TRUE, ... )
Y |
Observation matrix (P x N), OR a formula object for Formula Mode.
In Formula Mode, use |
A |
Covariate matrix. Default is |
rank |
Integer. The rank of the basis matrix |
data |
Optional. A data frame from which variables in the formula should be taken. |
epsilon |
Positive convergence tolerance. From the 10th iteration
the fit stops when |
maxit |
Maximum number of iterations. |
verbose |
Logical. If |
... |
Additional arguments passed for fine-tuning regularization, initialization, constraints,
and output control. This includes the backward-compatible arguments
|
A list with components:
call |
The matched call, as captured by |
dims |
A character string summarizing the matrix dimensions of the model. |
runtime |
A character string summarizing the computation time. |
X |
Basis matrix. Column normalization depends on |
B |
Coefficient matrix |
XB |
Fitted values for |
C |
Parameter matrix. |
B.prob |
Soft-clustering probabilities derived from columns of |
B.cluster |
Hard-clustering labels (argmax over |
X.prob |
Row-wise soft-clustering probabilities derived from |
X.cluster |
Hard-clustering labels (argmax over |
X.anchor |
Present only when |
X.restriction |
The constraint that was applied to the columns of
|
A.attr |
List of attributes of the input covariate matrix |
formula.meta |
If fitted via Formula Mode, a list with |
objfunc |
Final objective value. |
objfunc.iter |
Objective values by iteration. The trace is trimmed to
|
iter |
Number of iterations performed. |
maxit |
The iteration cap that was in force. |
epsilon |
The convergence tolerance that was in force. |
epsilon.iter |
Relative change of the objective at the last step —
the quantity the stopping rule compares with |
objfunc.increases |
Number of steps at which the objective rose,
counted on the recorded (trimmed) trace. The multiplicative update is
monotone by itself, so a positive count means something outside it — a
restriction imposed by projection, or a numerical problem — is pushing
the iterate back, in which case the fit can oscillate and exhaust
|
converged |
Logical; whether |
r.squared |
|
r.squared.uncentered |
Uncentered |
r.squared.centered |
Row-mean centered |
method |
Character string indicating the optimization method used ( |
n.missing |
Number of missing (or zero-weighted) elements in |
n.total |
Total number of elements in |
rank |
The rank |
sigma |
The residual standard error, representing the typical deviation of the observed values |
mae |
Mean Absolute Error between |
criterion |
A list with
With
|
Satoh, K. (2026). Applying Non-negative Matrix Factorization with Covariates to the Longitudinal Data as Growth Curve Model. American Journal of Mathematical and Management Sciences. In press. arXiv:2403.05359. https://arxiv.org/abs/2403.05359
Satoh, K. (2025). Applying non-negative matrix factorization with covariates to multivariate time series data as a vector autoregression model. Japanese Journal of Statistics and Data Science. arXiv:2501.17446. doi:10.1007/s42081-025-00314-0
Satoh, K. (2026). Applying non-negative matrix factorization with covariates to label matrix for classification. Japanese Journal of Statistics and Data Science. doi:10.1007/s42081-026-00349-x
Ding, C., Li, T., Peng, W., & Park, H. (2006). Orthogonal Nonnegative Matrix Tri-Factorizations for Clustering. In Proc. 12th ACM SIGKDD (pp. 126–135). doi:10.1145/1150402.1150420
Roy, O., & Vetterli, M. (2007). The effective rank: A measure of effective dimensionality. In 15th European Signal Processing Conference (EUSIPCO) (pp. 606–610).
Gillis, N., & Vavasis, S. A. (2014). Fast and robust recursive algorithms for separable nonnegative matrix factorization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 36(4), 698–714. doi:10.1109/TPAMI.2013.226
nmfkc.inference (standard errors and p-values for
), summary.nmfkc, plot.nmfkc,
nmfkc.cv, nmfkc.ecv, nmfkc.rank,
nmfkc.kernel, nmfkc.ar,
nmfkc.DOT, predict.nmfkc
# Example 1. Matrix Mode (Existing) X <- cbind(c(1,0,1),c(0,1,0)) B <- cbind(c(1,0),c(0,1),c(1,1)) Y <- X %*% B rownames(Y) <- paste0("P",1:nrow(Y)) colnames(Y) <- paste0("N",1:ncol(Y)) print(X); print(B); print(Y) res <- nmfkc(Y,rank=2,epsilon=1e-6) res$X res$B # Example 2. Formula Mode set.seed(1) dummy_data <- data.frame(Y1=rpois(10,5), Y2=rpois(10,10), A1=abs(rnorm(10,5)), A2=abs(rnorm(10,3))) res_f <- nmfkc(Y1 + Y2 ~ A1 + A2, data=dummy_data, rank=2) # Example 3. Anchors: in Example 1, rows P1 and P2 of X carry one basis # each. X.anchor imposes such rows (here found by SPA), which makes the # factorization unique; X.init = "spa" only starts from them. res_a <- nmfkc(Y, rank = 2, X.anchor = "spa", epsilon = 1e-6) res_a$X.anchor # For symmetric NMF (Y approximated by X X^T or X C X^T), # use \code{\link{nmfkc.net}()} instead.# Example 1. Matrix Mode (Existing) X <- cbind(c(1,0,1),c(0,1,0)) B <- cbind(c(1,0),c(0,1),c(1,1)) Y <- X %*% B rownames(Y) <- paste0("P",1:nrow(Y)) colnames(Y) <- paste0("N",1:ncol(Y)) print(X); print(B); print(Y) res <- nmfkc(Y,rank=2,epsilon=1e-6) res$X res$B # Example 2. Formula Mode set.seed(1) dummy_data <- data.frame(Y1=rpois(10,5), Y2=rpois(10,10), A1=abs(rnorm(10,5)), A2=abs(rnorm(10,3))) res_f <- nmfkc(Y1 + Y2 ~ A1 + A2, data=dummy_data, rank=2) # Example 3. Anchors: in Example 1, rows P1 and P2 of X carry one basis # each. X.anchor imposes such rows (here found by SPA), which makes the # factorization unique; X.init = "spa" only starts from them. res_a <- nmfkc(Y, rank = 2, X.anchor = "spa", epsilon = 1e-6) res_a$X.anchor # For symmetric NMF (Y approximated by X X^T or X C X^T), # use \code{\link{nmfkc.net}()} instead.
nmfkc.ar generates the observation matrix and covariate matrix
corresponding to a specified autoregressive lag order.
If the input Y is a ts object, its time properties are preserved
in the "tsp_info" attribute, adjusted for the lag.
Additionally, the column names of Y and A are set to the corresponding time points.
nmfkc.ar(Y, degree = 1, intercept = TRUE)nmfkc.ar(Y, degree = 1, intercept = TRUE)
Y |
An observation matrix (P x N) or a |
degree |
The lag order of the autoregressive model. The default is 1. |
intercept |
Logical. If TRUE (default), an intercept term is added to the covariate matrix. |
A list containing:
Y |
Observation matrix (P x N_A) used for NMF. Includes adjusted |
A |
Covariate matrix (R x N_A) constructed according to the specified lag order. Includes adjusted |
A.columns |
Index matrix used to generate |
degree.max |
Maximum lag order. |
Satoh, K. (2025). Applying non-negative matrix factorization with covariates to multivariate time series data as a vector autoregression model. Japanese Journal of Statistics and Data Science. arXiv:2501.17446. doi:10.1007/s42081-025-00314-0
nmfkc, nmfkc.ar.degree.cv, nmfkc.ar.stationarity, nmfkc.ar.DOT
# Example using AirPassengers (ts object) d <- AirPassengers ar_data <- nmfkc.ar(d, degree = 2) dim(ar_data$Y) dim(ar_data$A) # Example using matrix input Y <- matrix(1:20, nrow = 2) ar_data <- nmfkc.ar(Y, degree = 1) ar_data$degree.max# Example using AirPassengers (ts object) d <- AirPassengers ar_data <- nmfkc.ar(d, degree = 2) dim(ar_data$Y) dim(ar_data$A) # Example using matrix input Y <- matrix(1:20, nrow = 2) ar_data <- nmfkc.ar(Y, degree = 1) ar_data$degree.max
nmfkc.ar.degree.cv selects the optimal lag order for an autoregressive model
by applying cross-validation over candidate degrees.
This function accepts both standard matrices (Variables x Time) and ts objects
(Time x Variables). ts objects are automatically transposed internally.
nmfkc.ar.degree.cv( Y, rank = 1, degree = 1:2, intercept = TRUE, plot = TRUE, ... )nmfkc.ar.degree.cv( Y, rank = 1, degree = 1:2, intercept = TRUE, plot = TRUE, ... )
Y |
Observation matrix |
rank |
Rank of the basis matrix. For backward compatibility,
|
degree |
A vector of candidate lag orders to be evaluated. |
intercept |
Logical. If TRUE (default), an intercept is added to the covariate matrix. |
plot |
Logical. If TRUE (default), a plot of the objective function values is drawn. |
... |
Additional arguments passed to |
A list with components:
degree |
The lag order that minimizes the cross-validation objective function. |
degree.max |
Maximum recommended lag order, computed as |
objfunc |
Objective function values for each candidate lag order. |
# Example using ts object directly d <- AirPassengers # Selection of degree (using ts object) # Note: Y is automatically transposed if it is a ts object nmfkc.ar.degree.cv(Y=d, rank=1, degree=11:14)# Example using ts object directly d <- AirPassengers # Selection of degree (using ts object) # Note: Y is automatically transposed if it is a ts object nmfkc.ar.degree.cv(Y=d, rank=1, degree=11:14)
Produces a Graphviz DOT script for visualizing autoregressive
NMF-with-covariates models constructed via nmfkc.ar + nmfkc.
The diagram displays three types of directed relationships:
Lagged predictors: ,
Current latent factors: ,
Optional intercept effects: Const -> X.
Importantly, no direct edges from lagged variables to current outputs
() are drawn, in accordance with the NMF-AR
formulation.
Each block of lagged variables is displayed in its own DOT subgraph (e.g., “T-1”, “T-2”, ...), while latent factor nodes and current-time outputs are arranged in separate clusters.
nmfkc.ar.DOT( result, degree = 1, intercept = any(colnames(result$C) == "(Intercept)"), threshold = 0.1, rankdir = "RL", fill = TRUE, weight_scale_xy = 5, weight_scale_lag = 5, weight_scale_int = 3, hide.isolated = TRUE )nmfkc.ar.DOT( result, degree = 1, intercept = any(colnames(result$C) == "(Intercept)"), threshold = 0.1, rankdir = "RL", fill = TRUE, weight_scale_xy = 5, weight_scale_lag = 5, weight_scale_int = 3, hide.isolated = TRUE )
result |
A fitted |
degree |
Maximum AR lag to visualize. |
intercept |
Logical; if |
threshold |
Minimum coefficient magnitude required to draw an edge. |
rankdir |
Graphviz rank direction (e.g., |
fill |
Logical; whether nodes are filled with color. |
weight_scale_xy |
Scaling factor for edges |
weight_scale_lag |
Scaling factor for lagged edges |
weight_scale_int |
Scaling factor for intercept edges. |
hide.isolated |
Logical. If |
A character string representing a Graphviz DOT file.
nmfkc.ar, nmfkc, plot.nmfkc.DOT
d <- AirPassengers ar_data <- nmfkc.ar(d, degree = 2) result <- nmfkc(ar_data$Y, ar_data$A, rank = 1) dot <- nmfkc.ar.DOT(result, degree = 2) cat(dot)d <- AirPassengers ar_data <- nmfkc.ar(d, degree = 2) result <- nmfkc(ar_data$Y, ar_data$A, rank = 1) dot <- nmfkc.ar.DOT(result, degree = 2) cat(dot)
This function is experimental. The interface may change in future versions: argument names, defaults and the contents of the returned object are not yet stable.
Reduces a fitted NMF-VAR to the operators that drive its
coefficient vector. Writing for the fitted
scores, the residual-free recursion is
and the carry the model's whole deterministic structure —
transition weights, propagation, roots, the long-run mean. They stay
readable where the matrices do not:
instead of .
Use plot.nmfkc.ar.latent to draw them,
nmfkc.ar.latent.inference to attach standard errors, and
nmfkc.ar.stationarity if the stationarity verdict is all that
is wanted.
nmfkc.ar.latent(x)nmfkc.ar.latent(x)
x |
A fitted |
The latent process is VARMA, not VAR. The display above drops the
observation residual. With the exact
identity is
i.e. a VARMA in general (verified to machine precision; dropping
the moving-average term leaves an error of order the residuals themselves).
is therefore the autoregressive operator of the fitted VAR
reduced to the latent space, not the coefficient matrix of a latent VAR that
the scores obey exactly. In particular an entry of should not be
read as Granger causality between latent conditions: the moving-average term
is common to all coordinates and is not conditioned on.
Direction convention: is the effect of
on — row = effect (at ),
column = cause (at ).
Identification. The factorization is unique only up to
, under which .
The eigenvalues, , and are therefore
invariant, while individual entries of are not. Column
normalization removes only the scale part of , and a separable fitted
does not certify the rest: with ,
and
the
basis is perfectly separable and column-normalized, yet ,
, the product is unchanged and is no
permutation. Read the entries of as descriptive and the invariants
as inferential; see nmfkc.ar.latent.inference.
Oscillation. A complex pair means a cycle, but a negative
real root also oscillates — it alternates in sign each period, i.e. with
period 2. cycle.period covers the first case and alternating
the second; non.oscillatory is TRUE only when neither occurs.
Positive real roots rule out oscillation but not overshoot: a
non-normal can send a coordinate up before it comes down, which is
exactly what the off-diagonal impulse responses of the Canada fit do. Note
also that for , and the roots are real by
construction (the discriminant is ), so in that setting
the absence of a cycle is a property of the design, not a finding.
An object of class "nmfkc.ar.latent": a list with
G |
List of |
G.sum |
|
eigenvalues |
The |
spectral.radius |
|
spectral.radius.sum |
|
stationary |
Logical; |
cycle.period |
Period |
alternating |
|
non.oscillatory |
|
theta0 |
Latent intercept |
mu.b, mu.y
|
Long-run means |
p.star |
Composition of the long-run mean,
|
colsum |
Column sums |
colsum.max |
|
df |
Effective dimension |
separability |
Per-basis maximum of the row-normalised |
X, Theta
|
The basis and the full coefficient matrix the quantities
above were computed from, carried along so that
|
X.restriction |
The constraint |
dims |
Named vector of |
Satoh, K. (2025). Applying non-negative matrix factorization with covariates to multivariate time series data as a vector autoregression model. Japanese Journal of Statistics and Data Science. arXiv:2501.17446. doi:10.1007/s42081-025-00314-0
plot.nmfkc.ar.latent,
nmfkc.ar.latent.inference,
nmfkc.ar.stationarity, nmfkc.ar,
nmfkc.ar.DOT
d <- AirPassengers ar_data <- nmfkc.ar(d, degree = 2) result <- nmfkc(ar_data$Y, ar_data$A, rank = 1) lat <- nmfkc.ar.latent(result) lat$G$lag1 # latent transition matrix of lag 1 lat$eigenvalues # all roots of the model lat # full reportd <- AirPassengers ar_data <- nmfkc.ar(d, degree = 2) result <- nmfkc(ar_data$Y, ar_data$A, rank = 1) lat <- nmfkc.ar.latent(result) lat$G$lag1 # latent transition matrix of lag 1 lat$eigenvalues # all roots of the model lat # full report
This function is experimental. The interface may change in future versions: argument names, defaults and the contents of the returned object are not yet stable.
nmfkc.ar.latent.inference attaches a wild bootstrap to the quantities
reported by nmfkc.ar.latent: standard errors and intervals for
the spectral radius, the long-run mean, the long-run composition and the
entries of the latent transition matrices .
nmfkc.ar.latent.inference(object, Y, A, ...)nmfkc.ar.latent.inference(object, Y, A, ...)
object |
An object from |
Y, A
|
The observation and covariate matrices the model was fitted on,
i.e. the |
... |
Additional arguments:
|
What is identified, and why that decides what can be inferred. The
factorization admits for invertible
, so an entry of means nothing across replicates;
nmfkc.inference deals with this by conditioning on the
estimated . This function instead re-estimates the basis in every
replicate, propagating its uncertainty. That splits the output in two:
Inferential. , the eigenvalues and the observed
long-run mean are invariant under the whole
family ( leaves the spectrum alone), so their
replicates need no alignment and their intervals mean what they say.
Conditional on identification. The entries
and the composition need to be
a permutation, and three conditions together deliver that: the column
scale fixed (X.restriction = "colSums"), an anchor row for
every column of , and every row of pinned.
The anchor rows force . The rows are then pinned one at a
time: row joins the pinned set when some covariate
(column of ) is non-zero only for bases in ,
including , and, for every , such a covariate has a
zero coefficient for . The first row pinned needs a pure
column (a covariate that loads on that basis alone); later rows only need
zeros that separate them from the rows already pinned, possibly in
different columns. Each pinned row makes a column of ,
so pinning every row makes a permutation. (This used to require a
pure column for every row, which is stronger: every fit that passed before
passes now.) Either condition alone is not enough: with
,
and
the basis is perfectly separable and column-normalized, yet ,
, is unchanged, is no
permutation, and differs between the two representatives —
no row of can be pinned there. Nor are zeros that merely
separate every pair of rows enough: with anchors and the covariate columns
, , , is admissible.
A coefficient counts as zero when it is at most kkt.tol or at most
pin.tol times the coefficient of the same covariate for the basis
being pinned (that ratio is the mixing the covariate still allows), and a
covariate whose coefficients are all numerical noise is ignored.
The mirror route swaps the two factors and is checked as well,
with nothing to choose: a pure covariate for every basis (forcing
), and the bases pinned one at a time through the rows of
– the first by an anchor row, the later ones by variables whose
zero loadings separate them from the bases already pinned (an entry of
counts as zero when it is at most 1 - sep.tol times the
loading of the basis being pinned). Either route suffices; identification$route
names the one that passed. identified.approx reports
the verdict and identification its parts. What a pass
means, to the precision of the fit: multiplicative updates leave small
positive values where the argument wants exact zeros, so the check
treats values within the tolerances as zero and reports how exact they
are. The margin of the route that passed is the largest such value
relative to the entry it was compared with, i.e. the rotation slack:
every equivalent fit mixes one basis into
another by at most about that much, so the entries of are
determined to about that relative precision. print() states it
(e.g. "the zeros it rests on hold to 0.0024, so an equivalent fit mixes
the bases by at most about 0.24\
part that fails. Two things lie outside the check: whether the
estimated zeros are zeros in the population (a question for their
intervals), and other fits that are only nearly as good (different
starting values can reach one). Replicates are permutation-aligned to
the original basis regardless.
X.restriction = "fixed" is the one clean case for — the
basis is not estimated, so .
prob.nonstationary is not a unit-root test. It is the
bootstrap tail probability , with replicates drawn
around the estimated model. Nothing imposes , so it
does not have the interpretation of a p-value for that hypothesis; a proper
test would have to resample under the null.
The zeros of are only meaningful under
misspecification. The strict-complementarity diagnostic reports the dual
margin over the zero
coefficients and the primal margin
over the positive ones. Both must be
bounded away from zero for the support to be recovered exactly and for the
active coefficients to attain the distribution that knowing the support would
give. If the model is correctly specified, makes
exactly, so strict complementarity fails at
every zero: then tends to a constant in
(one half for an isolated zero) instead of to one, and no
resampling scheme is consistent there. A small kkt$ratio is therefore
not a numerical defect — it says the zeros are not estimating a sign
restriction. What makes them estimate something is misspecification: they
recover the coordinates whose unconstrained population coefficient is
strictly negative.
Bootstrap validity at the boundary. A naive nonparametric bootstrap
is inconsistent for a parameter on the boundary. The scheme used here
resamples the residuals and re-imposes the non-negativity constraint in every
re-fit, which is the practical remedy. wild.unit = "block" in
addition keeps the serial dependence of the residuals within a block; the
lagged design itself stays fixed (see below).
Bootstrap p-values invert the centred distribution. Every replicate
of a non-negative coefficient is , so comparing the raw replicates
with zero would return for every entry. The reported
G.p.value is the basic-bootstrap two-sided obtained from
.
Resampling scheme. With and residuals
, each replicate forms and
re-fits nmfkc at the same rank. carries one multiplier
per time point (wild.unit = "column", the default, which preserves the
contemporaneous covariance between series), one per block of consecutive
time points (wild.unit = "block", below), or one per cell.
With i.i.d. multipliers over time the intervals ignore residual
autocorrelation and may be too narrow when it is present. Negative entries
of are truncated at 0 to keep the response non-negative; the
truncation rate is reported, and a large one is a warning sign for the whole
procedure.
Block wild bootstrap (wild.unit = "block"). One multiplier
is drawn per block of consecutive time points and shared by all series, where
is the number of time points, ncol(Y) (after the lags are
removed): for , 5 for 100, 8 for 500, 10 for 1000,
and (the same draws as "column") for . The
block grid starts at a random offset in ,
drawn anew in each replicate, so the first and last blocks may be shorter
and no time point is always at a block edge. Residuals within a block keep
their serial dependence, which matters when the residuals are autocorrelated
(for example when the lag order is too small); for a correctly specified VAR
the default "column" suffices. is the usual rate
for the block length of block resampling. The length is fixed by ; it
cannot be set and is not stored in the result.
This is a fixed-design bootstrap: the covariate matrix is
held at the observed lags rather than rebuilt from , so the
intervals are conditional on the observed past and do not propagate the
randomness of the lagged design itself. They may therefore be somewhat
optimistic.
An object of class "nmfkc.ar.latent.inference": a list with
coefficients |
Data frame of the |
G, G.se, G.ci.lower, G.ci.upper, G.p.value
|
The same quantities as
lists of |
spectral.radius |
The radius of the |
prob.nonstationary |
Bootstrap tail probability
|
complex, complex.frac
|
Whether the point estimate has complex roots, and the fraction of replicates that do. |
mu.y |
With |
p.star |
With |
kkt |
Strict-complementarity diagnostic: |
identified.approx |
Logical: is the factorization unique up to a
permutation of the bases, to the precision of the fit? |
identification |
The parts of the diagnostic: for the first route
|
refit.args |
The arguments the bootstrap re-fits actually used, so the replicates can be checked against the original estimator. |
Plus the bookkeeping entries wild.B, wild.B.requested,
n.fail, truncate.rate, wild.unit, wild.dist,
wild.level, dims.
Satoh, K. (2025). Applying non-negative matrix factorization with covariates to multivariate time series data as a vector autoregression model. Japanese Journal of Statistics and Data Science. arXiv:2501.17446. doi:10.1007/s42081-025-00314-0
nmfkc.ar.latent, nmfkc.ar.stationarity,
nmfkc.ar, nmfkc.inference
set.seed(1) Y <- matrix(abs(rnorm(4 * 40)) + 1, 4, 40) ar <- nmfkc.ar(Y, degree = 1) fit <- nmfkc(ar$Y, ar$A, rank = 2, verbose = FALSE) lat <- nmfkc.ar.latent(fit) nmfkc.ar.latent.inference(lat, ar$Y, ar$A, wild.B = 50)set.seed(1) Y <- matrix(abs(rnorm(4 * 40)) + 1, 4, 40) ar <- nmfkc.ar(Y, degree = 1) fit <- nmfkc(ar$Y, ar$A, rank = 2, verbose = FALSE) lat <- nmfkc.ar.latent(fit) nmfkc.ar.latent.inference(lat, ar$Y, ar$A, wild.B = 50)
nmfkc.ar.predict computes multi-step-ahead forecasts for a fitted NMF-VAR model
using recursive forecasting.
If the fitted model contains time series property information (from nmfkc.ar),
the forecasted values will have appropriate time-based column names.
nmfkc.ar.predict(x, Y, degree = NULL, n.ahead = 1)nmfkc.ar.predict(x, Y, degree = NULL, n.ahead = 1)
x |
An object of class |
Y |
The historical observation matrix used for fitting (or at least the last |
degree |
Optional integer. Lag order (D). If |
n.ahead |
Integer (>=1). Number of steps ahead to forecast. |
A list with components:
pred |
A |
time |
A numeric vector of future time points corresponding to the columns of |
# Forecast AirPassengers d <- AirPassengers ar_data <- nmfkc.ar(d, degree = 2) result <- nmfkc(ar_data$Y, ar_data$A, rank = 1) pred <- nmfkc.ar.predict(result, Y = matrix(d, nrow = 1), degree = 2, n.ahead = 3) pred$pred# Forecast AirPassengers d <- AirPassengers ar_data <- nmfkc.ar(d, degree = 2) result <- nmfkc(ar_data$Y, ar_data$A, rank = 1) pred <- nmfkc.ar.predict(result, Y = matrix(d, nrow = 1), degree = 2, n.ahead = 3) pred$pred
Reports whether a fitted NMF-VAR is stationary, and how cheaply that could be
decided. Three numbers are returned, in increasing order of sharpness:
colsum.max , the largest column sum of
, certifies stationarity when below 1 with
no eigenvalue computation at all (a sufficient condition only, so a value
above 1 is merely inconclusive); spectral.radius.sum
is below 1 exactly when the model is stationary (and
is bracketed by and ); and
spectral.radius is the spectral radius of the companion matrix itself.
Note that the are weighted by , i.e.
colSums(X) %*% Lambda. Using the raw column sums of
is correct only when and can
certify stationarity for a non-stationary fit otherwise.
This function answers the yes/no question only. For the dynamics behind the
number — the transition matrices, the roots, the impulse responses, the
long-run mean — use nmfkc.ar.latent, and for the standard
error of use nmfkc.ar.latent.inference.
nmfkc.ar.stationarity(x, method = c("latent", "companion"))nmfkc.ar.stationarity(x, method = c("latent", "companion"))
x |
A fitted |
method |
Route used for |
An object of class "nmfkc.ar.stationarity": a list with
spectral.radius |
Spectral radius of the companion matrix (route
selected by |
stationary |
Logical; |
method |
The route actually used. |
spectral.radius.sum |
|
colsum |
Column sums |
colsum.max |
|
dims |
Named vector of |
Satoh, K. (2025). Applying non-negative matrix factorization with covariates to multivariate time series data as a vector autoregression model. Japanese Journal of Statistics and Data Science. arXiv:2501.17446. doi:10.1007/s42081-025-00314-0
nmfkc.ar.latent,
nmfkc.ar.latent.inference, nmfkc.ar
# Check stationarity of fitted AR model d <- AirPassengers ar_data <- nmfkc.ar(d, degree = 2) result <- nmfkc(ar_data$Y, ar_data$A, rank = 1) st <- nmfkc.ar.stationarity(result) st$spectral.radius st$stationary st # full report# Check stationarity of fitted AR model d <- AirPassengers ar_data <- nmfkc.ar(d, degree = 2) result <- nmfkc(ar_data$Y, ar_data$A, rank = 1) st <- nmfkc.ar.stationarity(result) st$spectral.radius st$stationary st # full report
Prototype of Tan & Fevotte's (2013) ARD-NMF (Euclidean /
). Unlike the cross-validation (nmfkc.ecv,
nmfkc.bicv) and stability (nmfkc.consensus)
engines that scan a range of ranks, ARD fits NMF once at
an over-complete rank and prunes automatically: each component
carries a relevance weight with an inverse-gamma prior;
the multiplicative updates gain a penalty (L2 /
half-normal prior) or (L1 / exponential), which
drives unsupported components to zero. The number of surviving
components is the estimated rank. Covariates are ignored (plain NMF).
This is a model-based point estimate: the result depends on the prior, the starting rank and the random initialization, and can vary run to run. Use it as a complement to the CV / consensus engines, not a sole criterion.
nmfkc.ard(Y, rank = NULL, nrun = 10, plot = FALSE, ...)nmfkc.ard(Y, rank = NULL, nrun = 10, plot = FALSE, ...)
Y |
Observation matrix ( |
rank |
Over-complete starting rank |
nrun |
Number of random-initialization restarts (default |
plot |
Logical; draw the relevance bar plot. |
... |
Advanced options, rarely needed (defaults in parentheses):
|
Relation to Tan & Fevotte (2013). The update equations reproduce
the paper's / ARD-NMF for the Euclidean
() case exactly: the multiplicative penalties, the
closed-form update and the
constant (their Eq. 33). Two deliberate simplifications keep it
practical: only is implemented (the paper covers the general
-divergence), and the default b is an empirical
per-component energy scale that reliably avoids
the winner-take-all collapse, rather than the paper's method-of-moments value
(their Eq. 38).
An object of class "nmfkc.ard": a list with
rank (estimated = mode over restarts), rank.runs (the
per-run estimates), relevance (representative run, descending),
lambda, W, H (ordered by relevance),
rank.init, prior, nrun, objfunc (final
objective value of the representative run) and objfunc.iter (its
per-iteration trajectory).
V. Y. F. Tan and C. Fevotte (2013). Automatic relevance determination in nonnegative matrix factorization with the beta-divergence. IEEE TPAMI 35(7):1592–1605. doi:10.1109/TPAMI.2012.240.
nmfkc.ecv, nmfkc.bicv,
nmfkc.consensus, nmfkc.rank
set.seed(1) X <- matrix(abs(rnorm(40 * 3)), 40, 3) B <- matrix(abs(rnorm(3 * 60)), 3, 60) ar <- nmfkc.ard(X %*% B, rank = 10) # over-complete start ar$rank # ~ 3 surviving components plot(ar)set.seed(1) X <- matrix(abs(rnorm(40 * 3)), 40, 3) B <- matrix(abs(rnorm(3 * 60)), 3, 60) ar <- nmfkc.ard(X %*% B, rank = 10) # over-complete start ar$rank # ~ 3 surviving components plot(ar)
Owen & Perry's (2009) bi-cross-validation (BCV) for choosing the NMF
rank. A lightweight CV engine in the spirit of nmfkc.ecv:
it returns the held-out error per rank and nothing more (no plot, no
table) – pass the result to which.min(sigma), or build your own
diagnostics.
Unlike the element-wise CV of nmfkc.ecv (which holds out
scattered entries and refits with weights), BCV holds out a whole
row-block and column-block simultaneously: the model is fit only
on the retained block , and the held-out block is predicted
by folding the held-out rows/columns onto the fixed -factors via
non-negative regression (). Because the held-out
rows and columns never enter the fit, there is no information leakage.
Covariates are ignored (plain NMF). The recommended setting is to leave
out roughly half the rows and half the columns (nfolds = 2).
nmfkc.bicv(Y, rank = 1:3, ...)nmfkc.bicv(Y, rank = 1:3, ...)
Y |
Observation matrix ( |
rank |
Integer vector of ranks to evaluate. |
... |
Advanced options, rarely needed (defaults in parentheses):
|
Each fold keeps about of the rows
and columns, so the retained block must have more than
rank rows and columns. The largest testable rank is
therefore about ; with
nfolds = 2 this is roughly . Ranks above
this return NA and trigger a warning that
names the limit and the nfolds (or nmfkc.ecv)
that would reach the requested ranks. Raising nfolds lifts the
limit at the cost of a smaller hold-out and more compute
( full fits per rank).
A list (cf.\ nmfkc.ecv) with:
objfunc |
Held-out mean squared error for each rank. |
sigma |
Its square root (RMSE) for each rank. |
rank |
The evaluated rank vector. |
nfolds |
The number of folds used. |
A. B. Owen and P. O. Perry (2009). Bi-cross-validation of the SVD and the nonnegative matrix factorization. Ann. Appl. Stat. 3(2):564–594. doi:10.1214/08-AOAS227.
nmfkc.rank, nmfkc.ecv,
nmfkc.consensus, nmfkc.ard
## rank-3 non-negative data; bi-CV needs enough kept rows/cols per ## fold (> rank), so use a matrix with ample dimensions. set.seed(1) X <- matrix(abs(rnorm(30 * 3)), 30, 3) B <- matrix(abs(rnorm(3 * 40)), 3, 40) bv <- nmfkc.bicv(X %*% B, rank = 1:6) # nfolds = 2 (Owen & Perry) by default bv$sigma # held-out RMSE per rank bv$rank[which.min(bv$sigma)]## rank-3 non-negative data; bi-CV needs enough kept rows/cols per ## fold (> rank), so use a matrix with ample dimensions. set.seed(1) X <- matrix(abs(rnorm(30 * 3)), 30, 3) B <- matrix(abs(rnorm(3 * 40)), 3, 40) bv <- nmfkc.bicv(X %*% B, rank = 1:6) # nfolds = 2 (Owen & Perry) by default bv$sigma # held-out RMSE per rank bv$rank[which.min(bv$sigma)]
nmfkc.class converts a categorical or factor vector into a class matrix
(one-hot encoded representation), where each row corresponds to a category
and each column corresponds to an observation.
nmfkc.class(x)nmfkc.class(x)
x |
A categorical vector or a factor. |
A binary matrix with one row per unique category and one column per observation. Each column has exactly one entry equal to 1, indicating the category of the observation.
# Example. Y <- nmfkc.class(iris$Species) Y[,1:6]# Example. Y <- nmfkc.class(iris$Species) Y[,1:6]
The bioinformatics-standard stability approach to choosing the NMF rank.
A lightweight engine like nmfkc.ecv / nmfkc.bicv:
it returns one stability score per rank and nothing more.
For each rank, NMF is run nrun times from different random
initializations (X.init = "runif"). Each run gives a hard
clustering of the samples (the column of the coefficient
matrix). Averaging the connectivity matrices (1 if two
samples share a cluster) over the runs yields the consensus
matrix; its crispness measures how reproducible the clustering is. Two
summaries are reported per rank:
cophenetic: the cophenetic correlation coefficient (CPCC)
of the consensus matrix (Brunet et al. 2004). Close to 1 = stable.
dispersion: the Kim & Park (2007) dispersion
; 1 when
every consensus entry is exactly 0 or 1 (perfectly crisp).
pac: the Proportion of Ambiguous Clustering
(Senbabaoglu et al. 2014) – the fraction of off-diagonal consensus
entries falling in the ambiguous interval pac.range
(default ). Lower is better (less
ambiguity); a more sensitive readout than cophenetic, which
tends to saturate.
Unlike the cross-validation engines (where the rank minimizes the error), here a good rank maximizes stability, or is the largest rank before it drops.
nmfkc.consensus( Y, A = NULL, rank = 2:4, nrun = 30, keep.consensus = FALSE, ... )nmfkc.consensus( Y, A = NULL, rank = 2:4, nrun = 30, keep.consensus = FALSE, ... )
Y |
Observation matrix ( |
A |
Optional covariate matrix passed to |
rank |
Integer vector of ranks to evaluate ( |
nrun |
Number of random-initialization runs per rank (default 30). |
keep.consensus |
Logical; if |
... |
Advanced options, rarely needed (defaults in parentheses):
|
An object of class "nmfkc.consensus" (a list) with:
cophenetic |
Cophenetic correlation coefficient for each rank. |
dispersion |
Dispersion coefficient ( |
pac |
Proportion of Ambiguous Clustering ( |
rank |
The evaluated rank vector. |
nrun |
Number of runs per rank. |
consensus |
List of consensus matrices, or |
It has print.nmfkc.consensus and
plot.nmfkc.consensus (type = "criteria" /
"heatmap") methods.
Brunet, J.-P., Tamayo, P., Golub, T. R., Mesirov, J. P. (2004). Metagenes and molecular pattern discovery using matrix factorization. PNAS 101(12):4164–4169. doi:10.1073/pnas.0308531101. Kim, H., Park, H. (2007). Sparse non-negative matrix factorizations ... Bioinformatics 23(12):1495–1502. Senbabaoglu, Y., Michailidis, G., Li, J. Z. (2014). Critical limitations of consensus clustering in class discovery. Sci. Rep. 4:6207. doi:10.1038/srep06207.
nmfkc.rank, nmfkc.ecv,
nmfkc.bicv, nmfkc.ard,
nmf.cluster.criteria
Y <- t(as.matrix(iris[, 1:4])) cs <- nmfkc.consensus(Y, rank = 2:5, nrun = 20, keep.consensus = TRUE) cs # stability table per rank plot(cs) # type = "criteria": stability curves plot(cs, type = "heatmap") # all ranks, n2mfrow grid plot(cs, type = "heatmap", rank = 3) # one rank, with labelsY <- t(as.matrix(iris[, 1:4])) cs <- nmfkc.consensus(Y, rank = 2:5, nrun = 20, keep.consensus = TRUE) cs # stability table per rank plot(cs) # type = "criteria": stability curves plot(cs, type = "heatmap") # all ranks, n2mfrow grid plot(cs, type = "heatmap", rank = 3) # one rank, with labels
nmfkc.criterion computes the effective rank (raw and as the
broken-stick index) and the sample-clustering quality measures
(silhouette, CPCC, dist.cor) from a fitted nmfkc model.
This function can be called on a model that was fitted with
detail = "fast" or detail = "minimal" to compute the
full set of criteria afterwards.
nmfkc.criterion(object, Y, detail = c("full", "fast", "minimal"), ...)nmfkc.criterion(object, Y, detail = c("full", "fast", "minimal"), ...)
object |
An object of class |
Y |
The original observation matrix (P x N) used for fitting. |
detail |
Character string controlling the level of computation:
|
... |
Additional arguments: |
A list with components:
R-squared between Y and XB.
Residual standard deviation.
Mean absolute error.
Column-normalized coefficient matrix (soft-clustering probabilities).
Hard clustering labels (argmax of B.prob per column).
Row-normalized basis matrix.
Hard clustering labels per row of X.
Named list: effective.rank, effective.rank.index, silhouette, CPCC, dist.cor.
Y <- t(iris[, -5]) res <- nmfkc(Y, rank = 3, detail = "fast") crit <- nmfkc.criterion(res, Y) crit$criterion$silhouetteY <- t(iris[, -5]) res <- nmfkc(Y, rank = 3, detail = "fast") crit <- nmfkc.criterion(res, Y) crit$criterion$silhouette
nmfkc.cv performs k-fold cross-validation for the tri-factorization model
, where
is the observation matrix,
is the covariate (or kernel) matrix,
is the basis matrix,
is the parameter matrix, and
is the coefficient matrix ().
Given (and optionally ), and are fitted on each
training split and predictive performance is evaluated on the held-out split.
nmfkc.cv(Y, A = NULL, rank = 2, data, ...)nmfkc.cv(Y, A = NULL, rank = 2, data, ...)
Y |
Observation matrix, or a formula (see |
A |
Covariate matrix. If |
rank |
Rank of the basis matrix |
data |
A data frame (required when |
... |
Additional arguments controlling CV and the internal
|
A list with components:
objfuncMean loss per valid entry over all folds (MSE for method="EU").
sigmaResidual standard error (RMSE). Available only if method="EU"; on the same scale as Y.
objfunc.blockLoss for each fold.
blockVector of fold indices (1, …, div) assigned to each column of .
nmfkc, nmfkc.kernel.beta.cv, nmfkc.ar.degree.cv
# Example 1 (with explicit covariates): Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) res <- nmfkc.cv(Y, A, rank = 1) res$objfunc # Example 2 (kernel A and beta sweep): Y <- matrix(cars$dist, nrow = 1) U <- matrix(c(5, 10, 15, 20, 25), nrow = 1) V <- matrix(cars$speed, nrow = 1) betas <- 25:35/1000 obj <- numeric(length(betas)) for (i in seq_along(betas)) { A <- nmfkc.kernel(U, V, beta = betas[i]) obj[i] <- nmfkc.cv(Y, A, rank = 1, nfolds = 10)$objfunc } betas[which.min(obj)]# Example 1 (with explicit covariates): Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) res <- nmfkc.cv(Y, A, rank = 1) res$objfunc # Example 2 (kernel A and beta sweep): Y <- matrix(cars$dist, nrow = 1) U <- matrix(c(5, 10, 15, 20, 25), nrow = 1) V <- matrix(cars$speed, nrow = 1) betas <- 25:35/1000 obj <- numeric(length(betas)) for (i in seq_along(betas)) { A <- nmfkc.kernel(U, V, beta = betas[i]) obj[i] <- nmfkc.cv(Y, A, rank = 1, nfolds = 10)$objfunc } betas[which.min(obj)]
Compact print and (for rank sweeps) a score-vs-rank plot with a
“Best (Min)” marker, for the objects returned by
nmfkc.ecv, nmfkc.cv and nmfkc.bicv.
## S3 method for class 'nmfkc.ecv' print(x, ...) ## S3 method for class 'nmfkc.ecv' plot(x, main = "nmfkc element-wise CV", ylab = NULL, ...) ## S3 method for class 'nmfkc.bicv' print(x, ...) ## S3 method for class 'nmfkc.bicv' plot( x, main = "nmfkc bi-cross-validation", ylab = "sigma (held-out RMSE)", ... ) ## S3 method for class 'nmfkc.cv' print(x, ...) ## S3 method for class 'nmfkc.cv' plot( x, main = "nmfkc CV: per-fold held-out loss", ylab = "objfunc (block)", ... )## S3 method for class 'nmfkc.ecv' print(x, ...) ## S3 method for class 'nmfkc.ecv' plot(x, main = "nmfkc element-wise CV", ylab = NULL, ...) ## S3 method for class 'nmfkc.bicv' print(x, ...) ## S3 method for class 'nmfkc.bicv' plot( x, main = "nmfkc bi-cross-validation", ylab = "sigma (held-out RMSE)", ... ) ## S3 method for class 'nmfkc.cv' print(x, ...) ## S3 method for class 'nmfkc.cv' plot( x, main = "nmfkc CV: per-fold held-out loss", ylab = "objfunc (block)", ... )
x |
A |
... |
Passed to |
main, ylab
|
Plot labels. |
x, invisibly.
nmfkc.denormalize rescales a matrix with values in back to its
original scale using the column-wise minima and maxima of a reference matrix.
nmfkc.denormalize(x, ref = x)nmfkc.denormalize(x, ref = x)
x |
A numeric matrix (or vector) with values in |
ref |
A reference matrix used to obtain the original column-wise minima
and maxima. Must have the same number of columns as |
A numeric matrix with values transformed back to the original scale.
x <- nmfkc.normalize(iris[, -5]) x_recovered <- nmfkc.denormalize(x, iris[, -5]) apply(x_recovered - iris[, -5], 2, max)x <- nmfkc.normalize(iris[, -5]) x_recovered <- nmfkc.denormalize(x, iris[, -5]) apply(x_recovered - iris[, -5], 2, max)
Produces a Graphviz DOT script visualizing the structure of an NMF model
() or its simplified forms.
Supported visualization types:
"YX" — Standard NMF view: latent factors map to observations .
"YA" — Direct regression view: covariates map directly to
using the combined coefficient matrix .
"YXA" — Full tri-factorization: .
Edge widths are scaled by coefficient magnitude, and nodes with no edges above the threshold are omitted from the visualization.
nmfkc.DOT( result, type = c("YX", "YA", "YXA"), threshold = 0.01, C.signed = NULL, sig.level = 0.1, rankdir = "LR", fill = TRUE, weight_scale = 5, weight_scale_ax = weight_scale, weight_scale_xy = weight_scale, weight_scale_ay = weight_scale, Y.label = NULL, X.label = NULL, A.label = NULL, Y.title = "Observation (Y)", X.title = "Basis (X)", A.title = "Covariates (A)", hide.isolated = TRUE )nmfkc.DOT( result, type = c("YX", "YA", "YXA"), threshold = 0.01, C.signed = NULL, sig.level = 0.1, rankdir = "LR", fill = TRUE, weight_scale = 5, weight_scale_ax = weight_scale, weight_scale_xy = weight_scale, weight_scale_ay = weight_scale, Y.label = NULL, X.label = NULL, A.label = NULL, Y.title = "Observation (Y)", X.title = "Basis (X)", A.title = "Covariates (A)", hide.isolated = TRUE )
result |
The return value from |
type |
Character string specifying the visualization style:
one of |
threshold |
Minimum coefficient magnitude to display an edge. When
|
C.signed |
Logical or |
sig.level |
Significance level for filtering C edges when inference
results are available (i.e., |
rankdir |
Graphviz rank direction (e.g., |
fill |
Logical; whether nodes should be drawn with filled shapes. |
weight_scale |
Base scaling factor for edge widths. |
weight_scale_ax |
Scaling factor for edges |
weight_scale_xy |
Scaling factor for edges |
weight_scale_ay |
Scaling factor for edges |
Y.label |
Optional character vector for labels of Y nodes. |
X.label |
Optional character vector for labels of X (latent factor) nodes. |
A.label |
Optional character vector for labels of A (covariate) nodes. |
Y.title |
Cluster title for Y nodes. |
X.title |
Cluster title for X nodes. |
A.title |
Cluster title for A nodes. |
hide.isolated |
Logical. If |
A character string representing a Graphviz DOT script.
nmfkc, nmf.rrr.DOT, nmf.ffb.DOT,
nmfkc.ar.DOT, plot.nmfkc.DOT
Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) dot <- nmfkc.DOT(result) cat(dot)Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) dot <- nmfkc.DOT(result) cat(dot)
nmfkc.ecv performs k-fold cross-validation by randomly holding out
individual elements of the data matrix (element-wise), assigning them a
weight of 0 via Y.weights, and evaluating the reconstruction error on
those held-out elements.
This method (also known as Wold's CV) is theoretically robust for determining
the optimal rank (Q) in NMF. This function supports vector input for Q,
allowing simultaneous evaluation of multiple ranks on the same folds.
For symmetric (network) data use nmfkc.net.ecv, which
creates upper-triangle folds to prevent information leakage through
the symmetric entries . Passing the old
Y.symmetric argument here is no longer supported and stops
with a redirect message.
nmfkc.ecv(Y, A = NULL, rank = 1:3, data, ...)nmfkc.ecv(Y, A = NULL, rank = 1:3, data, ...)
Y |
Observation matrix, or a formula (see |
A |
Covariate matrix. Ignored when |
rank |
Vector of ranks to evaluate (e.g., 1:5). For backward compatibility,
|
data |
A data frame (required when |
... |
Additional arguments passed to |
A list with components:
objfunc |
Numeric vector containing the Mean Squared Error (MSE) for
each Q, named |
sigma |
Numeric vector containing the Residual Standard Error (RMSE) for each Q. Only available if method="EU". |
objfunc.fold |
List of length equal to Q vector. Each element contains the MSE values for the k folds. |
folds |
A list of length |
Wold, S. (1978). Cross-validatory estimation of the number of
components in factor and principal components models.
Technometrics, 20(4), 397–405.
doi:10.1080/00401706.1978.10489693
Owen, A. B., & Perry, P. O. (2009). Bi-cross-validation of the SVD
and the nonnegative matrix factorization. Ann. Appl. Stat.
3(2), 564–594. doi:10.1214/08-AOAS227 (cross-validation of the
NMF rank; see also nmfkc.bicv).
nmfkc, nmfkc.cv; other rank-selection
criteria: nmfkc.rank, nmfkc.bicv,
nmfkc.consensus, nmfkc.ard.
# Element-wise CV to select rank Y <- t(iris[1:30, 1:4]) res <- nmfkc.ecv(Y, rank = 1:2, nfolds = 3) res$objfunc# Element-wise CV to select rank Y <- t(iris[1:30, 1:4]) res <- nmfkc.ecv(Y, rank = 1:2, nfolds = 3) res$objfunc
nmfkc.inference performs statistical inference on the parameter matrix
() from a fitted nmfkc model, conditional on
the estimated basis matrix .
Under the working model where
,
inference is conducted via sandwich covariance estimation and
one-step wild bootstrap with non-negative projection.
nmfkc.inference(object, Y, A = NULL, wild.bootstrap = TRUE, ...)nmfkc.inference(object, Y, A = NULL, wild.bootstrap = TRUE, ...)
object |
An object of class |
Y |
Observation matrix (P x N). Must match the data used in |
A |
Covariate matrix (K x N). Default is |
wild.bootstrap |
Logical. If |
... |
Additional arguments:
|
An object of class c("nmfkc.inference", "nmfkc"), inheriting all
components from the input object, with additional inference components:
sigma2.used |
Estimated |
C.se |
Sandwich standard errors for |
C.se.boot |
Bootstrap standard errors for |
C.ci.lower |
Lower CI bounds for |
C.ci.upper |
Upper CI bounds for |
coefficients |
Data frame with one row per element of |
C.p.side |
P-value type used. |
wild.unit |
Multiplier unit of the |
identified.approx |
Logical, an approximate diagnostic of whether the
factorization is unique up to a permutation of the bases. The intervals
are conditional on |
identification |
Its parts: |
Satoh, K. (2026). Wild Bootstrap Inference for Non-Negative Matrix Factorization with Random Effects. arXiv:2603.01468. https://arxiv.org/abs/2603.01468
nmfkc, summary.nmfkc.inference
Y <- matrix(cars$dist, nrow = 1) A <- rbind(intercept = 1, speed = cars$speed) result <- nmfkc(Y, A, rank = 1) result2 <- nmfkc.inference(result, Y, A) summary(result2)Y <- matrix(cars$dist, nrow = 1) A <- rbind(intercept = 1, speed = cars$speed) result <- nmfkc(Y, A, rank = 1) result2 <- nmfkc.inference(result, Y, A) summary(result2)
nmfkc.kernel constructs a kernel matrix from covariate matrices.
It supports Gaussian, Exponential, Periodic, Linear, Normalized Linear, and Polynomial kernels.
nmfkc.kernel( U, V = NULL, kernel = c("Gaussian", "Exponential", "Periodic", "Linear", "NormalizedLinear", "Polynomial"), ... )nmfkc.kernel( U, V = NULL, kernel = c("Gaussian", "Exponential", "Periodic", "Linear", "NormalizedLinear", "Polynomial"), ... )
U |
Covariate matrix |
V |
Covariate matrix |
kernel |
Kernel function to use. Default is |
... |
Additional arguments passed to the specific kernel function (e.g., |
Kernel matrix .
Satoh, K. (2026). Applying Non-negative Matrix Factorization with Covariates to the Longitudinal Data as Growth Curve Model. American Journal of Mathematical and Management Sciences. In press. arXiv:2403.05359. https://arxiv.org/abs/2403.05359
nmfkc.kernel.gaussian, nmfkc.cv
# Example. Y <- matrix(cars$dist,nrow=1) U <- matrix(c(5,10,15,20,25),nrow=1) V <- matrix(cars$speed,nrow=1) A <- nmfkc.kernel(U,V,beta=28/1000) dim(A) result <- nmfkc(Y,A,rank=1) plot(as.vector(V),as.vector(Y)) lines(as.vector(V),as.vector(result$XB),col=2,lwd=2)# Example. Y <- matrix(cars$dist,nrow=1) U <- matrix(c(5,10,15,20,25),nrow=1) V <- matrix(cars$speed,nrow=1) A <- nmfkc.kernel(U,V,beta=28/1000) dim(A) result <- nmfkc(Y,A,rank=1) plot(as.vector(V),as.vector(Y)) lines(as.vector(V),as.vector(result$XB),col=2,lwd=2)
nmfkc.kernel.beta.cv selects the optimal beta parameter of the kernel function by applying cross-validation over a set of candidate values.
nmfkc.kernel.beta.cv(Y, rank = 2, U, V = NULL, beta = NULL, plot = TRUE, ...)nmfkc.kernel.beta.cv(Y, rank = 2, U, V = NULL, beta = NULL, plot = TRUE, ...)
Y |
Observation matrix |
rank |
Rank of the basis matrix. |
U |
Covariate matrix |
V |
Covariate matrix |
beta |
A numeric vector of candidate kernel parameters to evaluate via cross-validation. |
plot |
Logical. If TRUE (default), plots the objective function values for each candidate |
... |
Additional arguments passed to |
A list with components:
beta |
The beta value that minimizes the cross-validation objective function. |
objfunc |
Objective function values for each candidate |
nmfkc.kernel.gaussian, nmfkc.kernel.beta.nearest.med,
nmfkc.kernel
# Example. Y <- matrix(cars$dist,nrow=1) U <- matrix(c(5,10,15,20,25),nrow=1) V <- matrix(cars$speed,nrow=1) nmfkc.kernel.beta.cv(Y,rank=1,U,V,beta=25:30/1000) A <- nmfkc.kernel(U,V,beta=28/1000) result <- nmfkc(Y,A,rank=1) plot(as.vector(V),as.vector(Y)) lines(as.vector(V),as.vector(result$XB),col=2,lwd=2)# Example. Y <- matrix(cars$dist,nrow=1) U <- matrix(c(5,10,15,20,25),nrow=1) V <- matrix(cars$speed,nrow=1) nmfkc.kernel.beta.cv(Y,rank=1,U,V,beta=25:30/1000) A <- nmfkc.kernel(U,V,beta=28/1000) result <- nmfkc(Y,A,rank=1) plot(as.vector(V),as.vector(Y)) lines(as.vector(V),as.vector(result$XB),col=2,lwd=2)
Computes a data-driven reference scale for the Gaussian/RBF kernel from covariates
using a robust "median nearest-neighbor (or nearest-landmark) distance" heuristic,
and returns the corresponding kernel parameter .
The Gaussian/RBF kernel is assumed to be written in the form
hence . This function first estimates a typical distance
scale by the median of distances, then sets .
If Uk is NULL, is estimated as the median of
nearest-neighbor distances within U (excluding self-distance).
If Uk is provided, is estimated as the median of
nearest-landmark distances from each sample in U to its closest landmark in Uk.
To control memory usage for large N (and M), distances are computed in blocks.
Optionally, columns of U can be randomly subsampled via sample.size to reduce cost.
nmfkc.kernel.beta.nearest.med( U, Uk = NULL, block.size = 1000, block.size.Uk = 2000, sample.size = NULL, ... )nmfkc.kernel.beta.nearest.med( U, Uk = NULL, block.size = 1000, block.size.Uk = 2000, sample.size = NULL, ... )
U |
A numeric matrix of covariates ( |
Uk |
An optional numeric matrix of landmarks ( |
block.size |
Integer. Number of columns of |
block.size.Uk |
Integer. Number of columns of |
sample.size |
Integer or |
... |
Additional arguments. Hidden option |
Candidate grid:
Along with beta, the function returns beta_candidates, a
logarithmic grid suitable for cross-validation. The grid is symmetric on
the bandwidth scale around :
and since , this corresponds to
.
The grid of values can be customized through the hidden argument
candidates (passed via ...):
"7points" (default):
(7 candidates spanning one decade, matches the grid used in the
RFF-NMF research memo).
"4points": yielding
(the legacy short grid).
A numeric vector: user-specified values. The grid returned
is .
Prior to version 0.6.8, the grid depended on whether Uk was
supplied (4 candidates for Uk = NULL, 7 for supplied Uk).
The current implementation unifies both branches via candidates.
Notes:
When Uk is identical to U, the function detects this case and excludes
self-distances (distance 0) to avoid .
sample.size performs random subsampling without setting a seed. For reproducible
results, set set.seed() before calling this function.
A list with elements:
beta: Estimated kernel parameter .
beta_candidates: Numeric vector of candidate values (logarithmic grid)
intended for cross-validation.
dist_median: The estimated distance scale (median of nearest-neighbor
or nearest-landmark distances).
block.size.used: The effective block size(s) used. Either a scalar (no Uk) or
a named vector c(U=..., Uk=...) when Uk is provided.
sample.size.used: The number of columns of U actually used (after subsampling).
uk_is_u: Logical flag indicating whether Uk was detected as identical to U
(only returned when Uk is provided).
nmfkc.kernel.gaussian, nmfkc.kernel.beta.cv
# Basic (nearest-neighbor within U) U <- matrix(runif(20), nrow = 2) beta_info <- nmfkc.kernel.beta.nearest.med(U) beta0 <- beta_info$beta betas <- beta_info$beta_candidates # With landmarks (nearest-landmark distances) Uk <- matrix(runif(10), nrow = 2) beta_info2 <- nmfkc.kernel.beta.nearest.med(U, Uk)# Basic (nearest-neighbor within U) U <- matrix(runif(20), nrow = 2) beta_info <- nmfkc.kernel.beta.nearest.med(U) beta0 <- beta_info$beta betas <- beta_info$beta_candidates # With landmarks (nearest-landmark distances) Uk <- matrix(runif(10), nrow = 2) beta_info2 <- nmfkc.kernel.beta.nearest.med(U, Uk)
nmfkc.kernel.gaussian constructs a Gaussian (RBF) kernel matrix from covariate matrices.
The kernel is defined as .
When V contains NA values, two methods are available via na.method:
"pds"Partial Distance Strategy. Computes the kernel using only observed (non-NA) rows,
with beta adjusted by where is the total number of rows
and is the number of observed rows.
"egk"Expected Gaussian Kernel (Mesquita et al., 2019). Uses a Gaussian Mixture Model (GMM)
to estimate the conditional distribution of missing values given observed values,
then computes the expected kernel value via a Gamma approximation.
Requires gmm.means, gmm.sigmas, and gmm.weights passed through ....
nmfkc.kernel.gaussian( U, V = NULL, beta = 0.5, na.method = c("pds", "egk"), ... )nmfkc.kernel.gaussian( U, V = NULL, beta = 0.5, na.method = c("pds", "egk"), ... )
U |
Covariate matrix |
V |
Covariate matrix |
beta |
Bandwidth parameter for the Gaussian kernel. Default is |
na.method |
Method for handling |
... |
Additional arguments for EGK method:
|
Kernel matrix .
Mesquita, D., Gomes, J. P., & Rodrigues, L. R. (2019). Gaussian kernels for incomplete data. Applied Soft Computing, 77, 356–365.
nmfkc.kernel, nmfkc.kernel.beta.cv, nmfkc.kernel.beta.nearest.med
U <- matrix(c(5,10,15,20,25),nrow=1) V <- matrix(1:25,nrow=1) A <- nmfkc.kernel.gaussian(U,V,beta=28/1000) dim(A) # PDS example: V with NA in first row U2 <- matrix(rnorm(20), nrow=2) V2 <- matrix(rnorm(10), nrow=2) V2[1, c(2,4)] <- NA A2 <- nmfkc.kernel.gaussian(U2, V2, beta=0.5, na.method="pds")U <- matrix(c(5,10,15,20,25),nrow=1) V <- matrix(1:25,nrow=1) A <- nmfkc.kernel.gaussian(U,V,beta=28/1000) dim(A) # PDS example: V with NA in first row U2 <- matrix(rnorm(20), nrow=2) V2 <- matrix(rnorm(10), nrow=2) V2[1, c(2,4)] <- NA A2 <- nmfkc.kernel.gaussian(U2, V2, beta=0.5, na.method="pds")
Builds the two summary matrices that nmfkc needs from
the Nystr\"om kernel covariate matrix
(, )
without ever forming itself:
The input U is processed in column blocks; for each block the
kernel values are computed with
nmfkc.kernel, added into and , and
discarded. Peak memory is
instead of
.
The returned object is passed to nmfkc as its A
argument. The multiplicative updates are unchanged, since they only
ever use and ; the fit equals the one obtained from the
explicit matrix nmfkc.kernel(V, U, beta = beta) up to
floating-point summation order. This is the large- form of the
NMF-LAB Kernel design with Nystr\"om approximation (Satoh 2026).
nmfkc.kernel.gram(Y, U, V, beta = NULL, block.size = NULL, ...)nmfkc.kernel.gram(Y, U, V, beta = NULL, block.size = NULL, ...)
Y |
Non-negative |
U |
A |
V |
Landmarks. Either a |
beta |
Positive scalar. Gaussian kernel bandwidth
|
block.size |
Integer. Number of columns of |
... |
Hidden options:
|
With Gram input the fold-based nmfkc.cv is not available
(it must split the columns of ). Select as in the
NMF-LAB paper: fix the landmarks, build one Gram object per candidate
on the training columns (pass the landmark matrix
g$landmarks as V so all candidates share it), fit with
nmfkc(), and score validation columns with
predict(fit, newA = nmfkc.kernel(g$landmarks, U.val, beta = beta)).
An object of class "nmfkc.gram": a list with
S matrix
(non-negative, positive semi-definite).
G0 matrix .
N, D
Number of samples and of covariates
(D = M, the number of landmarks).
type, signed
"nystrom" and FALSE
(non-negative features; accepted by nmfkc and
nmfkc.signed).
landmarksThe landmark matrix
. Pass it to nmfkc.kernel as its first
argument to build covariates for new data.
beta, kernel
Bandwidth and kernel used.
landmarks.method, sample.size
How the landmarks
were chosen ("supplied" when V was a matrix).
block.size, n.blocks
Blocking actually used.
A.blockA function A.block(idx) returning the
kernel block .
nmfkc uses it to rebuild block by
block after the fit. It closes over U and V, so the
object keeps U alive (no copy is made).
rownamesRow names of (colnames(V) if
any, else NULL).
This function is experimental. The interface may change in future versions.
Williams, C. K. I., & Seeger, M. (2001). Using the Nystr\"om method to speed up kernel machines. Advances in NIPS, 13.
Zhang, K., Tsang, I. W., & Kwok, J. T. (2008). Improved Nystr\"om low-rank approximation and error analysis. ICML.
Arthur, D., & Vassilvitskii, S. (2007). k-means++: the advantages of careful seeding. SODA.
Satoh, K. (2026). Applying non-negative matrix factorization with covariates to label matrix for classification. Japanese Journal of Statistics and Data Science. doi:10.1007/s42081-026-00349-x
nmfkc, nmfkc.kernel,
nmfkc.kernel.beta.nearest.med,
nmfkc.signed.rff.gram (the signed, Random-Fourier
counterpart for nmfkc.signed)
## Iris, 3 classes: Nystroem covariates with 12 k-means++ landmarks. ## The Gram route gives the same fit as the explicit M x N matrix. data(iris) set.seed(1) idx <- sample(nrow(iris), 100) mn <- colMeans(iris[idx, 1:4]); sc <- apply(iris[idx, 1:4], 2, sd) U.train <- t(scale(iris[idx, 1:4], center = mn, scale = sc)) # 4 x 100 U.test <- t(scale(iris[-idx, 1:4], center = mn, scale = sc)) # 4 x 50 levs <- levels(iris$Species) Y.train <- sapply(iris$Species[idx], function(s) as.integer(levs == s)) rownames(Y.train) <- levs ## 12 landmarks by k-means++, bandwidth by the nearest-landmark median ## heuristic, S and G0 accumulated over blocks of 25 columns g <- nmfkc.kernel.gram(Y.train, U.train, V = 12, block.size = 25, seed = 1) g res.g <- nmfkc(Y.train, A = g, rank = 3, verbose = FALSE) ## Same landmarks and beta, explicit 12 x 100 matrix: same fit A <- nmfkc.kernel(g$landmarks, U.train, beta = g$beta) res.m <- nmfkc(Y.train, A = A, rank = 3, verbose = FALSE) all.equal(res.g$C, res.m$C) ## Test-set prediction: kernel values against the stored landmarks A.test <- nmfkc.kernel(g$landmarks, U.test, beta = g$beta) pred <- predict(res.g, newA = A.test, type = "class") mean(pred == as.character(iris$Species[-idx]))## Iris, 3 classes: Nystroem covariates with 12 k-means++ landmarks. ## The Gram route gives the same fit as the explicit M x N matrix. data(iris) set.seed(1) idx <- sample(nrow(iris), 100) mn <- colMeans(iris[idx, 1:4]); sc <- apply(iris[idx, 1:4], 2, sd) U.train <- t(scale(iris[idx, 1:4], center = mn, scale = sc)) # 4 x 100 U.test <- t(scale(iris[-idx, 1:4], center = mn, scale = sc)) # 4 x 50 levs <- levels(iris$Species) Y.train <- sapply(iris$Species[idx], function(s) as.integer(levs == s)) rownames(Y.train) <- levs ## 12 landmarks by k-means++, bandwidth by the nearest-landmark median ## heuristic, S and G0 accumulated over blocks of 25 columns g <- nmfkc.kernel.gram(Y.train, U.train, V = 12, block.size = 25, seed = 1) g res.g <- nmfkc(Y.train, A = g, rank = 3, verbose = FALSE) ## Same landmarks and beta, explicit 12 x 100 matrix: same fit A <- nmfkc.kernel(g$landmarks, U.train, beta = g$beta) res.m <- nmfkc(Y.train, A = A, rank = 3, verbose = FALSE) all.equal(res.g$C, res.m$C) ## Test-set prediction: kernel values against the stored landmarks A.test <- nmfkc.kernel(g$landmarks, U.test, beta = g$beta) pred <- predict(res.g, newA = A.test, type = "class") mean(pred == as.character(iris$Species[-idx]))
Single entry point for symmetric NMF of network data with correct
multiplicative updates. Three model types are supported via type:
tri (type="tri", default):
with (both non-negative;
symmetric by design). Uses Frobenius-full bilateral gradient.
bi (type="bi"):
( fixed to ),
cube-root damping (He et al. 2011).
signed (type="signed"):
with and signed
. Preserves the soft-clustering interpretation
of while allowing negative off-diagonals of
(inter-cluster repulsion).
Non-negative adjacency matrix assumption. All three types
assume (a non-negative adjacency/affinity matrix).
The qualifier “signed” in type = "signed" refers to the
middle coefficient , not to itself. The underlying
Ding, Li & Jordan (2010) sign-splitting updates require
to guarantee monotone descent; supplying a signed triggers an
error. For a signed data matrix, see nmfkc.signed.
nmfkc.net( Y, rank = 2, type = c("tri", "bi", "signed"), epsilon = 1e-04, maxit = 5000, verbose = FALSE, ... )nmfkc.net( Y, rank = 2, type = c("tri", "bi", "signed"), epsilon = 1e-04, maxit = 5000, verbose = FALSE, ... )
Y |
Symmetric (N x N) non-negative adjacency matrix.
|
rank |
Integer Q. |
type |
|
epsilon, maxit, verbose
|
Standard. Note: the convergence test
uses the unpenalized squared error, not the penalized objective the
updates minimize, so with a large |
... |
Hidden options:
The anchor options of Multi-start recommendation. For |
Object of class c("nmfkc.net.<type>", "nmfkc.net", "nmfkc").
For type = "signed" the return also carries $Cp, $Cn.
This function is experimental. The interface may change in future versions; details are to be described in an upcoming paper.
Creates a Graphviz DOT script that visualizes the two-layer structure
of a symmetric NMF model ().
The resulting diagram displays:
outer nodes: the original network nodes (rows/columns of Y),
inner nodes: the latent basis/group nodes (columns of X),
directed edges from basis to node: membership weights from X,
undirected edges between basis nodes: inter-group interactions from C matrix (for tri-symmetric models only).
nmfkc.net.DOT( result, threshold = 0.01, sig.level = 0.1, weight_scale = 5, weight_scale_xy = 1, weight_scale_xx = weight_scale, rankdir = "TB", fill = TRUE, hide.isolated = TRUE, Y.label = NULL, X.label = NULL, Y.title = "Nodes (Y)", X.title = "Basis (X)", show.theta = NULL, signed = inherits(result, "nmfkc.net.signed"), cluster.box = c("none", "normal", "faint", "invisible"), layout = c("neato", "fdp", "twopi", "circo", "dot"), X.color = NULL, Y.cluster = c("soft", "hard") )nmfkc.net.DOT( result, threshold = 0.01, sig.level = 0.1, weight_scale = 5, weight_scale_xy = 1, weight_scale_xx = weight_scale, rankdir = "TB", fill = TRUE, hide.isolated = TRUE, Y.label = NULL, X.label = NULL, Y.title = "Nodes (Y)", X.title = "Basis (X)", show.theta = NULL, signed = inherits(result, "nmfkc.net.signed"), cluster.box = c("none", "normal", "faint", "invisible"), layout = c("neato", "fdp", "twopi", "circo", "dot"), X.color = NULL, Y.cluster = c("soft", "hard") )
result |
A list returned by |
threshold |
Minimum coefficient value to display an edge. It is read
on the membership scale for X edges (
The default |
sig.level |
Significance level for filtering C edges (if inference
results are present). Set to |
weight_scale |
Base scaling factor for edge widths. |
weight_scale_xy |
Scaling factor for X edges (basis -> node). |
weight_scale_xx |
Scaling factor for C edges (basis <-> basis). |
rankdir |
Graphviz rank direction ( |
fill |
Logical; whether nodes are filled with color. |
hide.isolated |
Logical; if TRUE, omit outer nodes with no X edge
above |
Y.label |
Character vector of labels for outer nodes. |
X.label |
Character vector of labels for basis nodes. |
Y.title |
Cluster title for outer nodes. |
X.title |
Cluster title for basis nodes. |
show.theta |
Logical or NULL. Whether to draw C edges between basis nodes. NULL = auto-detect (TRUE for tri, FALSE for bi). |
signed |
Logical. If |
cluster.box |
Style of cluster box: |
layout |
Graphviz layout engine, in recommended order:
|
X.color |
Color palette for basis nodes (length Q). |
Y.cluster |
Coloring mode for outer nodes: |
A character string of class c("nmfkc.net.DOT", "nmfkc.DOT")
representing a valid Graphviz DOT script.
Use plot() to render (requires DiagrammeR).
This function is experimental. The interface may change in future versions; details are to be described in an upcoming paper.
nmfkc.net.inference, nmfkc.DOT,
plot.nmfkc.DOT
library(nmfkc) Y <- matrix(c(0,1,1,0,0,0, 1,0,1,0,0,0, 1,1,0,1,0,0, 0,0,1,0,1,1, 0,0,0,1,0,1, 0,0,0,1,1,0), 6, 6) res <- nmfkc.net(Y, rank = 2, type = "tri", nstart = 20) dot <- nmfkc.net.DOT(res) plot(dot)library(nmfkc) Y <- matrix(c(0,1,1,0,0,0, 1,0,1,0,0,0, 1,1,0,1,0,0, 0,0,1,0,1,1, 0,0,0,1,0,1, 0,0,0,1,1,0), 6, 6) res <- nmfkc.net(Y, rank = 2, type = "tri", nstart = 20) dot <- nmfkc.net.DOT(res) plot(dot)
k-fold CV with folds taken over the upper triangle of the
symmetric (mirrored to the lower triangle) to prevent information
leakage through symmetry. A single entry point covers all three symmetric
NMF variants; type selects the fitting function:
nmfkc.net.ecv(Y, rank = 1:3, type = c("tri", "bi", "signed"), ...)nmfkc.net.ecv(Y, rank = 1:3, type = c("tri", "bi", "signed"), ...)
Y |
Symmetric N x N non-negative matrix. |
rank |
Integer vector of ranks to evaluate. Default |
type |
Model type: |
... |
Passed to the underlying fitter; also accepts:
|
A list with objfunc, sigma, objfunc.fold,
folds, Q.grid, type, and
r.squared.cvHeld-out per rank, pooled over
the folds: over
all held-out upper-triangle entries (weighted when Y.weights
is given). Unlike the r.squared of a fit (a squared
correlation), it can be negative.
predWith pred = TRUE:
(ranks) array of held-out predictions, mirrored to the lower triangle;
NA where an entry was never held out.
diag.excludeTRUE, present only when used.
seeds, objfunc.rep, objfunc.sd,
r.squared.cv.rep, r.squared.cv.sd, folds.rep
Present only with two or more seeds: the per-seed values
(ranks x seeds), their standard deviations, and each seed's folds.
objfunc.fold, folds and pred are those of the
first seed.
This function is experimental. The interface may change in future versions; details are to be described in an upcoming paper.
Performs statistical inference on the parameter matrix of a
symmetric NMF model ().
This is a wrapper around nmfkc.inference that automatically
sets the covariate matrix to , which is the defining
property of symmetric NMF.
nmfkc.net.inference(object, Y, wild.bootstrap = TRUE, ...)nmfkc.net.inference(object, Y, wild.bootstrap = TRUE, ...)
object |
A fitted model returned by |
Y |
The original symmetric matrix (P x P). |
wild.bootstrap |
Logical; if TRUE (default), perform wild bootstrap inference for the C matrix. |
... |
Additional arguments passed to |
The object augmented with inference fields
(same as nmfkc.inference), with class
c("nmfkc.net.inference", "nmfkc.inference", "nmf.inference", "nmfkc", "nmf").
This function is experimental. The interface may change in future versions; details are to be described in an upcoming paper.
nmfkc.inference, nmfkc.net.DOT,
summary.nmfkc.net.inference
library(nmfkc) Y <- matrix(c(0,1,1,0,0,0, 1,0,1,0,0,0, 1,1,0,1,0,0, 0,0,1,0,1,1, 0,0,0,1,0,1, 0,0,0,1,1,0), 6, 6) res <- nmfkc.net(Y, rank = 2, type = "tri", nstart = 20) res_inf <- nmfkc.net.inference(res, Y) summary(res_inf)library(nmfkc) Y <- matrix(c(0,1,1,0,0,0, 1,0,1,0,0,0, 1,1,0,1,0,0, 0,0,1,0,1,1, 0,0,0,1,0,1, 0,0,0,1,1,0), 6, 6) res <- nmfkc.net(Y, rank = 2, type = "tri", nstart = 20) res_inf <- nmfkc.net.inference(res, Y) summary(res_inf)
Fits nmfkc.net across a range of ranks and reports the
three rank-selection criteria – r.squared, the effective rank
(utilization), and the element-wise CV error sigma.ecv – with
the same concise diagnostics plot as nmfkc.rank.
nmfkc.net.rank( Y, rank = 1:5, type = c("tri", "bi", "signed"), detail = c("full", "fast"), plot = TRUE, ... )nmfkc.net.rank( Y, rank = 1:5, type = c("tri", "bi", "signed"), detail = c("full", "fast"), plot = TRUE, ... )
Y |
Symmetric (network) observation matrix. |
rank |
Integer vector of ranks to evaluate. |
type |
One of |
detail |
|
plot |
Logical; draw the diagnostics plot (default |
... |
Passed on to |
A list with rank.best (ECV minimum, or the R-squared
elbow under detail = "fast") and criteria (data
frame: rank, effective.rank, effective.rank.ratio,
r.squared, sigma.ecv).
Roy, O., & Vetterli, M. (2007). The effective rank: A measure of
effective dimensionality. Proc. EUSIPCO, 606–610.
(effective.rank)
Wold, S. (1978). Cross-validatory estimation of the number of
components in factor and principal components models.
Technometrics, 20(4), 397–405. (sigma.ecv)
nmfkc.net, nmfkc.net.ecv,
nmfkc.rank
Y <- matrix(c(0,1,1,0,0,0, 1,0,1,0,0,0, 1,1,0,1,0,0, 0,0,1,0,1,1, 0,0,0,1,0,1, 0,0,0,1,1,0), 6, 6) nmfkc.net.rank(Y, rank = 1:3, type = "tri", nstart = 5, nfolds = 3)Y <- matrix(c(0,1,1,0,0,0, 1,0,1,0,0,0, 1,1,0,1,0,0, 0,0,1,0,1,1, 0,0,0,1,0,1, 0,0,0,1,1,0), 6, 6) nmfkc.net.rank(Y, rank = 1:3, type = "tri", nstart = 5, nfolds = 3)
nmfkc.normalize rescales the values of a matrix to lie between 0 and 1
using the column-wise minimum and maximum values of a reference matrix.
nmfkc.normalize(x, ref = x)nmfkc.normalize(x, ref = x)
x |
A numeric matrix (or vector) to be normalized. |
ref |
A reference matrix from which the column-wise minima and maxima are taken.
Default is |
A matrix of the same dimensions as x, with each column rescaled to the range.
# Example. x <- nmfkc.normalize(iris[,-5]) apply(x,2,range)# Example. x <- nmfkc.normalize(iris[,-5]) apply(x,2,range)
nmfkc.rank provides diagnostic criteria for selecting the rank ()
in NMF with kernel covariates. Three rank-selection measures are computed
(R-squared, the effective rank, and the element-wise CV error), and results
can be visualized in a plot. Sample-clustering quality (silhouette / CPCC /
dist.cor) is no longer part of rank selection; use nmf.cluster.criteria
on a fitted model for those.
By default (save.time = FALSE), this function also computes the
Element-wise Cross-Validation error (Wold's CV Sigma) using nmfkc.ecv.
The plot explicitly marks the "BEST" rank based on two criteria:
Elbow Method (Red): Based on the curvature of the R-squared values (always computed if Q > 2).
Min RMSE (Blue): Based on the minimum Element-wise CV Sigma (only if detail="full").
nmfkc.rank(Y, A = NULL, rank = 1:2, detail = "full", plot = TRUE, data, ...)nmfkc.rank(Y, A = NULL, rank = 1:2, detail = "full", plot = TRUE, data, ...)
Y |
Observation matrix, or a formula (see |
A |
Covariate matrix. If |
rank |
A vector of candidate ranks to be evaluated. |
detail |
|
plot |
Logical. If |
data |
A data frame (required when |
... |
Additional arguments passed to
|
A list containing:
rank.best |
The estimated optimal rank. Prioritizes ECV minimum if available, otherwise R-squared Elbow. |
criteria |
A data frame containing diagnostic metrics for each rank.
The |
Roy, O., & Vetterli, M. (2007). The effective rank: A measure of
effective dimensionality. Proc. 15th European Signal Processing
Conf. (EUSIPCO), 606–610. (effective.rank)
Wold, S. (1978). Cross-validatory estimation of the number of
components in factor and principal components models.
Technometrics, 20(4), 397–405.
doi:10.1080/00401706.1978.10489693 (sigma.ecv)
nmfkc,
nmfkc.ecv (element-wise CV, used internally),
nmfkc.bicv (block bi-cross-validation),
nmfkc.consensus (stability) and
nmfkc.ard (Bayesian ARD) for alternative rank criteria.
# Example. Y <- t(iris[,-5]) # Full run (default) nmfkc.rank(Y, rank=1:4) # Fast run (skip ECV) nmfkc.rank(Y, rank=1:4, detail="fast")# Example. Y <- t(iris[,-5]) # Full run (default) nmfkc.rank(Y, rank=1:4) # Fast run (skip ECV) nmfkc.rank(Y, rank=1:4, detail="fast")
This function generates a side-by-side plot of three heatmaps: the original observation matrix Y, the fitted matrix XB (from NMF), and the residual matrix E (Y - XB). This visualization aids in diagnosing whether the chosen rank Q is adequate by assessing if the residual matrix E appears to be random noise.
The axis labels (X-axis: Samples, Y-axis: Features) are integrated into the main title of each plot to maximize the plot area, reflecting the compact layout settings.
nmfkc.residual.plot( Y, result, fitted.palette = (grDevices::colorRampPalette(c("white", "orange", "red")))(256), residual.palette = (grDevices::colorRampPalette(c("blue", "white", "red")))(256), ... )nmfkc.residual.plot( Y, result, fitted.palette = (grDevices::colorRampPalette(c("white", "orange", "red")))(256), residual.palette = (grDevices::colorRampPalette(c("blue", "white", "red")))(256), ... )
Y |
The original observation matrix (P x N). |
result |
The result object returned by the nmfkc function. |
fitted.palette |
A vector of colors for Y and XB heatmaps. Defaults to white-orange-red.
For backward compatibility, |
residual.palette |
A vector of colors for the residuals heatmap. Defaults to blue-white-red.
For backward compatibility, |
... |
Additional graphical parameters passed to the internal image calls. |
NULL. The function generates a plot.
Y <- t(iris[1:30, 1:4]) result <- nmfkc(Y, rank = 2) nmfkc.residual.plot(Y, result)Y <- t(iris[1:30, 1:4]) result <- nmfkc(Y, rank = 2) nmfkc.residual.plot(Y, result)
Chooses the Gaussian-kernel bandwidth (and optionally the
number of random features ) for Random Fourier Feature
covariates by cross-validation, in the same way that
nmfkc.kernel.beta.cv does for exact kernel covariates.
For every candidate the random features are regenerated on the columns
of U and the model is cross-validated column-wise:
type = "signed": cosine features from
nmfkc.signed.rff (signed), fitted with
nmfkc.signed through nmfkc.signed.cv;
type = "positive": positive random features from
nmfkc.rff.positive (non-negative), fitted with
nmfkc through nmfkc.cv.
The candidate with the smallest cross-validated prediction error is
returned. Since the same random map cannot be shared across folds
without leaking (the map is data-independent, so it can: only the fit
is cross-validated), one draw per candidate is used and its
seed is fixed so the comparison across is paired.
For large set sample.size: the selection then runs on a
random subsample of the columns (default: all columns), after which the
chosen is used with nmfkc.signed.rff.gram /
nmfkc.rff.positive.gram on the full data.
nmfkc.rff.beta.cv( Y, rank = 2, U, beta = NULL, D = 500L, type = c("signed", "positive"), plot = TRUE, ... )nmfkc.rff.beta.cv( Y, rank = 2, U, beta = NULL, D = 500L, type = c("signed", "positive"), plot = TRUE, ... )
Y |
Response matrix |
rank |
Number of bases |
U |
Input matrix |
beta |
Numeric vector of candidate bandwidths. Default
|
D |
Integer vector of candidate numbers of random features. A
single value (default |
type |
|
plot |
Logical. If |
... |
Passed to the CV function (
|
A list with
beta, D
The selected bandwidth and feature count.
objfuncMatrix of CV objective values,
length(beta) rows by length(D) columns.
beta.candidates, D.candidates
The grids.
type, sample.size, n.used
Settings used.
This function is experimental. The interface may change in future versions.
nmfkc.kernel.beta.cv (exact kernel; also works for
Nystr\"om covariates with U = landmarks, V = data),
nmfkc.signed.rff.gram, nmfkc.rff.positive.gram
data(iris) set.seed(1) idx <- sample(nrow(iris), 100) U <- t(scale(iris[idx, 1:4])) levs <- levels(iris$Species) Y <- sapply(iris$Species[idx], function(s) as.integer(levs == s)); rownames(Y) <- levs ## signed RFF: choose beta over the median-heuristic grid at D = 50 sel <- nmfkc.rff.beta.cv(Y, rank = 3, U, D = 50, plot = FALSE) sel$beta g <- nmfkc.signed.rff.gram(Y, U, beta = sel$beta, D = 50, seed = 1) fit <- nmfkc.signed(Y, A = g, rank = 3, verbose = FALSE) ## positive RFF: select over a (beta, D) grid sel2 <- nmfkc.rff.beta.cv(Y, rank = 3, U, D = c(100, 400), type = "positive", plot = FALSE, verbose = FALSE) sel2$objfuncdata(iris) set.seed(1) idx <- sample(nrow(iris), 100) U <- t(scale(iris[idx, 1:4])) levs <- levels(iris$Species) Y <- sapply(iris$Species[idx], function(s) as.integer(levs == s)); rownames(Y) <- levs ## signed RFF: choose beta over the median-heuristic grid at D = 50 sel <- nmfkc.rff.beta.cv(Y, rank = 3, U, D = 50, plot = FALSE) sel$beta g <- nmfkc.signed.rff.gram(Y, U, beta = sel$beta, D = 50, seed = 1) fit <- nmfkc.signed(Y, A = g, rank = 3, verbose = FALSE) ## positive RFF: select over a (beta, D) grid sel2 <- nmfkc.rff.beta.cv(Y, rank = 3, U, D = c(100, 400), type = "positive", plot = FALSE, verbose = FALSE) sel2$objfunc
Generates and applies the
non-negative feature map
for which
:
the Gaussian kernel with bandwidth beta, times a constant. The
matrix therefore satisfies
and is non-negative, so it
can be used as the covariate matrix of nmfkc directly –
unlike the cosine Random Fourier Features of
nmfkc.signed.rff, which take both signs and need
nmfkc.signed. This is the "positive random features"
construction of Choromanski et al. (2021), stated there for the softmax
kernel, transported to the Gaussian kernel by rescaling.
The constant is minus the largest exponent
observed
on the training U (stored in pars), so every training
feature is at most and exp() cannot overflow
(points far from the others underflow to zero at large , which
only means that bandwidth fits poorly). It scales the kernel estimate by
, which absorbs, so the fitted model and its
predictions do not depend on it. Note that depends on the
realized : the stabilized features are an unbiased estimator of
the randomly rescaled kernel , not of itself; the
unbiasedness statement above is for the unstabilized map.
Inputs should nevertheless be centred and scaled: the
estimator's variance grows with , so
positive features need a larger D than cosine features for the
same accuracy.
nmfkc.rff.positive(U, beta = NULL, D = ceiling(ncol(U)/2), seed = NULL, ...)nmfkc.rff.positive(U, beta = NULL, D = ceiling(ncol(U)/2), seed = NULL, ...)
U |
A |
beta |
Positive scalar. Gaussian kernel bandwidth parameter
|
D |
Integer. Number of random features (default
|
seed |
Optional integer passed to |
... |
Hidden options:
|
A list with Z, the non-negative feature
matrix, and pars = list(omega, D, beta, kappa, hyperbolic, type =
"positive"). Pass Z to nmfkc (or
nmfkc.signed) as A, and pars back to this
function to build features for new data.
This function is experimental. The interface may change in future versions.
Choromanski, K., Likhosherstov, V., Dohan, D., Song, X., Gane, A., Sarlos, T., Hawkins, P., Davis, J., Mohiuddin, A., Kaiser, L., Belanger, D., Colwell, L., & Weller, A. (2021). Rethinking attention with Performers. ICLR. arXiv:2009.14794.
Rahimi, A., & Recht, B. (2007). Random features for large-scale kernel machines. Advances in NIPS, 20.
nmfkc.rff.positive.gram (large ),
nmfkc.signed.rff (signed cosine features),
nmfkc.kernel, nmfkc
set.seed(1) U <- scale(matrix(stats::rnorm(3 * 40), 3, 40)) # centred inputs pr <- nmfkc.rff.positive(U, beta = 0.5, D = 20000, seed = 1) all(pr$Z >= 0) ## Phi' Phi / exp(2 kappa) approximates the Gaussian kernel K <- nmfkc.kernel(U, beta = 0.5) Kp <- crossprod(pr$Z) / exp(2 * pr$pars$kappa) max(abs(Kp - K)) ## The same random map on new data U.new <- matrix(stats::rnorm(3 * 5), 3, 5) dim(nmfkc.rff.positive(U.new, pars = pr$pars)$Z) # 20000 x 5set.seed(1) U <- scale(matrix(stats::rnorm(3 * 40), 3, 40)) # centred inputs pr <- nmfkc.rff.positive(U, beta = 0.5, D = 20000, seed = 1) all(pr$Z >= 0) ## Phi' Phi / exp(2 kappa) approximates the Gaussian kernel K <- nmfkc.kernel(U, beta = 0.5) Kp <- crossprod(pr$Z) / exp(2 * pr$pars$kappa) max(abs(Kp - K)) ## The same random map on new data U.new <- matrix(stats::rnorm(3 * 5), 3, 5) dim(nmfkc.rff.positive(U.new, pars = pr$pars)$Z) # 20000 x 5
The nmfkc.rff.positive counterpart of
nmfkc.signed.rff.gram: walks the columns of U in
blocks, generates the non-negative feature block
with the same random map, accumulates
() and
(), and discards the block, so the
feature matrix never exists in memory. Because the features are
non-negative, the returned "nmfkc.gram" object is accepted by
nmfkc (as well as by nmfkc.signed), and the
NMF-LAB membership interpretation of is retained.
nmfkc.rff.positive.gram( Y, U, beta = NULL, D = NULL, seed = NULL, block.size = NULL, ... )nmfkc.rff.positive.gram( Y, U, beta = NULL, D = NULL, seed = NULL, block.size = NULL, ... )
Y |
Non-negative |
U |
A |
beta |
Positive scalar. Gaussian kernel bandwidth. Ignored when
|
D |
Integer. Number of random features. Required (no |
seed |
Optional integer seed for |
block.size |
Integer. Columns of |
... |
Hidden options |
An object of class "nmfkc.gram" with S, G0,
N, D, type = "prf", signed = FALSE,
pars, beta, block.size, n.blocks,
A.block (block generator) and rownames = NULL. Features
for new data: nmfkc.rff.positive(U.new, pars = g$pars)$Z.
This function is experimental. The interface may change in future versions.
Choromanski, K. et al. (2021). Rethinking attention with Performers. ICLR. arXiv:2009.14794.
nmfkc.rff.positive, nmfkc,
nmfkc.kernel.gram (Nystr\"om), nmfkc.signed.rff.gram
data(iris) set.seed(1) idx <- sample(nrow(iris), 100) mn <- colMeans(iris[idx, 1:4]); sc <- apply(iris[idx, 1:4], 2, sd) U.train <- t(scale(iris[idx, 1:4], center = mn, scale = sc)) U.test <- t(scale(iris[-idx, 1:4], center = mn, scale = sc)) levs <- levels(iris$Species) Y.train <- sapply(iris$Species[idx], function(s) as.integer(levs == s)) rownames(Y.train) <- levs beta <- nmfkc.kernel.beta.nearest.med(U.train)$beta g <- nmfkc.rff.positive.gram(Y.train, U.train, beta = beta, D = 200, seed = 1, block.size = 25) g res <- nmfkc(Y.train, A = g, rank = 3, verbose = FALSE) # standard NMF-LAB Z.test <- nmfkc.rff.positive(U.test, pars = g$pars)$Z pred <- predict(res, newA = Z.test, type = "class") mean(pred == as.character(iris$Species[-idx]))data(iris) set.seed(1) idx <- sample(nrow(iris), 100) mn <- colMeans(iris[idx, 1:4]); sc <- apply(iris[idx, 1:4], 2, sd) U.train <- t(scale(iris[idx, 1:4], center = mn, scale = sc)) U.test <- t(scale(iris[-idx, 1:4], center = mn, scale = sc)) levs <- levels(iris$Species) Y.train <- sapply(iris$Species[idx], function(s) as.integer(levs == s)) rownames(Y.train) <- levs beta <- nmfkc.kernel.beta.nearest.med(U.train)$beta g <- nmfkc.rff.positive.gram(Y.train, U.train, beta = beta, D = 200, seed = 1, block.size = 25) g res <- nmfkc(Y.train, A = g, rank = 3, verbose = FALSE) # standard NMF-LAB Z.test <- nmfkc.rff.positive(U.test, pars = g$pars)$Z pred <- predict(res, newA = Z.test, type = "class") mean(pred == as.character(iris$Species[-idx]))
Solves
where the covariate matrix and the coefficient matrix
may be signed. Internally
and with
(sign-splitting trick, Ding et al. 2010), and the problem is solved
by a Direct Multiplicative Update algorithm whose iteration cost is
, independent of .
Only is structurally constrained to be non-negative (Semi-NMF
sense of Ding, Li, & Jordan 2010). In particular, may
contain negative entries, in which case the response is fit in the
least-squares sense without any non-negativity requirement on .
When (so ), the result reduces to
nmfkc(Y, A, rank) with Euclidean loss, up to reordering.
nmfkc.signed( Y, A, rank = NULL, epsilon = 1e-04, maxit = 5000, verbose = TRUE, ... )nmfkc.signed( Y, A, rank = NULL, epsilon = 1e-04, maxit = 5000, verbose = TRUE, ... )
Y |
Real-valued |
A |
Real-valued Large |
rank |
Integer. Number of latent components |
epsilon |
Relative convergence tolerance on the objective
(default |
maxit |
Maximum number of iterations (default |
verbose |
Logical. Print dimensions at start (default |
... |
Additional arguments:
|
An object of class c("nmfkc.signed", "nmfkc") with
X: basis matrix (non-negative,
column-normalized according to X.restriction).
Cp, Cn: non-negative parts of
, so that .
C: (= ), signed.
B: , (signed).
objfunc.iter: objective values per iteration.
objfunc: final objective.
r.squared: (Pearson; in ).
r.squared.uncentered: uncentered (baseline = zero matrix).
r.squared.centered: row-mean centered .
mae: mean absolute error.
iter: number of iterations performed.
converged: logical; whether the relative change of the
objective fell below epsilon before maxit.
epsilon.iter: relative change of the objective at the last
step (the quantity compared with epsilon).
objfunc.increases: number of steps at which the objective
rose. Every restriction on offer is a gauge fix and the unweighted
multiplicative update is monotone, so this should be zero; a positive
count means something outside the multiplicative form is pushing the
iterate back, and such a fit can oscillate and exhaust maxit.
It is what identified the two offenders removed in 0.9.8.
runtime: elapsed seconds.
Y.signed: logical; whether contained negative
entries during fitting.
pars: RFF generating parameters, if supplied.
restarts: present only when nstart.signed > 1; a data frame
with one row per start (seed, objfunc, iter,
converged, epsilon.iter, objfunc.increases).
A start whose objective is larger than the best one
has reached a different local minimum; such a start typically also stops
after far fewer iterations.
call: the matched call.
This function is experimental. The interface may change in future versions.
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE TPAMI, 32(1), 45–55.
Rahimi, A., & Recht, B. (2007). Random features for large-scale kernel machines. Advances in NIPS, 20.
set.seed(1) ## Example 1: signed A (e.g., hand-built RFF features), non-negative Y ## Build simple signed features Z = sqrt(2/D) * cos(omega^T U + b): U <- matrix(stats::rnorm(5 * 40), 5, 40) # raw input D <- 20 # feature dim omega <- matrix(stats::rnorm(5 * D), 5, D) # random freqs b <- stats::runif(D, 0, 2 * pi) # phase Z <- sqrt(2 / D) * cos(t(omega) %*% U + matrix(b, D, 40)) # D x 40, signed Y <- matrix(abs(stats::rnorm(8 * 40)), 8, 40) res1 <- nmfkc.signed(Y, A = Z, rank = 3, maxit = 200) ## Example 2: signed Y (regression) Y2 <- matrix(stats::rnorm(8 * 40), 8, 40) # signed response res2 <- nmfkc.signed(Y2, A = Z, rank = 3, maxit = 200, warm.start = FALSE)set.seed(1) ## Example 1: signed A (e.g., hand-built RFF features), non-negative Y ## Build simple signed features Z = sqrt(2/D) * cos(omega^T U + b): U <- matrix(stats::rnorm(5 * 40), 5, 40) # raw input D <- 20 # feature dim omega <- matrix(stats::rnorm(5 * D), 5, D) # random freqs b <- stats::runif(D, 0, 2 * pi) # phase Z <- sqrt(2 / D) * cos(t(omega) %*% U + matrix(b, D, 40)) # D x 40, signed Y <- matrix(abs(stats::rnorm(8 * 40)), 8, 40) res1 <- nmfkc.signed(Y, A = Z, rank = 3, maxit = 200) ## Example 2: signed Y (regression) Y2 <- matrix(stats::rnorm(8 * 40), 8, 40) # signed response res2 <- nmfkc.signed(Y2, A = Z, rank = 3, maxit = 200, warm.start = FALSE)
Column-wise k-fold CV by held-out samples: for each
fold, the model is fit on the training columns and evaluated on the
held-out columns by solving for new-sample coefficients via a
weighted refit with fixed.
nmfkc.signed.cv(Y, A, rank = 2, ...)nmfkc.signed.cv(Y, A, rank = 2, ...)
Y |
Real-valued |
A |
Real-valued |
rank |
Integer |
... |
Passed to |
A list with objfunc (mean squared prediction error),
sigma (RMSE), objfunc.block (per-fold MSE vector),
block (integer fold assignment of length ). Field
names match nmfkc.cv.
This function is experimental.
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1), 45–55.
nmfkc.signed, nmfkc.signed.ecv
Element-wise k-fold CV (Wold's CV): held-out elements
are masked via Y.weights = 0 during fitting, and the RMSE on
those elements is reported. Loops over candidate rank values.
nmfkc.signed.ecv(Y, A, rank = 1:3, ...)nmfkc.signed.ecv(Y, A, rank = 1:3, ...)
Y |
Real-valued |
A |
Real-valued |
rank |
Integer vector of candidate ranks (default |
... |
Passed to |
A list with objfunc (MSE per rank), sigma
(RMSE), objfunc.fold (per-fold per-rank), folds,
Q.grid.
This function is experimental.
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1), 45–55.
Fits nmfkc.signed across a range of ranks and reports
r.squared, the effective rank, and the element-wise CV error
sigma.ecv, with the same concise plot as
nmfkc.rank.
nmfkc.signed.rank( Y, A, rank = 1:5, detail = c("full", "fast"), plot = TRUE, ... )nmfkc.signed.rank( Y, A, rank = 1:5, detail = c("full", "fast"), plot = TRUE, ... )
Y |
Observation matrix (may contain negative entries). |
A |
Covariate matrix (may be signed). |
rank |
Integer vector of ranks to evaluate. |
detail |
|
plot |
Logical; draw the diagnostics plot (default |
... |
Passed on to |
A list with rank.best and criteria
(rank, effective.rank, effective.rank.ratio,
r.squared, sigma.ecv).
Roy, O., & Vetterli, M. (2007). The effective rank: A measure of
effective dimensionality. Proc. EUSIPCO, 606–610.
(effective.rank)
Wold, S. (1978). Cross-validatory estimation of the number of
components in factor and principal components models.
Technometrics, 20(4), 397–405. (sigma.ecv)
Generates RFF random parameters
,
(Rahimi & Recht, 2007) and
applies the RFF transform
to each column of U, yielding a sign-unrestricted
feature matrix such that
, the Gaussian kernel matrix with bandwidth
beta.
The return value is a list with the feature matrix Z and
the generating parameters pars = list(omega, b, D, beta) so
that the same random map can be re-applied to new data (by passing
pars back) and nmfkc.signed can record the
parameters for downstream summary.
nmfkc.signed.rff(U, beta = NULL, D = ceiling(ncol(U)/2), seed = NULL, ...)nmfkc.signed.rff(U, beta = NULL, D = ceiling(ncol(U)/2), seed = NULL, ...)
U |
A |
beta |
Positive scalar. Gaussian kernel bandwidth parameter.
Can be obtained via |
D |
Integer. Number of random features. Defaults to
|
seed |
Optional integer passed to |
... |
Hidden option |
A list with two elements:
ZA sign-unrestricted numeric
matrix. Pass this to nmfkc.signed as its
A argument.
parsA list list(omega, b, D, beta). Pass
this to nmfkc.signed via its pars
argument (for summary display) and to subsequent
nmfkc.signed.rff() calls (to reuse the same random map
on new data).
This function is experimental. The interface may change in future versions; details are to be described in an upcoming paper.
nmfkc.signed,
nmfkc.kernel.beta.nearest.med
## Iris 3-class classification with RFF + direct MU (Ding-Li-Jordan) data(iris) set.seed(1) idx <- sample(nrow(iris), 100) # 100 training, 50 test ## Scale features using TRAINING mean/sd; transpose to p x N layout mn <- colMeans(iris[idx, 1:4]) sc <- apply(iris[idx, 1:4], 2, sd) U.train <- t(scale(iris[idx, 1:4], center = mn, scale = sc)) # 4 x 100 U.test <- t(scale(iris[-idx, 1:4], center = mn, scale = sc)) # 4 x 50 ## One-hot encode training labels as a Q_obs x N target matrix levs <- levels(iris$Species) Y.train <- sapply(iris$Species[idx], function(s) as.integer(levs == s)) rownames(Y.train) <- levs # 3 x 100 lab.train <- iris$Species[idx] lab.test <- iris$Species[-idx] ## Beta candidates from nearest-neighbour median heuristic beta_info <- nmfkc.kernel.beta.nearest.med(U.train) betas <- beta_info$beta_candidates ## CV over beta candidates: for each beta, generate RFF, fit, and ## evaluate column-wise CV-MSE on training data cv_mse <- numeric(length(betas)) for (i in seq_along(betas)) { rff_i <- nmfkc.signed.rff(U.train, beta = betas[i], D = 50, seed = 1) cv_i <- nmfkc.signed.cv(Y.train, A = rff_i$Z, rank = 3, seed = 123) cv_mse[i] <- cv_i$objfunc } beta_best <- betas[which.min(cv_mse)] ## Generate signed RFF features with the best beta rff.train <- nmfkc.signed.rff(U.train, beta = beta_best, D = 50, seed = 1) rff.test <- nmfkc.signed.rff(U.test, pars = rff.train$pars) ## Fit on training data only res <- nmfkc.signed(Y.train, A = rff.train$Z, rank = 3, pars = rff.train$pars, verbose = FALSE) ## Predict on training and test separately pred.train <- predict(res, newA = rff.train$Z, type = "class") pred.test <- predict(res, newA = rff.test$Z, type = "class") mean(pred.train == as.character(lab.train)) mean(pred.test == as.character(lab.test))## Iris 3-class classification with RFF + direct MU (Ding-Li-Jordan) data(iris) set.seed(1) idx <- sample(nrow(iris), 100) # 100 training, 50 test ## Scale features using TRAINING mean/sd; transpose to p x N layout mn <- colMeans(iris[idx, 1:4]) sc <- apply(iris[idx, 1:4], 2, sd) U.train <- t(scale(iris[idx, 1:4], center = mn, scale = sc)) # 4 x 100 U.test <- t(scale(iris[-idx, 1:4], center = mn, scale = sc)) # 4 x 50 ## One-hot encode training labels as a Q_obs x N target matrix levs <- levels(iris$Species) Y.train <- sapply(iris$Species[idx], function(s) as.integer(levs == s)) rownames(Y.train) <- levs # 3 x 100 lab.train <- iris$Species[idx] lab.test <- iris$Species[-idx] ## Beta candidates from nearest-neighbour median heuristic beta_info <- nmfkc.kernel.beta.nearest.med(U.train) betas <- beta_info$beta_candidates ## CV over beta candidates: for each beta, generate RFF, fit, and ## evaluate column-wise CV-MSE on training data cv_mse <- numeric(length(betas)) for (i in seq_along(betas)) { rff_i <- nmfkc.signed.rff(U.train, beta = betas[i], D = 50, seed = 1) cv_i <- nmfkc.signed.cv(Y.train, A = rff_i$Z, rank = 3, seed = 123) cv_mse[i] <- cv_i$objfunc } beta_best <- betas[which.min(cv_mse)] ## Generate signed RFF features with the best beta rff.train <- nmfkc.signed.rff(U.train, beta = beta_best, D = 50, seed = 1) rff.test <- nmfkc.signed.rff(U.test, pars = rff.train$pars) ## Fit on training data only res <- nmfkc.signed(Y.train, A = rff.train$Z, rank = 3, pars = rff.train$pars, verbose = FALSE) ## Predict on training and test separately pred.train <- predict(res, newA = rff.train$Z, type = "class") pred.test <- predict(res, newA = rff.test$Z, type = "class") mean(pred.train == as.character(lab.train)) mean(pred.test == as.character(lab.test))
Builds the two summary matrices that nmfkc.signed needs
from a Random Fourier Feature (RFF) covariate matrix
() without ever forming itself:
The input U is processed in column blocks; for each block the
features are generated with
nmfkc.signed.rff, added into and , and
discarded. Peak memory is therefore
instead of , which is what makes the RFF route usable when
is in the hundreds of thousands and exceeds the input
dimension by a large factor (the matrix is
times the size of the data itself).
The returned object is passed to nmfkc.signed as its
A argument. The multiplicative updates are unchanged, since
they only ever use and ; the fit equals the one
obtained from the explicit matrix nmfkc.signed.rff(U, ...)$Z
with warm.start = FALSE up to floating-point summation order.
nmfkc.signed.rff.gram( Y, U, beta = NULL, D = NULL, seed = NULL, block.size = NULL, ... )nmfkc.signed.rff.gram( Y, U, beta = NULL, D = NULL, seed = NULL, block.size = NULL, ... )
Y |
Real-valued |
U |
A |
beta |
Positive scalar. Gaussian kernel bandwidth
|
D |
Integer. Number of random features. Required (there is no
|
seed |
Optional integer passed to |
block.size |
Integer. Number of columns of |
... |
Hidden option |
An object of class "nmfkc.gram": a list with
S matrix (signed,
positive semi-definite).
G0 matrix
(signed).
N, D
Number of samples and of features.
parsThe RFF parameters list(omega, b, D, beta).
Pass them to nmfkc.signed.rff to generate features
for new data, e.g. for predict.nmfkc.signed.
block.size, n.blocks
Blocking actually used.
A.blockA function A.block(idx) returning the
feature block for the columns
idx of U. nmfkc.signed uses it to
rebuild block by block after the fit. It closes over
U and pars, so the object keeps U alive
(no copy is made).
rownamesAlways NULL (RFF features are
unnamed); present so that nmfkc.signed() can treat the
object like a matrix.
With Gram input the fold-based nmfkc.signed.cv is not
available (it must split the columns of ). Select
as in the NMF-LAB paper: build one Gram object per candidate on the
training columns, fit, and score the validation columns with
predict() on features from nmfkc.signed.rff(U.val,
pars = g$pars).
This function is experimental. The interface may change in future versions.
Rahimi, A., & Recht, B. (2007). Random features for large-scale kernel machines. Advances in NIPS, 20.
Satoh, K. (2026). Applying non-negative matrix factorization with covariates to label matrix for classification. Japanese Journal of Statistics and Data Science. doi:10.1007/s42081-026-00349-x
nmfkc.signed, nmfkc.signed.rff,
predict.nmfkc.signed,
nmfkc.kernel.beta.nearest.med
## Iris, 3 classes: the Gram route gives the same fit as the explicit ## feature matrix, without ever holding the D x N matrix. data(iris) set.seed(1) idx <- sample(nrow(iris), 100) mn <- colMeans(iris[idx, 1:4]); sc <- apply(iris[idx, 1:4], 2, sd) U.train <- t(scale(iris[idx, 1:4], center = mn, scale = sc)) # 4 x 100 U.test <- t(scale(iris[-idx, 1:4], center = mn, scale = sc)) # 4 x 50 levs <- levels(iris$Species) Y.train <- sapply(iris$Species[idx], function(s) as.integer(levs == s)) rownames(Y.train) <- levs beta <- nmfkc.kernel.beta.nearest.med(U.train)$beta ## Accumulate S, G0 over blocks of 25 columns (4 blocks here) g <- nmfkc.signed.rff.gram(Y.train, U.train, beta = beta, D = 50, seed = 1, block.size = 25) g res.g <- nmfkc.signed(Y.train, A = g, rank = 3, verbose = FALSE) ## Same random map, explicit matrix, no warm-start: same fit Z <- nmfkc.signed.rff(U.train, pars = g$pars)$Z # 50 x 100 res.m <- nmfkc.signed(Y.train, A = Z, rank = 3, warm.start = FALSE, verbose = FALSE) all.equal(res.g$C, res.m$C) ## Test-set prediction: regenerate features with the stored pars Z.test <- nmfkc.signed.rff(U.test, pars = g$pars)$Z pred <- predict(res.g, newA = Z.test, type = "class") mean(pred == as.character(iris$Species[-idx]))## Iris, 3 classes: the Gram route gives the same fit as the explicit ## feature matrix, without ever holding the D x N matrix. data(iris) set.seed(1) idx <- sample(nrow(iris), 100) mn <- colMeans(iris[idx, 1:4]); sc <- apply(iris[idx, 1:4], 2, sd) U.train <- t(scale(iris[idx, 1:4], center = mn, scale = sc)) # 4 x 100 U.test <- t(scale(iris[-idx, 1:4], center = mn, scale = sc)) # 4 x 50 levs <- levels(iris$Species) Y.train <- sapply(iris$Species[idx], function(s) as.integer(levs == s)) rownames(Y.train) <- levs beta <- nmfkc.kernel.beta.nearest.med(U.train)$beta ## Accumulate S, G0 over blocks of 25 columns (4 blocks here) g <- nmfkc.signed.rff.gram(Y.train, U.train, beta = beta, D = 50, seed = 1, block.size = 25) g res.g <- nmfkc.signed(Y.train, A = g, rank = 3, verbose = FALSE) ## Same random map, explicit matrix, no warm-start: same fit Z <- nmfkc.signed.rff(U.train, pars = g$pars)$Z # 50 x 100 res.m <- nmfkc.signed(Y.train, A = Z, rank = 3, warm.start = FALSE, verbose = FALSE) all.equal(res.g$C, res.m$C) ## Test-set prediction: regenerate features with the stored pars Z.test <- nmfkc.signed.rff(U.test, pars = g$pars)$Z pred <- predict(res.g, newA = Z.test, type = "class") mean(pred == as.character(iris$Species[-idx]))
Estimates the NMF-RE model
where () is a non-negative observation matrix,
() is a non-negative basis matrix learned from the data,
() is a non-negative coefficient matrix capturing
systematic covariate effects on latent scores,
() is a covariate matrix, and
() is a random effects matrix capturing
unit-specific deviations in the latent score space.
NMF-RE can be viewed as a mixed-effects latent-variable model defined on a
reconstruction (mean) structure. The non-negativity constraint on
induces sparse, parts-based loadings, achieving measurement-side variable
selection without an explicit sparsity penalty. Inference on
provides covariate-side variable selection by identifying which covariates
significantly affect which components.
Estimation is an outer–inner ECM: an inner loop (fixed ) runs a
ridge/BLUP update for , a complete-EM semi-NMF update for , and
a fixed-effect update for ; an outer loop runs the EM M-steps for
. The variance components are estimated from the data;
the effective degrees of freedom are reported only as a diagnostic.
nmfre() performs optimization only. Hypothesis tests and
standard errors for are obtained separately with
nmfre.inference (sandwich SE + wild bootstrap), mirroring the
nmfkc / nmfkc.inference split.
nmfre( Y, A = NULL, rank = 2, C.signed = TRUE, epsilon = 1e-05, maxit = 5000, ... )nmfre( Y, A = NULL, rank = 2, C.signed = TRUE, epsilon = 1e-05, maxit = 5000, ... )
Y |
Observation matrix (P x N), non-negative. |
A |
Covariate matrix (K x N). Default is a row of ones (intercept only). |
rank |
Integer. Rank of the basis matrix |
C.signed |
Logical. Whether the fixed-effect coefficients |
epsilon |
Convergence tolerance for relative change in objective (default 1e-5). |
maxit |
Maximum number of iterations. Default |
... |
Additional arguments for initialization, variance control, dfU control, optimization, and inference settings.
|
A list of class "nmfre" with components.
The model is .
Core matrices
XBasis matrix (), columns normalized to sum to 1.
X.probRow-wise soft-clustering probabilities from the
non-negative (each row normalized to sum to 1), as in nmfkc.
X.clusterHard-clustering label for each row of
(argmax over X.prob).
CCoefficient matrix ().
URandom effects matrix ().
Variance components
sigma2Residual variance .
tau2Random effect variance .
lambdaRidge penalty .
Convergence diagnostics
convergedLogical. Whether the algorithm converged.
stop.reasonCharacter string describing why iteration stopped.
iterNumber of iterations performed.
maxitMaximum iterations setting used.
epsilonConvergence tolerance used.
objfuncFinal objective function value
.
rel.change.finalFinal relative change in objective.
objfunc.iterNumeric vector of the fixed-
penalized objective
per iteration. Monotone within an inner loop but not across
outer iterations (the penalty jumps when is updated).
rss.traceNumeric vector of per iteration.
nll.traceNumeric vector of the marginal negative
log-likelihood per iteration
(random effects integrated out). This is the ECM-monotone quantity and
is what plot.nmfre displays.
Effective degrees of freedom (dfU) diagnostics
dfUFinal effective degrees of freedom
,
where are eigenvalues of .
dfU.capUpper bound imposed on .
dfU.cap.rateRate used to compute the cap.
lambda.enforcedFinal enforced to satisfy the cap.
dfU.hit.capLogical. Whether the cap was binding.
dfU.hit.iterIteration at which the cap first bound.
dfU.frac, fraction of maximum df.
dfU.cap.frac.
Fitted matrices
BFixed-effect scores ().
B.probColumn-normalized probabilities from
.
B.blupBLUP scores ().
B.blup.posNon-negative BLUP scores
().
B.blup.probColumn-normalized probabilities from
.
XBFitted values from fixed effects
().
XB.blupFitted values including random effects
().
Fit statistics
r.squaredPearson
(BLUP prediction).
r.squared.uncenteredUncentered
(BLUP;
baseline = zero matrix).
r.squared.centeredRow-mean centered
(BLUP; baseline = per-row mean).
r.squared.fixedPearson
(fixed-only prediction).
r.squared.fixed.uncentered, r.squared.fixed.centered
Uncentered
and centered for the fixed-only prediction.
ICCTrace-based Intraclass Correlation Coefficient.
In the NMF-RE model, the conditional covariance of the -th
observation column is
,
a matrix. Unlike a standard random intercept model
where the design matrix is a simple indicator (so the ICC
reduces to ), the basis matrix
plays the role of in a random slopes model,
making the variance contribution of depend on .
To obtain a scalar summary, we take the trace of each component:
This equals the average (over dimensions) proportion of
per-column variance attributable to the random effects.
Sign convention
C.signedLogical. Whether was sign-free (TRUE) or non-negative (FALSE).
Standard errors, z-values, p-values, and confidence intervals for
are not computed here; obtain them by passing the fit
to nmfre.inference.
Satoh, K. (2026). Wild Bootstrap Inference for Non-Negative Matrix Factorization with Random Effects. arXiv:2603.01468. https://arxiv.org/abs/2603.01468
nmfre.inference,
nmfkc.DOT, summary.nmfre
# Example 1. cars data Y <- matrix(cars$dist, nrow = 1) A <- rbind(intercept = 1, speed = cars$speed) res <- nmfre(Y, A, rank = 1, maxit = 5000) summary(res) # Example 2. Orthodont data (nlme) if (requireNamespace("nlme", quietly = TRUE)) { Y <- matrix(nlme::Orthodont$distance, 4, 27) male <- ifelse(nlme::Orthodont$Sex[seq(1, 108, 4)] == "Male", 1, 0) A <- rbind(intercept = 1, male = male) # Fit (sign-free Theta by default; variances estimated by EM/ECM) res <- nmfre(Y, A, rank = 1) summary(res) }# Example 1. cars data Y <- matrix(cars$dist, nrow = 1) A <- rbind(intercept = 1, speed = cars$speed) res <- nmfre(Y, A, rank = 1, maxit = 5000) summary(res) # Example 2. Orthodont data (nlme) if (requireNamespace("nlme", quietly = TRUE)) { Y <- matrix(nlme::Orthodont$distance, 4, 27) male <- ifelse(nlme::Orthodont$Sex[seq(1, 108, 4)] == "Male", 1, 0) A <- rbind(intercept = 1, male = male) # Fit (sign-free Theta by default; variances estimated by EM/ECM) res <- nmfre(Y, A, rank = 1) summary(res) }
Selects the basis rank for nmfre by Wold-style
element-wise (entry-holdout) cross-validation that exercises the
random effects. For each fold the held-out entries are hidden and filled by
iterative imputation: fit the NMF-RE model on the current matrix, replace the
held-out entries with the BLUP prediction , and repeat.
The score is the held-out prediction RMSE (sigma); the selected
rank minimizes it.
Because the held-out entries of a column are predicted using the random
effect fitted from that column's retained entries, this
evaluates the full NMF-RE model (including ) — unlike
nmfkc.ecv, which masks entries by zero weight and predicts from
the fixed-effect fit only. The two scores are therefore
not directly comparable.
nmfre.ecv(Y, A = NULL, rank = 1:3, C.signed = TRUE, ...)nmfre.ecv(Y, A = NULL, rank = 1:3, C.signed = TRUE, ...)
Y |
Observation matrix (P x N), non-negative. |
A |
Covariate matrix (K x N). Default is a row of ones (intercept only). |
rank |
Integer vector of ranks |
C.signed |
Logical. Sign convention for |
... |
Additional arguments:
|
A list of class "nmfre.ecv" with components:
rankThe ranks evaluated.
sigmaNamed numeric vector of held-out RMSE per rank
(same field name as nmfkc.ecv).
bestThe rank minimizing sigma.
nfolds, rounds, C.signed
Settings used.
nmfre, nmfre.inference, nmfkc.ecv
if (requireNamespace("nlme", quietly = TRUE)) { Y <- matrix(nlme::Orthodont$distance, 4, 27) male <- ifelse(nlme::Orthodont$Sex[seq(1, 108, 4)] == "Male", 1, 0) A <- rbind(intercept = 1, male = male) cv <- nmfre.ecv(Y, A, rank = 1:3) cv$best plot(cv) }if (requireNamespace("nlme", quietly = TRUE)) { Y <- matrix(nlme::Orthodont$distance, 4, 27) male <- ifelse(nlme::Orthodont$Sex[seq(1, 108, 4)] == "Male", 1, 0) A <- rbind(intercept = 1, male = male) cv <- nmfre.ecv(Y, A, rank = 1:3) cv$best plot(cv) }
nmfre.inference performs statistical inference on the coefficient
matrix () from a fitted nmfre model,
conditional on the estimated basis matrix and random
effects .
Under the working model ,
inference is conducted via sandwich covariance estimation and
one-step wild bootstrap with non-negative projection.
The result is compatible with nmfkc.DOT for visualization
(pass the result directly as x with type = "YXA").
nmfre.inference(object, Y, A = NULL, wild.bootstrap = TRUE, ...)nmfre.inference(object, Y, A = NULL, wild.bootstrap = TRUE, ...)
object |
An object of class |
Y |
Observation matrix (P x N). Must match the data used in
|
A |
Covariate matrix (K x N). Default is |
wild.bootstrap |
Logical. If |
... |
Additional arguments:
|
The input object with additional inference components:
sigma2.used |
Estimated |
C.vec.cov |
Full covariance matrix for |
C.se |
Sandwich standard errors for |
C.se.boot |
Bootstrap standard errors for |
C.ci.lower |
Lower CI bounds for |
C.ci.upper |
Upper CI bounds for |
coefficients |
Data frame with Basis, Covariate, Estimate, SE, BSE, z_value, p_value, CI_low, CI_high. |
C.p.side |
P-value type used. |
identified.approx |
Whether |
identification |
The parts of the diagnostic, plus |
Satoh, K. (2026). Wild Bootstrap Inference for Non-Negative Matrix Factorization with Random Effects. arXiv:2603.01468. https://arxiv.org/abs/2603.01468
nmfre, nmfkc.DOT,
summary.nmfre
Y <- matrix(cars$dist, nrow = 1) A <- rbind(intercept = 1, speed = cars$speed) res <- nmfre(Y, A, rank = 1, wild.bootstrap = FALSE) res2 <- nmfre.inference(res, Y, A) res2$coefficientsY <- matrix(cars$dist, nrow = 1) A <- rbind(intercept = 1, speed = cars$speed) res <- nmfre(Y, A, rank = 1, wild.bootstrap = FALSE) res2 <- nmfre.inference(res, Y, A) res2$coefficients
Line plot of silhouette, CPCC, and dist.cor
against the result index, for an object from
nmf.cluster.criteria. X-axis ticks default to each
result's $rank (overridable via the names argument of
nmf.cluster.criteria).
## S3 method for class 'nmf.cluster.criteria' plot( x, main = "Clustering quality across rank", xlab = "rank (Q)", ylab = "criterion", lwd = 2, ... )## S3 method for class 'nmf.cluster.criteria' plot( x, main = "Clustering quality across rank", xlab = "rank (Q)", ylab = "criterion", lwd = 2, ... )
x |
An object of class |
main |
Plot title. |
xlab, ylab
|
Axis labels. |
lwd |
Line width. |
... |
Further arguments passed to the initial
|
x, invisibly.
Draws the alluvial / Sankey-style cluster-flow diagram for an object
created by nmf.cluster.flow.
## S3 method for class 'nmf.cluster.flow' plot( x, col = NULL, lwd = 1, xlab = "rank (Q)", ylab = "individuals", main = "Cluster flow across rank", ... )## S3 method for class 'nmf.cluster.flow' plot( x, col = NULL, lwd = 1, xlab = "rank (Q)", ylab = "individuals", main = "Cluster flow across rank", ... )
x |
An object of class |
col |
Optional colour vector indexed by reference
cluster id ( |
lwd |
Line width of the flow segments. |
xlab, ylab
|
Axis labels. |
main |
Plot title. |
... |
Further arguments passed to the initial
|
x, invisibly.
type = "convergence" (default) plots the EM objective
() over iterations. type = "adjusted.scores" draws the
cluster scatterplot, coloured by the hard cluster assignment (or by
group).
## S3 method for class 'nmf.gmm' plot( x, type = c("convergence", "adjusted.scores", "scores"), group = NULL, ... )## S3 method for class 'nmf.gmm' plot( x, type = c("convergence", "adjusted.scores", "scores"), group = NULL, ... )
x |
An object of class |
type |
|
group |
Optional length-N factor/vector to colour the scatterplot by
(e.g. known labels). Default |
... |
Additional graphical parameters passed to |
The scatterplot shows the adjusted scores, not the scores. What is
drawn is , the least-squares scores with the
covariate effect removed — which is the space the mixture actually acts in,
since the model puts the covariate effect in the component means. The
distinction matters: on the crabs example the raw scores
separate the four true groups no better with a size covariate than without
(0.635 versus 0.648, between/within on the first two principal components),
while the adjusted space separates them at 4.628. The visible separation is
produced by the adjustment, so a panel labelled "scores" would claim
something the picture does not show. type = "scores" is accepted as
an alias.
The projection depends on the rank: a rank-2 fit is drawn in its own two
coordinates (rotating it through prcomp would trade two interpretable
part axes for arbitrary ones), a higher rank is projected onto its first two
principal components, and a rank-1 fit becomes a strip plot with
deterministic stacking.
Invisible NULL.
Part of the experimental NMF-GMM family: see nmf.gmm.
The interface may change in future versions — argument names, defaults
and the contents of the returned object are not yet stable.
nmf.gmm, nmf.gmm.inference, nmf.gmm.select
set.seed(1) Y <- matrix(abs(rnorm(20 * 80)) + 1, 20, 80); A <- rbind(1, rnorm(80)) fit <- nmf.gmm(Y, A, rank = 2, K = 3) plot(fit, type = "adjusted.scores")set.seed(1) Y <- matrix(abs(rnorm(20 * 80)) + 1, 20, 80); A <- rbind(1, rnorm(80)) fit <- nmf.gmm(Y, A, rank = 2, K = 3) plot(fit, type = "adjusted.scores")
Draws BIC and ICL against , with the selected marked. Every
other selector in the package (nmfkc.rank,
nmfkc.ard, nmfkc.consensus,
nmfkc.bicv) has a plot method; without this one
plot(nmf.gmm.select(...)) fell through to plot.default and
errored with "'x' is a list, but does not have components 'x' and 'y'".
## S3 method for class 'nmf.gmm.select' plot(x, ...)## S3 method for class 'nmf.gmm.select' plot(x, ...)
x |
An object of class |
... |
Passed to the underlying plot call. |
x, invisibly.
Part of the experimental NMF-GMM family: see nmf.gmm.
The interface may change in future versions — argument names, defaults
and the contents of the returned object are not yet stable.
Draws the concise three-criterion rank-selection plot for an object
returned by nmfkc.rank (or nmfkc.net.rank,
nmfkc.signed.rank, nmfae.rank, nmfae.signed.rank):
r.squared (red) and the broken-stick-corrected effective-rank
index effective.rank.index (green, legend label
eff.rank; not the raw effective.rank.ratio)
on the left axis, and the
cross-validation error sigma.ecv (blue) on the right axis, each
with points, rank-number labels and a highlighted best marker.
## S3 method for class 'nmf.rank' plot( x, main = NULL, xlab = "Rank (Q)", ylab = "R-squared / eff.rank (0-1)", lwd = 3, ... )## S3 method for class 'nmf.rank' plot( x, main = NULL, xlab = "Rank (Q)", ylab = "R-squared / eff.rank (0-1)", lwd = 3, ... )
x |
An object of class |
main |
Plot title (defaults to the title stored in |
xlab, ylab
|
Axis labels. |
lwd |
Line width of the criterion curves. |
... |
Further arguments passed to the initial
|
x, invisibly.
plot.nmfae displays the convergence trajectory of the objective function
across iterations. The title shows the achieved plot.nmfae displays the convergence trajectory of the objective function
across iterations. The title shows the achieved .
## S3 method for class 'nmfae' plot(x, ...)## S3 method for class 'nmfae' plot(x, ...)
x |
An object of class |
... |
Additional graphical parameters passed to |
Invisible NULL. Called for its side effect (plot).
set.seed(1) Y <- matrix(runif(20), nrow = 4) res <- nmf.rrr(Y, rank1 = 2) plot(res)set.seed(1) Y <- matrix(runif(20), nrow = 4) res <- nmf.rrr(Y, rank1 = 2) plot(res)
nmfkc
plot.nmfkc produces a diagnostic plot for the return value of
nmfkc, showing the objective function across iterations.
## S3 method for class 'nmfkc' plot(x, ...)## S3 method for class 'nmfkc' plot(x, ...)
x |
An object of class |
... |
Additional arguments passed to the base |
Called for its side effect (a plot). Returns NULL invisibly.
Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) plot(result)Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) plot(result)
Draws the model's own dynamics — what the fit does — rather than
the realized series. All of it lives in the latent
transition matrices returned by
nmfkc.ar.latent, so every panel stays readable at a size the
observed space cannot: a 47-region, 7-lag, rank-4 fit has
observed impulse responses but only latent ones.
## S3 method for class 'nmfkc.ar.latent' plot( x, type = c("roots", "graph", "irf", "lag", "phase"), lag = NULL, horizon = 16L, ... )## S3 method for class 'nmfkc.ar.latent' plot( x, type = c("roots", "graph", "irf", "lag", "phase"), lag = NULL, horizon = 16L, ... )
x |
An object from |
type |
Which panel to draw; see Details. |
lag |
Which lag to use for |
horizon |
Number of periods for |
... |
Passed to the underlying plot call. |
The available panels are
"roots"(default) The eigenvalues of the latent
companion matrix on the unit circle. Complex roots are drawn in a
different colour and their implied period is reported: a fit can sit near
the boundary either because it cycles or because it is seasonal, and the
spectral radius alone does not say which.
"graph"The transition graph of : self-loops are
persistence, arrows are spillover, widths are proportional
to the entries. Unlike a bipartite picture this graph
is closed, so an asymmetric feedback loop is visible as such.
"irf"Latent impulse responses , a
panel. Since the responses never change sign, so the
shape is a decay or a hump and its half-life can be read off.
"lag" as a function of : how many
periods a shock takes to arrive. With this is a single bar
plot of the lag profile.
"phase"Phase portrait of
with the fixed point and the eigen-directions. Requires
and . This is where a monotone relaxation and a true
cycle look different; two coefficient series plotted against time do not
distinguish them.
invisible(NULL); called for the plot. "irf" and
"lag" set up their own multi-panel layout and restore the previous
one on exit.
nmfkc.ar.latent,
nmfkc.ar.latent.inference,
plot.nmfkc.ar.stationarity, nmfkc.ar.DOT
set.seed(1) Y <- matrix(abs(rnorm(4 * 60)) + 1, 4, 60) ar <- nmfkc.ar(Y, degree = 3) fit <- nmfkc(ar$Y, ar$A, rank = 2, verbose = FALSE) lat <- nmfkc.ar.latent(fit) plot(lat) # roots on the unit circle plot(lat, type = "graph") # latent transition graph plot(lat, type = "irf") # latent impulse responsesset.seed(1) Y <- matrix(abs(rnorm(4 * 60)) + 1, 4, 60) ar <- nmfkc.ar(Y, degree = 3) fit <- nmfkc(ar$Y, ar$A, rank = 2, verbose = FALSE) lat <- nmfkc.ar.latent(fit) plot(lat) # roots on the unit circle plot(lat, type = "graph") # latent transition graph plot(lat, type = "irf") # latent impulse responses
Draws the Perron-Frobenius bracket against the unit boundary, where are the column
sums of . It shows at a glance whether the cheap
sufficient condition settles stationarity or an eigenvalue computation is
needed. For the model's dynamics see plot.nmfkc.ar.latent.
## S3 method for class 'nmfkc.ar.stationarity' plot(x, ...)## S3 method for class 'nmfkc.ar.stationarity' plot(x, ...)
x |
An object from |
... |
Passed to the underlying plot call. |
invisible(NULL); called for the plot.
nmfkc.ar.stationarity,
plot.nmfkc.ar.latent
set.seed(1) Y <- matrix(abs(rnorm(4 * 60)) + 1, 4, 60) ar <- nmfkc.ar(Y, degree = 2) fit <- nmfkc(ar$Y, ar$A, rank = 2, verbose = FALSE) plot(nmfkc.ar.stationarity(fit))set.seed(1) Y <- matrix(abs(rnorm(4 * 60)) + 1, 4, 60) ar <- nmfkc.ar(Y, degree = 2) fit <- nmfkc(ar$Y, ar$A, rank = 2, verbose = FALSE) plot(nmfkc.ar.stationarity(fit))
Bar plot of the per-component relevance (descending), with a line at the pruning threshold; bars above it (the estimated rank) are highlighted.
## S3 method for class 'nmfkc.ard' plot(x, main = NULL, ...)## S3 method for class 'nmfkc.ard' plot(x, main = NULL, ...)
x |
An object of class |
main |
Plot title. |
... |
Passed to |
x, invisibly.
Two views of a nmfkc.consensus result:
type = "criteria" (default): the stability curves
cophenetic (blue) and dispersion (red) against rank.
No "best" marker is drawn – a stable rank tends to be a coarse one,
so the curves are shown for inspection rather than as an optimum.
type = "heatmap": the consensus matrix of each requested
rank, reordered by average-linkage hierarchical clustering of
(blue = 0 / different cluster, red = 1 / same
cluster). Requires keep.consensus = TRUE at compute time.
## S3 method for class 'nmfkc.consensus' plot( x, type = c("criteria", "heatmap"), rank = NULL, mfrow = NULL, mar = NULL, col = grDevices::hcl.colors(50, "Blue-Red"), main = NULL, ... )## S3 method for class 'nmfkc.consensus' plot( x, type = c("criteria", "heatmap"), rank = NULL, mfrow = NULL, mar = NULL, col = grDevices::hcl.colors(50, "Blue-Red"), main = NULL, ... )
x |
An object of class |
type |
|
rank |
For |
mfrow |
Panel layout for multiple heatmaps, as
|
mar |
Per-panel margins |
col |
Heatmap colour palette (length-50 blue-to-red by default). |
main |
Plot title for the |
... |
x, invisibly.
Renders a DOT graph string using DiagrammeR::grViz.
If the DiagrammeR package is not installed, prints the DOT source
to the console instead.
This method handles all DOT objects produced by the nmfkc package:
nmfkc.DOT, nmf.rrr.DOT, nmf.ffb.DOT,
and nmfkc.ar.DOT.
## S3 method for class 'nmfkc.DOT' plot(x, ...)## S3 method for class 'nmfkc.DOT' plot(x, ...)
x |
An object of class |
... |
Not used. |
Called for its side effect (rendering). Returns x invisibly.
nmfkc.DOT, nmf.rrr.DOT,
nmf.ffb.DOT, nmfkc.ar.DOT,
print.nmfkc.DOT
Plot method for nmfkc.signed (convergence)
## S3 method for class 'nmfkc.signed' plot(x, ...)## S3 method for class 'nmfkc.signed' plot(x, ...)
x |
An |
... |
Passed to |
Invisible x.
This function is experimental. The interface may change in future versions.
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1), 45–55.
nmfkc.signed, nmfkc, nmfkc.signed.cv
Plots the objective function value over iterations for nmfre and
nmf.ffb objects. (For nmfkc and nmfae, plot methods
are defined in their respective source files.)
## S3 method for class 'nmfre' plot(x, ...) ## S3 method for class 'nmf.ffb' plot(x, ..., which = c("penalized", "reconstruction", "both"))## S3 method for class 'nmfre' plot(x, ...) ## S3 method for class 'nmf.ffb' plot(x, ..., which = c("penalized", "reconstruction", "both"))
x |
A fitted model object. |
... |
Additional graphical arguments passed to |
which |
For |
Invisible NULL.
set.seed(1) Y <- matrix(runif(20), nrow = 4) A <- diag(5) res <- nmfre(Y, A, rank = 2) plot(res)set.seed(1) Y <- matrix(runif(20), nrow = 4) A <- diag(5) res <- nmfre(Y, A, rank = 2) plot(res)
Draws the element-wise CV score (held-out RMSE) against the rank, marking the minimizing rank as “Best (Min)”.
## S3 method for class 'nmfre.ecv' plot( x, main = "NMF-RE element-wise CV", xlab = "Rank (Q)", ylab = "sigma.ecv (held-out RMSE)", ... )## S3 method for class 'nmfre.ecv' plot( x, main = "NMF-RE element-wise CV", xlab = "Rank (Q)", ylab = "sigma.ecv (held-out RMSE)", ... )
x |
An object of class |
main |
Plot title. |
xlab, ylab
|
Axis labels. |
... |
Passed to |
The input object, invisibly.
Class assignments / responsibilities from an NMF-GMM fit
## S3 method for class 'nmf.gmm' predict(object, type = c("class", "responsibility"), ...)## S3 method for class 'nmf.gmm' predict(object, type = c("class", "responsibility"), ...)
object |
An object of class |
type |
|
... |
Ignored. |
A vector of labels or the responsibility matrix.
Part of the experimental NMF-GMM family: see nmf.gmm.
The interface may change in future versions — argument names, defaults
and the contents of the returned object are not yet stable.
nmf.gmm, nmf.gmm.inference, nmf.gmm.select
nmfkc
predict.nmfkc generates predictions from an object of class nmfkc,
either using the fitted covariates or a new covariate matrix.
When the model was fitted using a formula (Formula Mode), a newdata
data frame can be supplied instead of newA; the covariate matrix is
then constructed automatically from the stored formula metadata.
## S3 method for class 'nmfkc' predict( object, newA = NULL, newdata = NULL, type = c("response", "prob", "class"), ... )## S3 method for class 'nmfkc' predict( object, newA = NULL, newdata = NULL, type = c("response", "prob", "class"), ... )
object |
An object of class |
newA |
Optional. A new covariate matrix to be used for prediction. |
newdata |
Optional data frame. Only available when the model was fitted
using a formula. Covariate columns are extracted automatically using the
stored formula metadata. If both |
type |
Type of prediction to return. Options are "response" (fitted values matrix), "prob" (soft-clustering probabilities), or "class" (hard-clustering labels based on row names of X). |
... |
Further arguments passed to or from other methods. |
Depending on type: a numeric matrix ("response" or "prob")
or a character vector of class labels ("class").
# Prediction with newA Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) newA <- rbind(1, c(10, 20, 30)) predict(result, newA = newA)# Prediction with newA Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) newA <- rbind(1, c(10, 20, 30)) predict(result, newA = newA)
Computes
().
For type = "response" the raw prediction is returned
(possibly signed).
With signed covariates the scores can be
negative, so the NMF-LAB recipe of normalizing to a
membership vector does not apply. Since is column-stochastic,
still sums to
(close to one for one-hot targets) and is the least-squares estimate of
the class indicator; for type = "prob" it is mapped to the
probability simplex by the Euclidean projection
the closest probability vector in the same Frobenius geometry the fit
minimizes (Duchi et al. 2008; Wang & Carreira-Perpinan 2013;
per column by sorting). The former rule – clip
negatives to zero and renormalize – is kept as
prob.method = "clip". type = "class" is
, which both rules preserve.
## S3 method for class 'nmfkc.signed' predict(object, newA = NULL, type = c("response", "prob", "class"), ...)## S3 method for class 'nmfkc.signed' predict(object, newA = NULL, type = c("response", "prob", "class"), ...)
object |
A fitted |
newA |
Real-valued |
type |
Output: |
... |
Hidden option |
A numeric matrix ("response" or "prob") or a
character vector ("class").
This function is experimental. The interface may change in future versions.
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1), 45–55.
Duchi, J., Shalev-Shwartz, S., Singer, Y., & Chandra, T. (2008). Efficient projections onto the l1-ball for learning in high dimensions. Proceedings of the 25th International Conference on Machine Learning (ICML), 272–279.
Wang, W., & Carreira-Perpinan, M. A. (2013). Projection onto the probability simplex: an efficient algorithm with a simple proof, and an application. arXiv:1309.1541.
nmfkc.signed, nmfkc, nmfkc.signed.cv
Predicts the response from a fitted NMF-RE model. With newA supplied,
returns the fixed-effect prediction (new
units carry no estimated random effect , so only the population mean
is predicted). Without newA, returns the in-sample BLUP fit
(object$XB.blup), matching fitted().
## S3 method for class 'nmfre' predict(object, newA = NULL, ...)## S3 method for class 'nmfre' predict(object, newA = NULL, ...)
object |
An object of class |
newA |
Optional covariate matrix (K x M) for new units. If |
... |
Not used. |
A matrix of predicted (P x M) values.
Prints the call, the model shape and whether the fit converged, in the
spirit of print.lm: enough to see what was fitted and whether to
trust it, with everything else left to summary(). Registered on the
shared "nmf" class, so it covers every optimizer that returns one
(nmfkc, nmf.rrr / nmf.rrr,
nmfre, nmf.ffb and the signed variants);
nmf.gmm has its own methods.
Without it, printing a fit dumps the whole list — over 350 lines for a
problem, and unusable for a large one.
Unlike lm, the coefficient matrix is not shown: it is
with the number of covariates — the number of
samples when there is no covariate matrix — so it is routinely far too
large for a print method. Use coef for it.
## S3 method for class 'nmf' print(x, digits = max(3L, getOption("digits") - 3L), ...)## S3 method for class 'nmf' print(x, digits = max(3L, getOption("digits") - 3L), ...)
x |
A fitted model object inheriting class |
digits |
Number of significant digits. |
... |
Additional arguments (currently unused). |
x, invisibly.
Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) result coef(result) # the coefficient matrix CY <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) result coef(result) # the coefficient matrix C
Print method for nmf.cluster.criteria objects
## S3 method for class 'nmf.cluster.criteria' print(x, ...)## S3 method for class 'nmf.cluster.criteria' print(x, ...)
x |
An object of class |
... |
Passed to the criteria table's |
x, invisibly.
Prints a one-line header, the adjacent-result ARI, and the
cluster table (rows = individuals, columns =
results, entries = cluster number). Use plot for the
diagram.
## S3 method for class 'nmf.cluster.flow' print(x, ...)## S3 method for class 'nmf.cluster.flow' print(x, ...)
x |
An object of class |
... |
Passed to the table's |
x, invisibly.
Reports the observed likelihood ratio, its bootstrap -value, and the
two quantities that say whether the test means anything for these data: the
null false-selection rate and, for the default calibration, how often the
exclusion restriction moved.
## S3 method for class 'nmf.ffb.test' print(x, ...)## S3 method for class 'nmf.ffb.test' print(x, ...)
x |
An object from |
... |
Ignored. |
x, invisibly.
Print an NMF-GMM fit
## S3 method for class 'nmf.gmm' print(x, ...)## S3 method for class 'nmf.gmm' print(x, ...)
x |
An object of class |
... |
Ignored. |
x, invisibly.
Part of the experimental NMF-GMM family: see nmf.gmm.
The interface may change in future versions — argument names, defaults
and the contents of the returned object are not yet stable.
nmf.gmm, nmf.gmm.inference, nmf.gmm.select
Print method for nmf.gmm.select objects
## S3 method for class 'nmf.gmm.select' print(x, ...)## S3 method for class 'nmf.gmm.select' print(x, ...)
x |
An object of class |
... |
Ignored. |
x, invisibly.
Part of the experimental NMF-GMM family: see nmf.gmm.
The interface may change in future versions — argument names, defaults
and the contents of the returned object are not yet stable.
nmf.gmm, nmf.gmm.inference, nmf.gmm.select
Prints a summary of any NMF inference result object
("nmfkc.inference" or "nmfae.inference").
## S3 method for class 'nmf.inference' print(x, ...)## S3 method for class 'nmf.inference' print(x, ...)
x |
An object of class |
... |
Additional arguments passed to the corresponding
|
Called for its side effect (printing). Returns x invisibly.
nmfkc.inference, nmf.rrr.inference
Print method for rank-selection (nmf.rank) objects
## S3 method for class 'nmf.rank' print(x, ...)## S3 method for class 'nmf.rank' print(x, ...)
x |
An object of class |
... |
Passed to the criteria table's |
x, invisibly.
A two-line summary, matching print.nmfkc.ecv. Without it the object
printed as a raw list (33 lines of $objfunc / $folds).
## S3 method for class 'nmfae.ecv' print(x, ...)## S3 method for class 'nmfae.ecv' print(x, ...)
x |
An object of class |
... |
Ignored. |
x, invisibly.
Print method for nmfkc.ar.latent objects
## S3 method for class 'nmfkc.ar.latent' print(x, digits = 3, ...)## S3 method for class 'nmfkc.ar.latent' print(x, digits = 3, ...)
x |
An object of class |
digits |
Number of digits used when printing the matrices. |
... |
Ignored. |
x, invisibly.
Print method for nmfkc.ar.latent.inference objects
## S3 method for class 'nmfkc.ar.latent.inference' print(x, digits = 4, ...)## S3 method for class 'nmfkc.ar.latent.inference' print(x, digits = 4, ...)
x |
An object of class |
digits |
Number of digits used when printing. |
... |
Ignored. |
x, invisibly.
Print method for nmfkc.ar.stationarity objects
## S3 method for class 'nmfkc.ar.stationarity' print(x, digits = 3, ...)## S3 method for class 'nmfkc.ar.stationarity' print(x, digits = 3, ...)
x |
An object of class |
digits |
Ignored; kept for compatibility with the print generic. |
... |
Ignored. |
x, invisibly.
Print method for nmfkc.ard objects
## S3 method for class 'nmfkc.ard' print(x, ...)## S3 method for class 'nmfkc.ard' print(x, ...)
x |
An object of class |
... |
Unused. |
x, invisibly.
Print method for nmfkc.consensus objects
## S3 method for class 'nmfkc.consensus' print(x, ...)## S3 method for class 'nmfkc.consensus' print(x, ...)
x |
An object of class |
... |
Unused. |
x, invisibly.
Writes the Graphviz source out as text. Without it, printing a DOT object
fell to print.default, which showed the whole graph as one escaped
string ([1] "digraph NMF \{\n ...") — unreadable. Shared by every
*.DOT function through the common "nmfkc.DOT" class.
## S3 method for class 'nmfkc.DOT' print(x, ...)## S3 method for class 'nmfkc.DOT' print(x, ...)
x |
An object of class |
... |
Not used. |
x, invisibly.
Print method for nmfkc.gram
## S3 method for class 'nmfkc.gram' print(x, ...)## S3 method for class 'nmfkc.gram' print(x, ...)
x |
Object of class |
... |
Unused. |
Invisible x.
This function is experimental.
nmfkc.kernel.gram, nmfkc.signed.rff.gram
Print method for nmfre.ecv objects
## S3 method for class 'nmfre.ecv' print(x, ...)## S3 method for class 'nmfre.ecv' print(x, ...)
x |
An object of class |
... |
Not used. |
The input object, invisibly.
Prints the NMF-FFB model summary (dimensions, convergence, stability diagnostics, fit statistics, and inference results if available).
## S3 method for class 'summary.nmf.ffb' print(x, ...)## S3 method for class 'summary.nmf.ffb' print(x, ...)
x |
An object of class |
... |
Not used. |
Invisible x.
Print method for summary.nmf.gmm objects
## S3 method for class 'summary.nmf.gmm' print(x, digits = max(3L, getOption("digits") - 3L), ...)## S3 method for class 'summary.nmf.gmm' print(x, digits = max(3L, getOption("digits") - 3L), ...)
x |
A |
digits |
Minimum significant digits. |
... |
Ignored. |
x, invisibly.
Part of the experimental NMF-GMM family: see nmf.gmm.
The interface may change in future versions — argument names, defaults
and the contents of the returned object are not yet stable.
nmf.gmm, nmf.gmm.inference, nmf.gmm.select
summary.nmfkc objectsPrints a formatted summary of an nmfkc model fit.
## S3 method for class 'summary.nmfkc' print(x, digits = max(3L, getOption("digits") - 3L), ...)## S3 method for class 'summary.nmfkc' print(x, digits = max(3L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Minimum number of significant digits to be used. |
... |
Additional arguments (currently unused). |
Called for its side effect (printing). Returns x invisibly.
Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) print(summary(result))Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) print(summary(result))
Prints a formatted summary including the coefficients table.
## S3 method for class 'summary.nmfkc.inference' print( x, digits = max(3L, getOption("digits") - 3L), max.coef = 20, by = c("covariate", "basis"), ... )## S3 method for class 'summary.nmfkc.inference' print( x, digits = max(3L, getOption("digits") - 3L), max.coef = 20, by = c("covariate", "basis"), ... )
x |
An object of class |
digits |
Minimum number of significant digits. |
max.coef |
Maximum coefficient rows to display. Default is 20. |
by |
Character; grouping order of the coefficients table.
|
... |
Additional arguments (currently unused). |
Called for its side effect (printing). Returns x invisibly.
Print method for summary.nmfkc.net objects
## S3 method for class 'summary.nmfkc.net' print(x, digits = max(3L, getOption("digits") - 3L), ...)## S3 method for class 'summary.nmfkc.net' print(x, digits = max(3L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Minimum number of significant digits. |
... |
Unused. |
Invisible x.
Prints a formatted summary of a symmetric NMF model with inference results.
The coefficients table uses Basis.row and Basis.col labels
reflecting the symmetric structure .
## S3 method for class 'summary.nmfkc.net.inference' print( x, digits = max(3L, getOption("digits") - 3L), max.coef = 20, by = c("covariate", "basis"), ... )## S3 method for class 'summary.nmfkc.net.inference' print( x, digits = max(3L, getOption("digits") - 3L), max.coef = 20, by = c("covariate", "basis"), ... )
x |
An object of class |
digits |
Minimum number of significant digits. |
max.coef |
Largest number of coefficient rows to print (default 20).
The table has one row per basis pair, so it grows as |
by |
Character; grouping order of the coefficients table. Because the
symmetric model sets |
... |
Additional arguments (currently unused). |
Called for its side effect (printing). Returns x invisibly.
Print method for summary.nmfkc.net.signed objects
## S3 method for class 'summary.nmfkc.net.signed' print(x, digits = max(3L, getOption("digits") - 3L), ...)## S3 method for class 'summary.nmfkc.net.signed' print(x, digits = max(3L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Minimum number of significant digits. |
... |
Unused. |
Invisible x.
Print method for summary.nmfkc.signed
## S3 method for class 'summary.nmfkc.signed' print(x, digits = max(3L, getOption("digits") - 3L), ...)## S3 method for class 'summary.nmfkc.signed' print(x, digits = max(3L, getOption("digits") - 3L), ...)
x |
Object of class |
digits |
Number of significant digits. |
... |
Unused. |
Invisible x.
This function is experimental. The interface may change in future versions.
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1), 45–55.
nmfkc.signed, nmfkc, nmfkc.signed.cv
Returns the residual matrix from a fitted NMF model.
Requires the original observation matrix Y to be supplied.
For nmfre objects, residuals are computed from the BLUP
reconstruction () by default. Set
type = "fixed" to use fixed-effects only.
## S3 method for class 'nmf' residuals(object, Y, ...) ## S3 method for class 'nmfae' residuals(object, Y, ...) ## S3 method for class 'nmfre' residuals(object, Y, type = c("blup", "fixed"), ...) ## S3 method for class 'nmf.ffb' residuals(object, Y, ...)## S3 method for class 'nmf' residuals(object, Y, ...) ## S3 method for class 'nmfae' residuals(object, Y, ...) ## S3 method for class 'nmfre' residuals(object, Y, type = c("blup", "fixed"), ...) ## S3 method for class 'nmf.ffb' residuals(object, Y, ...)
object |
A fitted model object. |
Y |
The original observation matrix used for fitting. |
... |
Not used. |
type |
For |
The residual matrix.
nmfkc, nmf.rrr, nmfre,
nmf.ffb, fitted.nmf
Y <- matrix(runif(50), 5, 10) result <- nmfkc(Y, rank = 2) residuals(result, Y)Y <- matrix(runif(50), 5, 10) result <- nmfkc(Y, rank = 2) residuals(result, Y)
, matching residuals.nmf. Without this method
residuals() fell through to stats::residuals.default and
returned NULL, even though fitted() worked.
## S3 method for class 'nmf.gmm' residuals(object, Y, ...)## S3 method for class 'nmf.gmm' residuals(object, Y, ...)
object |
An object of class |
Y |
The observation matrix the model was fitted on. |
... |
Ignored. |
The residual matrix.
Part of the experimental NMF-GMM family: see nmf.gmm.
The interface may change in future versions — argument names, defaults
and the contents of the returned object are not yet stable.
Produces a formatted summary of a fitted NMF-FFB model, including
matrix dimensions, convergence, stability diagnostics, fit statistics,
and inference results (if available). For nmf.ffb(method =
"fiml") objects the fit statistics are the log-likelihoods, parameter
counts and BIC of the feed-forward null, the unpenalized feedback fit and
the BIC-selected model, the two likelihood-ratio statistics, and – after
nmf.ffb.inference – their bootstrap p-values and the
false-selection rate under the null.
## S3 method for class 'nmf.ffb' summary(object, ...)## S3 method for class 'nmf.ffb' summary(object, ...)
object |
An object of class |
... |
Not used. |
An object of class "summary.nmf.ffb" (the fitted model
tagged for printing); printed by print.summary.nmf.ffb.
Y <- t(iris[, -5]) Y1 <- Y[1:2, ]; Y2 <- Y[3:4, ] result <- nmf.ffb(Y1, Y2, rank = 2, maxit = 500) summary(result)Y <- t(iris[, -5]) Y1 <- Y[1:2, ]; Y2 <- Y[3:4, ] result <- nmf.ffb(Y1, Y2, rank = 2, maxit = 500) summary(result)
Summary of an NMF-GMM fit
## S3 method for class 'nmf.gmm' summary(object, ...)## S3 method for class 'nmf.gmm' summary(object, ...)
object |
An object of class |
... |
Ignored. |
An object of class "summary.nmf.gmm".
Part of the experimental NMF-GMM family: see nmf.gmm.
The interface may change in future versions — argument names, defaults
and the contents of the returned object are not yet stable.
nmf.gmm, nmf.gmm.inference, nmf.gmm.select
nmfkc
Produces a summary of an nmfkc object, including matrix dimensions,
runtime, fit statistics, and diagnostics.
## S3 method for class 'nmfkc' summary(object, ...)## S3 method for class 'nmfkc' summary(object, ...)
object |
An object of class |
... |
Additional arguments (currently unused). |
Separability measures how close the basis is to having an anchor
row for every basis – a row of that carries that basis alone. For
each basis take the largest share it has in any row,
, with the
basis rescaled so that its columns sum to one (so the value does not depend
on X.restriction); these are anchor.purity, attained at rows
anchor.row. The separability is their minimum: 1 when every basis
has a row of its own, when every row mixes all bases equally, and
NA for rank = 1. The same three values are stored on the fit,
in fit$criterion.
Anchors for every basis, together with , make the
factorization unique up to scale, so a value near 1 says the solution is
nearly unique. A fit with X.anchor scores exactly 1. The measure
looks at only, however: zeros in or in can pin
the solution as well, so a low value does not show that it is not unique.
An object of class summary.nmfkc, containing summary statistics.
Among them, separability, anchor.purity and
anchor.row; see Details.
nmfkc, nmfkc.inference, plot.nmfkc
Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) summary(result)Y <- matrix(cars$dist, nrow = 1) A <- rbind(1, cars$speed) result <- nmfkc(Y, A, rank = 1) summary(result)
Produces a summary of a fitted NMF model with inference results,
including the coefficients table for ().
## S3 method for class 'nmfkc.inference' summary(object, ...)## S3 method for class 'nmfkc.inference' summary(object, ...)
object |
An object of class |
... |
Additional arguments (currently unused). |
An object of class "summary.nmfkc.inference".
nmfkc.inference, summary.nmfkc
Summary method for nmfkc.net objects
## S3 method for class 'nmfkc.net' summary(object, ...)## S3 method for class 'nmfkc.net' summary(object, ...)
object |
An |
... |
Unused. |
An object of class "summary.nmfkc.net".
Produces a summary of a symmetric NMF model with inference results,
including the coefficients table for .
## S3 method for class 'nmfkc.net.inference' summary(object, ...)## S3 method for class 'nmfkc.net.inference' summary(object, ...)
object |
An object of class |
... |
Additional arguments passed to |
An object of class "summary.nmfkc.net.inference".
Summary method for nmfkc.net.signed objects
## S3 method for class 'nmfkc.net.signed' summary(object, ...)## S3 method for class 'nmfkc.net.signed' summary(object, ...)
object |
An |
... |
Unused. |
An object of class "summary.nmfkc.net.signed".
Summary method for nmfkc.signed
## S3 method for class 'nmfkc.signed' summary(object, ...)## S3 method for class 'nmfkc.signed' summary(object, ...)
object |
An |
... |
Unused. |
An object of class "summary.nmfkc.signed".
This function is experimental. The interface may change in future versions.
Ding, C. H. Q., Li, T., & Jordan, M. I. (2010). Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1), 45–55.
nmfkc.signed, nmfkc, nmfkc.signed.cv
nmfre
Displays a concise summary of an NMF-RE model fit, including dimensions, convergence, variance components, and a coefficient table following standard R regression output conventions.
## S3 method for class 'nmfre' summary( object, ci.show = FALSE, max.coef = 20, by = c("covariate", "basis"), ... ) ## S3 method for class 'summary.nmfre' print(x, ...)## S3 method for class 'nmfre' summary( object, ci.show = FALSE, max.coef = 20, by = c("covariate", "basis"), ... ) ## S3 method for class 'summary.nmfre' print(x, ...)
object |
An object of class |
ci.show |
Logical. If |
max.coef |
Largest number of coefficient rows to print (default 20).
The table has one row per (basis, covariate) pair, so it grows as
|
by |
Grouping order of the coefficient rows: |
... |
Additional arguments (currently unused). |
x |
An object of class |
The input object, invisibly.
Y <- matrix(cars$dist, nrow = 1) A <- rbind(intercept = 1, speed = cars$speed) res <- nmfre(Y, A, rank = 1, maxit = 5000) summary(res)Y <- matrix(cars$dist, nrow = 1) A <- rbind(intercept = 1, speed = cars$speed) res <- nmfre(Y, A, rank = 1, maxit = 5000) summary(res)