bug: graph penalty does not affect gp loss

[travis-leith]

min rep (generated by AI)

# Minimal reproducible example: graph-Laplacian penalty (graph_W / graph_lambda) is a no-op
# for IRLS losses (loss = "gp"), but works for the Gaussian solver (loss = "mse").
#
# RcppML build: GitHub zdebruine/RcppML, version 1.0.0, RemoteSha df69dddbe37962cd0179b6b6c6ace94e400b4f60
#
# Expectation: a nonzero graph_lambda over a nonempty feature graph should change the fitted
# factor W. It does under loss = "mse", but leaves W BIT-FOR-BIT IDENTICAL under loss = "gp",
# even at an extreme graph_lambda. The `angular` penalty, by contrast, DOES take effect under
# loss = "gp", so it is the graph penalty specifically that is not wired into the IRLS path.

library(RcppML)

cat("RcppML version:", as.character(packageVersion("RcppML")), "\n\n")

# Built-in dense count matrix (824 features x 135 samples); no external data needed.
data(aml)
A <- as.matrix(aml)

K <- 4L
SEED <- 1L
LAMBDA <- 50 # deliberately large so any effect is unmissable

# A chain graph over ADJACENT FEATURES (rows of A) -- a valid nonzero feature adjacency.
n_feat <- nrow(A)
feat_chain <- (abs(row(diag(n_feat)) - col(diag(n_feat))) == 1) * 1
cat("feature graph:", nrow(feat_chain), "x", ncol(feat_chain),
    "with", sum(feat_chain), "nonzeros\n\n")

fit_nmf <- function(loss, ...) {
    RcppML::nmf(A, k = K, loss = loss, seed = SEED, verbose = FALSE, ...)
}

# --- Graph penalty under the Gaussian (MSE) solver: EXPECTED to change W --------------------
mse_plain <- fit_nmf("mse")
mse_graph <- fit_nmf("mse", graph_W = feat_chain, graph_lambda = LAMBDA)
cat(sprintf("mse loss: max|W_plain - W_graph(%.0f)| = %.6g   <- graph penalty active\n",
            LAMBDA, max(abs(mse_plain@w - mse_graph@w))))

# --- Graph penalty under the IRLS (gp) solver: NO EFFECT (the bug) --------------------------
gp_plain <- fit_nmf("gp")
gp_graph <- fit_nmf("gp", graph_W = feat_chain, graph_lambda = LAMBDA)
cat(sprintf("gp  loss: max|W_plain - W_graph(%.0f)| = %.6g   <- graph penalty IGNORED\n",
            LAMBDA, max(abs(gp_plain@w - gp_graph@w))))

# --- Control: `angular` DOES take effect under gp, so the IRLS path is reached ---------------
gp_angular <- fit_nmf("gp", angular = c(0.5, 0))
cat(sprintf("gp  loss: max|W_plain - W_angular(0.5)| = %.6g   <- angular penalty active\n",
            max(abs(gp_plain@w - gp_angular@w))))

Either this should be supported, or an explicit error should be raised indicating that it is not supported.