Classification: a level for the logit and exponential losses, and two fixes for tree families without a split

NEWS.md

@@ -69,6 +69,13 @@
 
 ## New features
 
+* Classification forests with `loss = "logit"` or `"exponential"` get a level: after the fit, one constant per score column, the maximum-likelihood intercept of the log-loss given the forest's out-of-bag scores, so that every class's mean out-of-bag probability is its frequency in the training rows.
+  It is stored in every tree family's null leaf, so it is part of the intercept of the forest, of its purification and of its pruned model; `$level` holds it and `print()` shows it (`?rpf`, "The level of the classification losses").
+  The families are averaged on the score scale, and where they disagree the probability of their mean score is below the mean of their probabilities for a class rarer than 1/2, so the forest predicted a rare class too rarely, the more so the more splits it made.
+  With the level the forest predicts each class about as often as it occurs, at any number of splits.
+  The out-of-bag choice of `n_splits` scores every number of splits at its own level, and the pruning starts from the forest with its level.
+  This changes every fit with these losses, binary and multiclass; `options(randomPlantedForest.level = FALSE)` fits without the level.
+
 * `rpf()` gains `keep_level_order`, to keep the level order of factor predictors instead of ordering their levels by the outcome.
   `TRUE` keeps it for every factor, ordered or not, and a character vector for the named factors only; splits then separate neighbouring levels in the order of `levels()`.
   The default, `FALSE`, orders every factor by the outcome as before, and fits are unchanged.
@@ -153,6 +160,12 @@
 
 ## Bug fixes
 
+* `prune_components()` (and so `rpf()` with pruning) left the value of a tree family that made no split out of the pruned model's intercept, so that rule `"D"` at `tau = 0` was not the forest where a family had found no split. The intercept now includes it, as `purify()`'s does.
+
+* The out-of-bag path (`oob_path`, and the out-of-bag choice of `n_splits`) counted a tree family that had made no split at the score 0 its rows start from, instead of at the fit of its rows that the family predicts (the mean for regression and the L1 / L2 losses, the log-odds of its sample for `"logit"` and `"exponential"`, so a probability of 1/2 instead).
+  It now counts such a family at its fit.
+  Families without a split are common only with a rare class and few rows, where they made the out-of-bag log-loss of few splits too high, so that the out-of-bag choice of `n_splits` could take far more splits than it should.
+
 * `split_structure = "tree_then_leaf"` samples leaves in proportion to their available thresholds, as documented.
   After a split, the sampling weights of new leaves counted rows instead of distinct values, and the cached prefix sums were built incorrectly.
   Leaves with tied values were then drawn too often, and draws were wasted on leaves that cannot be split.

R/criterion.R

@@ -118,6 +118,247 @@ crit_subset <- function(crit, idx) {
   crit
 }
 
