From 757d82df5720ead13ce0071f20296c3c2e96cd1a Mon Sep 17 00:00:00 2001 From: Lukas Burk Date: Wed, 15 Jul 2026 16:15:30 +0200 Subject: [PATCH] fix R interface to multiclass logit loss --- NEWS.md | 15 +++++++++++ R/predict_components.R | 9 ++++--- R/predict_rpf.R | 13 ++++++--- R/rpf.R | 24 ++++++++++------- R/utils.R | 19 +++++++++++++ man/rpf.Rd | 22 ++++++++------- tests/testthat/test-predict-multiclass.R | 27 ++++++++++++++++++- .../test-rpf-classification-multiclass.R | 14 ++++++++++ 8 files changed, 117 insertions(+), 26 deletions(-) diff --git a/NEWS.md b/NEWS.md index 29a971f..e7fdcde 100644 --- a/NEWS.md +++ b/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 diff --git a/R/predict_components.R b/R/predict_components.R index bf6be47..92fe647 100644 --- a/R/predict_components.R +++ b/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 diff --git a/R/predict_rpf.R b/R/predict_rpf.R index 54f2988..85ff747 100644 --- a/R/predict_rpf.R +++ b/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")) { diff --git a/R/rpf.R b/R/rpf.R index d02d13a..0703223 100644 --- a/R/rpf.R +++ b/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 diff --git a/R/utils.R b/R/utils.R index 259edd7..82ff269 100644 --- a/R/utils.R +++ b/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 diff --git a/man/rpf.Rd b/man/rpf.Rd index 93373f7..aec1a0c 100644 --- a/man/rpf.Rd +++ b/man/rpf.Rd @@ -26,7 +26,7 @@ rpf(x, ...) purify = FALSE, cv = FALSE, loss = "L2", - delta = 0, + delta = 0.001, epsilon = 0.1, split_structure = "leaves", export_forest = FALSE, @@ -49,7 +49,7 @@ rpf(x, ...) purify = FALSE, cv = FALSE, loss = "L2", - delta = 0, + delta = 0.001, epsilon = 0.1, split_structure = "leaves", export_forest = FALSE, @@ -72,7 +72,7 @@ rpf(x, ...) purify = FALSE, cv = FALSE, loss = "L2", - delta = 0, + delta = 0.001, epsilon = 0.1, split_structure = "leaves", export_forest = FALSE, @@ -95,7 +95,7 @@ rpf(x, ...) purify = FALSE, cv = FALSE, loss = "L2", - delta = 0, + delta = 0.001, epsilon = 0.1, split_structure = "leaves", export_forest = FALSE, @@ -145,13 +145,17 @@ Can be achieved after fitting with \code{\link[=purify]{purify()}}.} classification, \code{"L1"}, \code{"logit"} and \code{"exponential"} are also available. \code{"exponential"} yields similar results as \code{"logit"} while being significantly faster.} -\item{delta}{\verb{[0]}: Only used if \code{loss} is \code{"logit"} or \code{"exponential"}. -Proportion of class membership is truncated to be smaller 1-delta when calculating -the loss to determine the optimal split.} +\item{delta}{\verb{[0.001]}: Only used if \code{loss} is \code{"logit"} or \code{"exponential"}. +Proportion of class membership is truncated to be within \verb{[delta, 1-delta]} +when calculating the loss to determine the optimal split. Should be positive +for \code{"logit"}: with \code{delta = 0}, nodes containing only a single class produce +an infinite loss and the corresponding splits are always rejected.} \item{epsilon}{\verb{[0.1]}: Only used if loss = \code{"logit"} or \code{"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 \verb{[epsilon, 1-epsilon]} +when calculating the fit in a leaf. Unlike \code{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.} \item{split_structure}{\verb{["leaves"]}: Defines the structure of a possible split and how to choose split_candidates. Can be one of "leaves", "hist", "cur_trees_1", "cur_trees_2", or "res_trees". Further details are given below.} diff --git a/tests/testthat/test-predict-multiclass.R b/tests/testthat/test-predict-multiclass.R index 2fc5c94..54955a6 100644 --- a/tests/testthat/test-predict-multiclass.R +++ b/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", { diff --git a/tests/testthat/test-rpf-classification-multiclass.R b/tests/testthat/test-rpf-classification-multiclass.R index 3d0b357..94accea 100644 --- a/tests/testthat/test-rpf-classification-multiclass.R +++ b/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)