NMF-KC with signed covariate matrix

Solves $$Y \approx X\,\Theta\,A,\qquad X \ge 0,\;\Theta\in\R^{Q\times D}, \;A\in\R^{D\times N},$$ where the covariate matrix $A$ and the coefficient matrix $\Theta$ may be signed. Internally $A = A_{+} - A_{-}$ and $\Theta = C_{+} - C_{-}$ with $A_{\pm}, C_{\pm} \ge 0$ (sign-splitting trick, Ding et al. 2010), and the problem is solved by a Direct Multiplicative Update algorithm whose iteration cost is $O(Q D^2)$, independent of $N$.

Only $X$ is structurally constrained to be non-negative (Semi-NMF sense of Ding, Li, & Jordan 2010). In particular, $Y$ may contain negative entries, in which case the response is fit in the least-squares sense without any non-negativity requirement on $Y$.

When $A \ge 0$ (so $A_{-} = 0$), the result reduces to nmfkc(Y, A, rank) with Euclidean loss, up to reordering.

Usage

nmfkc.signed(
  Y,
  A,
  rank = NULL,
  epsilon = 1e-04,
  maxit = 5000,
  verbose = TRUE,
  ...
)

Arguments

Y

Real-valued $Q_{\mathrm{obs}} \times N$ response matrix. Unlike nmfkc, negative entries are allowed.

A

Real-valued $D \times N$ covariate matrix (signed). A single matrix is passed; its positive and negative parts $A_{+} = \max(A, 0)$ and $A_{-} = \max(-A, 0)$ are computed internally. When using Random Fourier Features (Rahimi & Recht 2007) as $A$, supply the RFF parameters via the hidden pars argument so that predict() can regenerate features for new data (see pars entry in ... below).

Large $N$: A may instead be a Gram object of class "nmfkc.gram" returned by nmfkc.signed.rff.gram, which holds only $S = AA^\top$ ($D \times D$) and $G_0 = YA^\top$ ($Q_{\mathrm{obs}} \times D$) accumulated over column blocks, so the $D \times N$ feature matrix never has to exist in memory. The multiplicative updates are unchanged (they only ever use $S$ and $G_0$), and the fit is the same as with the matrix up to floating-point summation order. Three things differ with Gram input: warm.start is unavailable (the posneg warm-start needs the $2D \times N$ split matrix) and the direct initialization is used with a message; Y.weights and NA in Y are errors (the weighted loss needs $A$ on every iteration); and the fold-based helpers nmfkc.signed.cv / nmfkc.signed.ecv / nmfkc.signed.rank do not accept it – select tuning parameters on a held-out validation set, scoring with predict.nmfkc.signed on features regenerated by nmfkc.signed.rff with the stored pars. B, XB and the residual statistics are rebuilt block by block from the object's feature generator.

rank

Integer. Number of latent components $Q$ in $X$.

epsilon

Relative convergence tolerance on the objective (default 1e-4).

maxit

Maximum number of iterations (default 5000).

verbose

Logical. Print dimensions at start (default TRUE).

...

Additional arguments:

