Fix multiclass logit

NEWS.md

@@ -1,5 +1,20 @@
 # randomPlantedForest 0.3.0.9000
 
+* Fixed multiclass classification with `loss = "logit"`, which effectively never
+  split and predicted near-uniform class probabilities (#40). The C++ logit loss
+  is a reference-class multinomial formulation expecting `K-1` indicator columns,
+  but the R wrapper passed a full `K`-column one-hot matrix, pinning the implicit
+  reference-class probability to zero. Multiclass logit outcomes are now encoded
+  with the first factor level as reference class:
+  * `predict(type = "prob")` reconstructs all `K` class probabilities from the
+    `K-1` logits; rows sum to 1 exactly.
+  * `predict(type = "numeric"/"link")` now returns `K-1` columns named after the
+    non-reference levels (previously `K` columns).
+  * `predict_components()` on multiclass logit fits returns per-class components
+    for the `K-1` non-reference levels; `target_levels` reflects this.
+* The default for `delta` changed from 0 to 0.001. With `delta = 0`, splits
+  producing single-class nodes have infinite logit loss and are always rejected,
+  which also degraded binary logit fits.
 * New `rpf_marshal()` / `rpf_unmarshal()` serialize a fitted forest to a plain
   R list and back, making `saveRDS()`-based storage of rpf models possible (#52).
   Purified forests restore their purified state directly; training data is only

R/predict_components.R

@@ -130,9 +130,11 @@ predict_components <- function(object, new_data, max_interaction = NULL, predict
   )
 
   # Get outcome levels for multiclass handling
+  # (model levels: multiclass logit has no column for the reference level)
   outcome_levels <- levels(object$blueprint$ptypes$outcomes[[1]])
+  model_levels <- model_outcome_levels(object)
   if (length(outcome_levels) > 2) {
-    ret$target_levels <- outcome_levels
+    ret$target_levels <- model_levels
   }
 
   # If max_interaction here is smaller than that of model fit, we calculate a remainder term
@@ -143,7 +145,7 @@ predict_components <- function(object, new_data, max_interaction = NULL, predict
 
     # handling differs if multiclass
     if (length(outcome_levels) > 2) {
-      ret$remainder <- calc_remainders_multiclass(all_components, outcome_levels, pred, intercept)
+      ret$remainder <- calc_remainders_multiclass(all_components, model_levels, pred, intercept)
     } else {
       # regression and binary classif: straight forward.
       ret$remainder <- pred[[1]] - ret$intercept - rowSums(all_components)
@@ -207,7 +209,8 @@ predict_components <- function(object, new_data, max_interaction = NULL, predict
 
   if (length(outcome_levels) > 2) {
     # Multiclass needs disambiguation with one column for each predicted class
-    colnames(ret) <- paste0(out_names, "__class:", outcome_levels)
+    # (model levels: multiclass logit has no column for the reference level)
+    colnames(ret) <- paste0(out_names, "__class:", model_outcome_levels(object))
   } else {
     # Regression and binary classif: Single-column component
     colnames(ret) <- out_names

R/predict_rpf.R

@@ -86,8 +86,8 @@ predict_rpf_numeric <- function(object, new_data, ...) {
   } else {
     # multiclass case for 'link' type: get prediction matrix, clean up
     # via spruce_prob, but no transformation
-    outcome_levels <- levels(object$blueprint$ptypes$outcomes[[1]])
-    out <- hardhat::spruce_prob(outcome_levels, pred)
+    # Multiclass logit: one logit per non-reference level
+    out <- hardhat::spruce_prob(model_outcome_levels(object), pred)
   }
 
   out
@@ -120,7 +120,14 @@ predict_rpf_prob <- function(object, new_data, ...) {
   } else {
     # Multiclass
 
-    if (object$params$loss %in% c("logit", "exponential")) {
+    if (object$params$loss == "logit") {
+      # Raw predictions are K-1 logits relative to the reference class
+      # (first outcome level), so p_ref = 1 / (1 + sum(exp(W))) and
+      # p_k = exp(W_k) * p_ref. Rows sum to 1 by construction.
+      exp_raw <- exp(pred_raw)
+      denom <- 1 + rowSums(exp_raw)
+      pred_prob <- cbind(1 / denom, exp_raw / denom)
+    } else if (object$params$loss == "exponential") {
       # softmax for multi-class
       pred_prob <- softmax(pred_raw)
     } else if (object$params$loss %in% c("L1", "L2")) {

R/rpf.R

@@ -27,12 +27,16 @@
 #' @param loss `["L2"]`: For regression, only `"L2"` is supported. For
 #'   classification, `"L1"`, `"logit"` and `"exponential"` are also available.
 #'   `"exponential"` yields similar results as `"logit"` while being significantly faster.
-#' @param delta `[0]`: Only used if `loss` is `"logit"` or `"exponential"`.
-#'   Proportion of class membership is truncated to be smaller 1-delta when calculating
-#'   the loss to determine the optimal split.
+#' @param delta `[0.001]`: Only used if `loss` is `"logit"` or `"exponential"`.
+#'   Proportion of class membership is truncated to be within `[delta, 1-delta]`
+#'   when calculating the loss to determine the optimal split. Should be positive
+#'   for `"logit"`: with `delta = 0`, nodes containing only a single class produce
+#'   an infinite loss and the corresponding splits are always rejected.
 #' @param epsilon `[0.1]`: Only used if loss = `"logit"` or `"exponential"`.
-#'   Proportion of class membership is truncated to be smaller 1-epsilon when calculating
-#'   the fit in a leaf.
+#'   Proportion of class membership is truncated to be within `[epsilon, 1-epsilon]`
+#'   when calculating the fit in a leaf. Unlike `delta` (a numerical guard for the
+#'   split criterion), this caps the magnitude of individual leaf updates and acts
+#'   as regularization: smaller values permit larger per-leaf jumps on the link scale.
 #' @param export_forest `[FALSE]`: Whether to store the flattened forest in the
 #'   returned object as `$forest`. If `FALSE`, `$forest` is `NULL`, reducing
 #'   memory use of the returned object.
@@ -118,7 +122,7 @@ rpf.data.frame <- function(
   purify = FALSE,
   cv = FALSE,
   loss = "L2",
-  delta = 0,
+  delta = 0.001,
   epsilon = 0.1,
   split_structure = "leaves",
   export_forest = FALSE,
@@ -168,7 +172,7 @@ rpf.matrix <- function(
   purify = FALSE,
   cv = FALSE,
   loss = "L2",
-  delta = 0,
+  delta = 0.001,
   epsilon = 0.1,
   split_structure = "leaves",
   export_forest = FALSE,
@@ -218,7 +222,7 @@ rpf.formula <- function(
   purify = FALSE,
   cv = FALSE,
   loss = "L2",
-  delta = 0,
+  delta = 0.001,
   epsilon = 0.1,
   split_structure = "leaves",
   export_forest = FALSE,
@@ -268,7 +272,7 @@ rpf.recipe <- function(
   purify = FALSE,
   cv = FALSE,
   loss = "L2",
-  delta = 0,
+  delta = 0.001,
   epsilon = 0.1,
   split_structure = "leaves",
   export_forest = FALSE,
@@ -318,7 +322,7 @@ rpf_bridge <- function(
   purify = FALSE,
   cv = FALSE,
   loss = "L2",
-  delta = 0,
+  delta = 0.001,
   epsilon = 0.1,
   split_structure = "leaves",
   export_forest = FALSE

R/utils.R

@@ -210,6 +210,15 @@ preprocess_outcome <- function(processed, loss) {
       outcomes_mat <- stats::model.matrix(~ -1 + outcomes)
       colnames(outcomes_mat) <- levels(outcomes)
 
+      if (loss == "logit") {
+        # The C++ logit loss is a reference-class multinomial formulation:
+        # it models K-1 logits and treats 1 - rowSums(Y) as the implicit
+        # reference class (first factor level), so it needs K-1 indicator
+        # columns. Full one-hot pins the reference-class probability to
+        # zero and no split ever improves the loss (#40).
+        outcomes_mat <- outcomes_mat[, -1, drop = FALSE]
+      }
+
       outcomes <- outcomes_mat
     }
 
@@ -252,6 +261,16 @@ softmax <- function(x) {
 #' @param levels Outcome levels as stored in `rpf$blueprint$ptypes$outcomes`.
 #' @param pred Regular model predictions as returned by `predict.rpf`.
 #' @param intercept Intercept as stored in output of `predict_components`.
+# Outcome levels represented by columns of the C++ prediction matrix.
+# Multiclass logit uses reference-class encoding: the first level has no column.
+model_outcome_levels <- function(object) {
+  outcome_levels <- levels(object$blueprint$ptypes$outcomes[[1]])
+  if (object$params$loss == "logit" && length(outcome_levels) > 2) {
+    outcome_levels <- outcome_levels[-1]
+  }
+  outcome_levels
+}
+
 calc_remainders_multiclass <- function(m, levels, pred, intercept) {
   # data.table NSE warnings
   term <- remainder <- m_sum <- NULL

tests/testthat/test-predict-multiclass.R

@@ -150,7 +150,23 @@ test_that("logit: Numeric/link prediction", {
   classif_pred_lnk <- predict(classif_fit, new_data = xdat, type = "link")
 
   expect_identical(classif_pred, classif_pred_lnk)
-  expect_equal(dim(classif_pred), c(nrow(xdat), nlevels(xdat$yfact)))
+  # Reference-class encoding: raw predictions are K-1 logits relative to the
+  # first factor level, so the reference level has no link column
+  expect_equal(dim(classif_pred), c(nrow(xdat), nlevels(xdat$yfact) - 1))
+  expect_named(classif_pred, paste0(".pred_", levels(xdat$yfact)[-1]))
+})
+
+test_that("logit: prob is consistent with link", {
+  classif_fit <- rpf(yfact ~ ., data = xdat, loss = "logit")
+
+  pred_link <- as.matrix(predict(classif_fit, new_data = xdat, type = "link"))
+  pred_prob <- as.matrix(predict(classif_fit, new_data = xdat, type = "prob"))
+
+  denom <- 1 + rowSums(exp(pred_link))
+  expected <- cbind(1 / denom, exp(pred_link) / denom)
+
+  expect_equal(unname(pred_prob), unname(expected))
+  expect_equal(rowSums(pred_prob), rep(1, nrow(xdat)))
 })
 
 # exponential loss --------------------------------------------------------
@@ -195,6 +211,15 @@ test_that("exponential: Numeric/link prediction", {
   expect_equal(dim(classif_pred), c(nrow(xdat), nlevels(xdat$yfact)))
 })
 
+test_that("logit: predict_components uses non-reference levels", {
+  classif_fit <- rpf(yfact ~ ., data = xdat, loss = "logit", max_interaction = 2)
+
+  components <- predict_components(classif_fit, xdat, max_interaction = 1)
+
+  expect_identical(components$target_levels, levels(xdat$yfact)[-1])
+  expect_named(components$remainder, levels(xdat$yfact)[-1])
+})
+
 # Classif and prob agree with each other ----------------------------------
 
 test_that("prob and classif yield same result", {

tests/testthat/test-rpf-classification-multiclass.R

@@ -31,6 +31,20 @@ test_that("Multiclass: Detection works", {
   expect_error(rpf(ychar ~ x3 + x4, xdat), regexp = "Ordering of factor columns only implemented")
 })
 
+# Multiclass logit learns signal (#40) -------------------------------------
+test_that("Multiclass logit learns signal (#40)", {
+  set.seed(42)
+  n <- 400
+  dat <- data.frame(x1 = runif(n), x2 = runif(n))
+  dat$y <- factor(ifelse(dat$x1 > 0.5, "a", ifelse(dat$x2 > 0.5, "b", "c")))
+  idx <- sample(n, 280)
+
+  fit <- rpf(y ~ ., data = dat[idx, ], loss = "logit", ntrees = 50, max_interaction = 2)
+  pred <- predict(fit, dat[-idx, ], type = "class")
+
+  expect_gt(mean(pred$.pred_class == dat$y[-idx]), 0.85)
+})
+
 test_that("Remainder is calculcated correctly", {
   classif_fit <- rpf(yfact ~ ., data = xdat, max_interaction = 3)