+# The classification level. The tree families are averaged on the score (logit) scale, and the mean of their
+# probabilities lies above the probability of their mean score where the probability is convex in the score (a rare
+# class: Jensen), more so the more the families disagree, i.e. with more splits. So a forest predicts a rare class too
+# rarely. Its level is the constant c (one value per score column) that the forest's score gets after the fit (the C++
+# shift_level(), in every family's null leaf, so part of the intercept): the maximum-likelihood intercept of the
+# log-loss given the forest's out-of-bag scores, whose equations are the balance property
+#   sum_i p_k(c + g_i) - Delta_ik = n_k   (every free class k: the mean probability is the class's frequency)
+# over the rows `idx` with an out-of-bag family, g the OOB mean scores (G; the half means G0, G1, half counts m0, m1),
+# B the families. Delta is the second-order term of the OOB correction (criterion.R, above): an OOB mean of m families
+# scatters more than the forest's mean of B, which moves its expected probability by
+#   Delta_k = (1/m - 1/B) / 2 * sum_uv d2 p_k / dg_u dg_v V_uv = (1/m - 1/B) / 2 * w p_k ((d_k - p'd)^2 - (p'd^2 - (p'd)^2))
+# with V = w d d' the across-family covariance from the halves (d over all K classes, the reference class's 0); rows
+# without both halves get none. Solved by Newton's method from c = 0 with the Jacobian of the first term,
+# sum_i diag(p_i) - p_i p_i' (positive definite), halving the step while the residual's norm does not fall; it stops at
+# a largest residual below 1e-10 of a row (usually after 3 to 5 steps). Free classes: binary the second class;
+# multiclass logit the K - 1 classes after the first (the reference, score 0); multiclass exponential (K scores, the
+# softmax) the classes after the first, the first's score then fixed at 0, and c is centred at the end (the softmax
+# ignores what all classes share; SAMME's scores are centred). 0 where a free class or the reference class has no
+# out-of-bag row (no finite level). Returns V values, with the attribute `iterations`.
+#
+# Its arithmetic is plain double arithmetic in a fixed order (level_colsums(), level_rowsums(), level_solve()), so the
+# level's bits depend on its inputs and on exp() and log() only: no BLAS or LAPACK, whose crossprod() and solve() round
+# differently by library (OpenBLAS's triangular solve multiplies by the pivot's reciprocal where the reference BLAS
+# divides by it), and no sum(), colSums() or rowSums(), which add in long double (80 bits on x86_64, 64 on arm64 macOS).
+# Then the stopping test and the step halving take the same branches on every machine too.
+classif_level <- function(crit, idx, G, G0, G1, m0, m1, B) {
+  G <- as.matrix(G)
+  V <- crit$V
+  n <- nrow(G)
+  cls <- crit$cls[idx]
+  K <- crit$K
+  zero <- structure(numeric(V), iterations = 0L)
+  if (crit$type != "logloss" || n == 0) {
+    return(zero)
+  }
+  n_k <- tabulate(cls, K)
+  if (any(n_k == 0)) {
+    return(zero)
+  }
+  # the free classes and their score columns (exponential's softmax: the first score fixed)
+  free <- seq_len(K)[-1]
+  free_col <- if (crit$link == "softmax") free else free - 1L
+  D <- as.matrix(G0) - as.matrix(G1)
+  m <- m0 + m1
+  both <- m0 > 0 & m1 > 0
+  coef <- (1 / m - 1 / B) * 0.5 / (1 / pmax(m0, 1) + 1 / pmax(m1, 1))
+  coef[!both] <- 0
+  # the full probabilities (n x K) and the score differences of the halves over all K classes (reference 0)
+  full <- function(M) {
+    if (crit$link == "softmax") M else cbind(0, M)
+  }
+  Df <- full(D)
+  probs <- function(cv) {
+    Gc <- G + rep(cv, each = n)
+    if (crit$link == "sigmoid") {
+      p <- stats::plogis(Gc[, 1])
+      return(cbind(1 - p, p))
+    }
+    Gf <- full(Gc)
+    # log-sum-exp of the rows (crit_lse()), the sum in a fixed order
+    mx <- apply(Gf, 1, max)
+    exp(Gf - (mx + log(level_rowsums(exp(Gf - mx)))))
+  }
+  resid <- function(Pf) {
+    pd <- level_rowsums(Pf * Df)
+    vp <- level_rowsums(Pf * Df^2) - pd^2
+    Delta <- coef * Pf * ((Df - pd)^2 - vp)
+    n_k[free] - level_colsums(Pf[, free, drop = FALSE] - Delta[, free, drop = FALSE])
+  }
+  # the Jacobian sum_i diag(p_i) - p_i p_i' of the free classes' probabilities P: entry (u, v) sums p_u (1[u = v] - p_v)
+  jacobian <- function(P) {
+    nf <- ncol(P)
+    uv <- which(upper.tri(diag(nf), diag = TRUE), arr.ind = TRUE)
+    terms <- vapply(
+      seq_len(nrow(uv)),
+      function(j) {
+        u <- uv[j, 1]
+        v <- uv[j, 2]
+        P[, u] * (as.numeric(u == v) - P[, v])
+      },
+      numeric(nrow(P))
+    )
+    s <- level_colsums(matrix(terms, nrow(P)))
+    J <- matrix(0, nf, nf)
+    J[uv] <- s
+    J[uv[, 2:1, drop = FALSE]] <- s
+    J
+  }
+  cv <- numeric(V)
+  Pf <- probs(cv)
+  r <- resid(Pf)
+  it <- 0L
+  while (max(abs(r)) > 1e-10 * n && it < 50L) {
+    it <- it + 1L
+    step <- level_solve(jacobian(Pf[, free, drop = FALSE]), r)
+    if (is.null(step) || !all(is.finite(step))) {
+      break
+    }
+    t <- 1
+    improved <- FALSE
+    for (h in seq_len(30)) {
+      cv_new <- cv
+      cv_new[free_col] <- cv[free_col] + t * step
+      Pf_new <- probs(cv_new)
+      r_new <- resid(Pf_new)
+      if (all(is.finite(r_new)) && level_sum(r_new^2) < level_sum(r^2)) {
+        improved <- TRUE
+        break
+      }
+      t <- t / 2
+    }
+    if (!improved) {
+      break
+    }
+    cv <- cv_new
+    Pf <- Pf_new
+    r <- r_new
+  }
+  if (!all(is.finite(cv))) {
+    return(zero)
+  }
+  if (crit$link == "softmax") {
+    cv <- cv - level_sum(cv) / length(cv)
+  }
+  structure(cv, iterations = it)
+}
+
+# The level's sums (classif_level()), in double arithmetic in a fixed order, the same on every machine. Column sums
+# of an n x k matrix: pairwise, the first half of the rows plus the second half (the odd row out kept last) until one
+# row is left.
+level_colsums <- function(M) {
+  M <- as.matrix(M)
+  if (nrow(M) == 0L) {
+    return(numeric(ncol(M)))
+  }
+  while (nrow(M) > 1L) {
+    h <- nrow(M) %/% 2L
+    top <- M[seq_len(h), , drop = FALSE] + M[h + seq_len(h), , drop = FALSE]
+    M <- if (nrow(M) %% 2L == 1L) rbind(top, M[nrow(M), , drop = FALSE]) else top
+  }
+  M[1, ]
+}
+
+# The sum of a vector (level_colsums())
+level_sum <- function(x) {
+  level_colsums(matrix(as.numeric(x), ncol = 1L))
+}
+
+# Row sums of an n x k matrix: the columns added from the first
+level_rowsums <- function(M) {
+  s <- M[, 1]
+  for (j in seq_len(ncol(M))[-1]) {
+    s <- s + M[, j]
+  }
+  s
+}
+
+# The solution of J x = r for the level's Newton step (J the k x k Jacobian, positive definite): Gaussian elimination
+# with the diagonal as pivots in their order, in double arithmetic. NULL if a pivot is not positive and finite (J
+# numerically singular), as solve() fails.
+level_solve <- function(J, r) {
+  k <- length(r)
+  A <- J
+  b <- as.numeric(r)
+  for (p in seq_len(k)) {
+    if (!is.finite(A[p, p]) || A[p, p] <= 0) {
+      return(NULL)
+    }
+    for (i in seq_len(k - p) + p) {
+      l <- A[i, p] / A[p, p]
+      for (j in seq_len(k - p) + p) {
+        A[i, j] <- A[i, j] - l * A[p, j]
+      }
+      b[i] <- b[i] - l * b[p]
+    }
+  }
+  x <- numeric(k)
+  for (i in rev(seq_len(k))) {
+    s <- b[i]
+    for (j in seq_len(k - i) + i) {
+      s <- s - A[i, j] * x[j]
+    }
+    x[i] <- s / A[i, i]
+  }
+  x
+}
+
+# Whether a fit with these params gets the classification level (classif_level()): rpf()'s option
+# randomPlantedForest.level (TRUE by default) is stored in params$level; a forest fitted before the level existed has
+# no params$level and none. Classification with loss = "logit" or "exponential".
+level_on <- function(params, mode) {
+  mode == "classification" && isTRUE(params$level) && isTRUE(params$loss %in% c("logit", "exponential"))
+}
+
+# The level of the fitted forest `fit` (the C++ object, before shift_level()) with the classification criterion
+# `crit` of its training rows: from its OOB path's last checkpoint, at its last split (the fit tracks it,
+# rpf_impl()); for a forest without out-of-bag rows (deterministic) the balance on its own predictions of the
+# training rows X. B: the families.
+classif_fit_level <- function(fit, crit, X, B, deterministic) {
+  in_sample <- function() {
+    G <- matrix(fit$predict_matrix(X, 0, 0L), ncol = crit$V)
+    n <- nrow(G)
+    classif_level(crit, seq_len(n), G, G, G, rep(1L, n), rep(0L, n), 1)
+  }
+  if (deterministic) {
+    return(in_sample())
+  }
+  path <- fit$oob_path()
+  m0 <- path$count[, 1]
+  m1 <- path$count[, 2]
+  ok <- m0 + m1 > 0
+  if (!any(ok)) {
+    return(in_sample())
+  }
+  V <- crit$V
+  cols <- (length(path$checkpoints) - 1L) * V + seq_len(V)
+  s0 <- path$sum0[ok, cols, drop = FALSE]
+  s1 <- path$sum1[ok, cols, drop = FALSE]
+  m0 <- m0[ok]
+  m1 <- m1[ok]
+  classif_level(crit, which(ok), (s0 + s1) / (m0 + m1), s0 / pmax(m0, 1), s1 / pmax(m1, 1), m0, m1, B)
+}
+
+# The level of a forest with fixed splits, after its fit (rpf_impl() tracked its OOB path at its last split):
+# classif_fit_level(), then classif_level_finish(). Returns the level.
+classif_level_fixed <- function(fit, params, crit, X) {
+  level <- classif_fit_level(fit, crit, X, params$ntrees, isTRUE(params$deterministic))
+  classif_level_finish(fit, params, level)
+  level
+}
+
+# A forest after its fit with the level: shift_level(), then the purification that rpf(purify = TRUE) asked for (the
+# C++ fit does not purify such a forest itself, as the level comes after it; rpf_impl())
+classif_level_finish <- function(fit, params, level) {
+  fit$shift_level(as.numeric(level))
+  if (isTRUE(params$purify)) {
+    fit$purify_threads(0L, 0L, 2L)
+  }
+  invisible(fit)
+}
+
 # Log-sum-exp of the rows of the score columns G (n x V), with the reference score 0 joined if `ref`
 crit_lse <- function(G, ref) {
   mx <- apply(G, 1, max)

R/oob_path.R

@@ -90,7 +90,7 @@ oob_path_point <- function(s0, s1, count, y, ntrees, rows = FALSE) {
 # oob_path_point() for a classification criterion `crit` (oob_score_y()): the mean corrected and plain criterion over
 # the rows with an OOB family and their standard errors, from s0 / s1 (n x V: the sums of the families of even / odd
 # index) and count. Rows whose OOB families are all even or all odd are not corrected (as for regression). rows: as
-# for oob_path_point().
+# for oob_path_point(). With the classification level the point has the attribute `level`, the forest's level there.
 oob_path_point_crit <- function(s0, s1, count, crit, ntrees, rows = FALSE) {
   m0 <- count[, 1]
   m1 <- count[, 2]
@@ -99,18 +99,23 @@ oob_path_point_crit <- function(s0, s1, count, crit, ntrees, rows = FALSE) {
   n_ok <- sum(ok)
   s0 <- s0[ok, , drop = FALSE]
   s1 <- s1[ok, , drop = FALSE]
-  r <- crit_rows(
-    crit,
-    which(ok),
-    (s0 + s1) / m[ok],
-    s0 / pmax(m0[ok], 1),
-    s1 / pmax(m1[ok], 1),
-    m0[ok],
-    m1[ok],
-    ntrees,
-    impute = FALSE
-  )
+  G <- (s0 + s1) / m[ok]
+  G0 <- s0 / pmax(m0[ok], 1)
+  G1 <- s1 / pmax(m1[ok], 1)
+  # with the classification level (crit$level, rpf_fit()) the forest is scored at its level here (classif_level()):
+  # the forest with these splits is returned with it
+  level <- NULL
+  if (isTRUE(crit$level)) {
+    level <- classif_level(crit, which(ok), G, G0, G1, m0[ok], m1[ok], ntrees)
+    shift <- rep(as.numeric(level), each = n_ok)
+    G <- G + shift
+    G0 <- G0 + shift
+    G1 <- G1 + shift
+  }
+  r <- crit_rows(crit, which(ok), G, G0, G1, m0[ok], m1[ok], ntrees, impute = FALSE)
   point <- c(mean(r$corr), mean(r$e), stats::sd(r$corr) / sqrt(n_ok), stats::sd(r$e) / sqrt(n_ok))
+  # (the level, for the adaptive fit's forest)
+  attr(point, "level") <- level
   if (rows) list(point = point, rows = r$corr) else point
 }
 

R/oob_stop.R

@@ -65,6 +65,11 @@ tune_splits_oob <- function(
     )
   }
   res <- st$result()
+  # with the classification level the forest gets its level at the chosen checkpoint, as the fit with that many splits
+  # gets it (rpf_impl())
+  if (!is.null(res$level)) {
+    classif_level_finish(fit, params, res$level)
+  }
   if (res$capped) {
     cli::cli_warn(
       c(
@@ -85,7 +90,8 @@ tune_splits_oob <- function(
 # the lowest error) to stop: once the patience or the plateau and the floor allow it, or at the last checkpoint. The
 # error is the corrected OOB mean squared error, or the mean corrected criterion. result() gives the chosen number of
 # splits, the OOB path's table of the checkpoints grown (as oob_path_table() makes it), where growing stopped, and
-# whether it stopped at the cap with the error still falling (neither the patience nor the plateau reached).
+# whether it stopped at the cap with the error still falling (neither the patience nor the plateau reached), and with
+# the classification level (a criterion with `level`) the level at the chosen checkpoint (otherwise NULL).
 #
 # The plateau (tolerance > 0, from checkpoint patience + 1 on): with b the lowest error up to c and a the lowest up to
 # c - patience (the first of ties, as oob_select()), the gain e[a] - e[b] over the last `patience` checkpoints is at
@@ -107,6 +113,8 @@ oob_stopper <- function(y, ntrees, checkpoints, patience = 4L, tolerance = 0, fl
   st$no_oob <- FALSE
   st$points <- matrix(NA_real_, C, 4)
   st$rows <- vector("list", C)
+  # (the classification level: the forest's level at each checkpoint, classif_level())
+  st$levels <- vector("list", C)
   st$step <- function(c, path) {
     if (c != st$c + 1L) {
       cli::cli_abort("The adaptive fit's checkpoints came out of order.", .internal = TRUE)
@@ -130,6 +138,9 @@ oob_stopper <- function(y, ntrees, checkpoints, patience = 4L, tolerance = 0, fl
       st$rows[[c]] <- pt$rows
       pt <- pt$point
     }
+    if (!is.null(attr(pt, "level"))) {
+      st$levels[[c]] <- attr(pt, "level")
+    }
     st$points[c, ] <- pt
     e <- st$points[seq_len(c), 1]
     best <- oob_select(e)
@@ -161,7 +172,8 @@ oob_stopper <- function(y, ntrees, checkpoints, patience = 4L, tolerance = 0, fl
       splits = as.integer(checkpoints[st$sel]),
       path = oob_path_frame(checkpoints[g], st$points[g, , drop = FALSE], st$n_rows, crit),
       stopped_at = as.integer(checkpoints[st$c]),
-      capped = st$capped
+      capped = st$capped,
+      level = st$levels[[st$sel]]
     )
   }
   st

R/print.R

@@ -136,7 +136,21 @@ format.rpf <- function(x, ...) {
         if (identical(params$loss, "logit") && !is.null(params$exact_step)) "exact_step",
         if (params$loss %in% c("logit", "exponential") && !is.null(params$prior_rows)) "prior_rows"
       )
-      print_params(params[loss_params])
+      # the level (classif_level()), on the score scale; multiclass per class, centred as the scores are reported
+      level <- NULL
+      if (!is.null(x$level)) {
+        lev <- x$level
+        if (length(lev) > 1) {
+          if (identical(params$loss, "logit")) {
+            lev <- c(0, lev)
+          }
+          lev <- stats::setNames(lev - mean(lev), levels(x$blueprint$ptypes$outcomes[[1]]))
+          level <- list(level = paste0(names(lev), " ", sprintf("%+.3f", lev), collapse = ", "))
+        } else {
+          level <- list(level = sprintf("%+.3f", lev))
+        }
+      }
+      print_params(c(params[loss_params], level))
     }
 
     cli::cli_text("")