Q: alias for rank.
X.restriction: normalization applied to $X$ after every update. One of "colSums" (default, $\mathrm{colSums}(X) = 1$), "colSqSums", "totalSum", "none", "fixed". All of them are gauge fixes: the scale is divided out of $X$ and multiplied into $C_{+}, C_{-}$, so $XC$ and hence the fit is unchanged and only the parametrization is pinned down. "rowSums" was removed in 0.9.8 and is now refused: scaling the rows of $X$ changes $XCA$, so it was a restriction acting rather than a reparametrization, and the objective oscillated instead of descending.
X.L2.ortho: non-negative L2 orthogonality penalty on the columns of $X$ (default 0), penalizing $(\lambda/2)\lVert \mathrm{offdiag}(X^\top X)\rVert^2$. Same convention as nmfkc; skipped when X.restriction = "fixed".
X.L2.smooth: non-negative L2 row-smoothness penalty on $X$ (default 0), penalizing $(\lambda/2)\,\mathrm{tr}(X^\top L X)$ with $L$ the path-graph Laplacian over rows (adjacent-row differences). Useful for ordered rows (e.g. time / space); skipped when X.restriction = "fixed".
C.L2: non-negative ridge penalty on the signed coefficient matrix $C = C_{+} - C_{-}$ (default 0), adding $\lambda\lVert C_{+} - C_{-}\rVert^2$. Shrinks $\Theta$ toward zero (with zero gradient on the unidentified common mode $C_{+} + C_{-}$), injected into both the fast unweighted and the weighted $C_{+}$/$C_{-}$ updates.
C.L1: non-negative lasso penalty on the signed coefficient matrix (default 0), adding $\lambda\lVert C\rVert_1$. Under the split $C = C_{+} - C_{-}$ this is $\lambda\sum(C_{+}+C_{-})$: at a minimizer $C_{+}$ and $C_{-}$ have disjoint supports, so the sum equals the absolute value, and the constant gradient enters only the denominators of both updates (same $\lambda/2$ convention as nmfkc). Multiplicative updates approach zero geometrically rather than reaching it, so the effect is shrinkage; threshold the result if a sparse support is wanted.
update.power: exponent applied to the multiplicative ratio in the unweighted sweep (default 1). With 1 the updates are of the Lee–Seung form and coincide with those of nmfkc when $A \ge 0$; with 0.5 they are the square-root form of Ding, Li and Jordan (2010), which is the form their monotonicity proof for the nonnegative factor under a signed Gram matrix covers. Both have the same fixed points; the root takes smaller steps and needs more iterations. Ignored when Y.weights is supplied (the weighted sweep has its own damping with backtracking).
X.init: initialization strategy for the basis matrix $X$ ($Q_{\mathrm{obs}} \times Q$). Accepts the same menu as nmfkc: "kmeans" (default), "kmeansar", "nndsvd", "runif", or a user-supplied $Q_{\mathrm{obs}} \times Q$ non-negative numeric matrix. String methods delegate to the shared internal helper .init_X_method() (see nmfkc for the definitions of each method). For signed $Y$, "kmeans" cluster centers may contain negative entries; they are clipped to zero to satisfy $X \ge 0$, and any column that collapses to all-zeros is re-filled with small $\mathrm{Uniform}(0, 0.1)$ noise.
C.init: explicit initial $Q \times D$ coefficient matrix $\Theta$ (signed). Split internally.
warm.start: logical (default TRUE). If TRUE and $Y \ge 0$, runs nmfkc(Y, A = rbind(A_+, A_-), rank = Q) internally to seed $X, C_{+}, C_{-}$. The user's X.init, seed, nstart, and X.restriction are forwarded to the internal nmfkc call so that initialization choices propagate consistently between the warm-start and the signed MU loop. Ignored when $Y$ has negative entries (warm-start is disabled; X.init is used directly by the signed branch instead).
seed: RNG seed for random initialization (default 123).
prefix: name prefix for rows of $C$ and columns of $X$ (default "Basis").
pars: optional list list(omega, b, D, beta) of Random Fourier Feature parameters (Rahimi & Recht 2007; omega: frequency matrix, b: phase offset, D: feature dimension, beta: bandwidth). When supplied, it is stored in the returned object so that summary() can report $\beta$ and downstream predict() calls can regenerate RFF features for new data. If A is not RFF features, leave this NULL.
Y.weights: Optional non-negative weight matrix ($Q_{\mathrm{obs}} \times N$) or vector (length $N$), analogous to the weights argument of lm. Loss becomes $\sum W_{ij} \, (Y_{ij} - (XCA)_{ij})^2$ (lm()-style, linear in $W$). Logical matrices (TRUE / FALSE) are also accepted. Typical usage by nmfkc.signed.cv / nmfkc.signed.ecv passes a binary mask $W \in \{0,1\}$ to hold out test elements; real-valued weights for observation-level importance weighting are also supported. Default NULL: if Y has NA, a binary mask is auto-constructed (0 for NA, 1 elsewhere); otherwise no weighting.
nstart.signed: number of restarts of the signed fit (default 1). Signed models have more local minima than non-negative ones because $\Theta = C_{+} - C_{-}$ can take both positive and negative values. With nstart.signed > 1 the whole fit is repeated from that many consecutive seeds (seed, seed + 1, ...) and the fit with the smallest $objfunc is returned; the returned object then also carries $restarts, a data frame of the seed, objective, iteration count and convergence flag of every start. Note that this is a different quantity from nstart, which – here as in nmfkc – only governs the initialization of $X$ and is forwarded to the non-negative warm-start fit. A budget of 10-50 restarts is recommended for publication-grade runs, especially when the number of classes is large. Restarts are cheap for Gram input, where $S$ and $G_0$ are already accumulated.
cores: number of workers used for the restarts when nstart.signed > 1 (default 1, i.e. sequential).

Value

An object of class c("nmfkc.signed", "nmfkc") with

X: $Q_{\mathrm{obs}} \times Q$ basis matrix (non-negative, column-normalized according to X.restriction).
Cp, Cn: $Q \times D$ non-negative parts of $\Theta$, so that $\Theta = C_{+} - C_{-}$.
C: $C_{+} - C_{-}$ (= $\Theta$), signed.
B: $C \, A$, $Q \times N$ (signed).
objfunc.iter: objective values per iteration.
objfunc: final objective.
r.squared: $\mathrm{cor}(Y, \widehat Y)^2$ (Pearson; in $[0,1]$).
r.squared.uncentered: uncentered $R^2 = 1 - \|Y - \widehat Y\|_F^2 / \|Y\|_F^2$ (baseline = zero matrix).
r.squared.centered: row-mean centered $1 - \|Y - \widehat Y\|_F^2 / \|Y - \bar Y_{p\cdot}\|_F^2$.
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 $Y$ 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.

Lifecycle

This function is experimental. The interface may change in future versions.

References

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.

Examples

# \donttest{
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)
#> Y(8,40) ~ X(8,3) %*% C(3,20) %*% A(20,40)  [signed covariate]

## 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)
#> Y(8,40) ~ X(8,3) %*% C(3,20) %*% A(20,40)  [signed covariate], Y signed
# }