FastRiskScore: Fitting sparse integer risk scores via Autoresearch

Here is the fitted model, fit in 0.02 seconds:

The fitted model stores the score and its risk curve:

model.points_           # {"worst area <= 906.6": 5, ...}, the lines of the score
model.total_score(X)    # the total points of each row
model.risk(score)       # 1 / (1 + exp(-(model.scale_ * score + model.intercept_)))
model.train_loss_       # the mean training log loss of that risk

The points determine the total score of each row. The risk for a total comes from a logistic curve, risk(s) = 1 / (1 + exp(-(a·s + b))), with a = model.scale_ and b = model.intercept_; both are fit to the training data once the points are chosen. The points are chosen to make the log loss of this curve as small as possible (Sec 3), so the risks in the diagram are calibrated on the training data. RiskSLIM and FasterRisk show the same kind of table; they write b as an integer offset and 1/a as a multiplier.

The first six are the files distributed with RiskSLIM, which FasterRisk also uses; breastcancer and spambase keep real-valued columns, and mammo has one small-integer column. The other eight are binarized here: each numeric column becomes indicators x ≤ t at its deciles, each categorical column one indicator per level. Datasets above 12,500 rows are subsampled to 12,500, and every dataset is split 80/20 into training and test rows.

Every binary classification dataset of TabArena (OpenML suite tabarena-v0.1) except bank-marketing, heloc and diabetes, which the visible set already contains. They were binarized by the same rule as the visible set, with at most 40 columns kept (those with the highest mutual information with the label on the training split), and never used during development. Builder: evolve_slim/data/build_data.py.

The baselines came from two searches for methods that can be run without a commercial solver: one of packages for risk scores and sparse logistic regression, and a later one of 2023 to 2026 papers that cite FasterRisk, RiskSLIM, fastSparse, OKRidge, scoring lists or the certifier of Liu et al. (ICML 2025). Each "+ rounding" method gives a sparse real-valued solution with at most k features, which FasterRisk's rounding step turns into points. Penalized methods (riskscores, L0Learn, the L1 path, imodels' SLIMClassifier) have their penalty searched until at most k points are nonzero, and points are clipped to [−5, 5]. scikit-psl and okridge pin numpy < 2 and are vendored with small fixes (a numpy 2 incompatibility in scikit-psl, a call to a removed method in okridge), and scikit-psl gains a stop after k stages. AutoScore is only available in R, so its recipe (rank by random-forest importance, fit a logistic regression, divide by the smallest coefficient and round) is reimplemented on the binary features. RiskSLIM has no open-solver version, so its cutting-plane algorithm is reimplemented with HiGHS, re-solving the integer program each round and starting from rounded heuristic solutions as RiskSLIM does; it proves optimality on 6 of the 70 visible problems. L0Learn runs from its R package, since the Python package has no build past Python 3.10, and the R methods run in one persistent R process per worker so that R's start-up is not timed. OKGLM has no intercept, so a constant column is added. imodels' previous SLIMClassifier solves its integer program only with a commercial solver installed; without one, which is the default, it rounds an L2 logistic regression. Probabilistic scoring lists are killed at 3× the time limit on 10 visible and 51 hidden problems, and riskscores often passes the time limit, since it checks the clock only between fits; a solver without a score is charged the base rate. Left out: scorepyo (abandoned, pinned to 2022 packages), weight-of-evidence scorecards (optbinning, scorecardpy; not sparse integer points), GroupFasterRisk (the same as FasterRisk without feature groups), recent methods without public code (a GPU branch and bound for k-sparse GLMs, net-benefit risk scores) and those that optimise a different objective (boosted risk scores, SHAPoint, AgentScore).

Scoring. The benchmark rejects any answer whose points are not integers, fall outside [−5, 5] or use more than k features; no solver in the figure had such an answer. Given the points, the risk curve is fit by Newton's method on the two numbers a and b, a convex problem. Each solver process first fits a tiny warm-up problem, untimed, so numba compilation is not counted; FastRiskScore's first fit on a new machine takes about 2 minutes longer for that reason.

Best known losses. The excess loss is measured against the lowest loss any integer solver reached on the problem in any of these runs, including FasterRisk with every search width raised (a beam of 50 × 50, a pool of 200 and 100 multipliers) and a 20-minute limit, which takes about 12× its default time.

RiskSLIM. The free CPLEX Community Edition refuses problems above 1,000 variables or constraints, and RiskSLIM's cutting planes pass that limit on wider datasets; it returned no score on 25 of 70 visible problems and 102 of 135 hidden ones, and is charged the base rate there. With a full CPLEX licence it would do better, but within these time limits it still certifies optimality only on small problems, as the FasterRisk paper also found. RiskSLIM and SLIM run until their time limit on most problems, so their times sit near 60 s.

C. One run per solver at k = 5 on each of the 27 datasets at full size, 27 at a time, with a 10-minute limit (killed at 30 minutes). RiskSLIM stops at the CPLEX Community Edition's size limit on 22 datasets; SLIM passes its limit on 25, since HiGHS finishes its current step before stopping, and on the remaining 2 returns a proven optimum of its 0-1 loss.

Table 1. Mean training log loss, mean test AUC and geometric-mean fit time (s) of every solver on the three suites of Fig 1. Eight of the baselines (cutting planes, abess, fastSparse, OKRidge, the L1 path, scoring lists, AutoScore and unit weighting) were added after a search for baselines that can be run, and five more (L0Learn, OKGLM, skscope and two versions of riskscores) after a search of papers that cite them (details). The last row is not a risk score: it is FasterRisk's real-valued solution before rounding, shown as a reference.

The suites above split each numeric column at its 9 deciles. We also built versions of the visible and hidden suites that split each numeric column at its 99 percentiles, which gives the search many more features to choose from (990 instead of 90 on magic). The six RiskSLIM datasets were rebuilt from their raw OpenML sources for this. Six more loops ran on the fine visible suite with the same keep rule. Loops 1 and 2 both started from FastRiskScore; each later loop started from the best version of the one before. From loop 3 on, a change was also rejected if it lowered the mean test AUC by more than 0.001. The fine hidden suite (27 datasets × 5 values of k) was not used during the loops. Table 2 gives the best version of each loop that kept a change, on both.

Table 2. Mean training log loss, mean test AUC and geometric-mean fit time (s) on the fine suites (99 thresholds per numeric column). Run 2's version is the first release of FastRiskScore; each loop row is that loop's best version and its main change; the last row is the shipped solver with n_thresholds > 19.

Versions: loop 1 f31_fastsame, loop 2 g24_fastcol, loop 3 h8_fastbeam, loop 4 i30_final, loop 5 j18_poliso. Hidden-suite AUCs to four decimals: FasterRisk 0.7917, run 2's version 0.7924, loops 1 and 2 0.7927, loop 3 0.7930, loops 4 and 5 0.7932, shipped 0.7937. Times vary by about 5% between runs; the loop 1 and 2 hidden rows were run while another loop used the same machine.

• Loss. Every loop beats FastRiskScore on the fine visible suite, but on the fine hidden suite the gains are small: 0.0002 at most. Loop 1 has the lowest loss on both suites, at 6 to 15 times the hidden-suite time of the later loops.
• Feature handling. Loop 2 found that the 99 indicators of a numeric column are nested (x ≤ t for increasing t), and handled each column's indicators as one ordered variable. This made the search faster and let the beam pick the best threshold of a column instead of filling with neighbouring thresholds of one column.
• AUC. Loop 3 required each chosen indicator to cover at least √N training rows on each side. Without it, the search picks thresholds that fit a handful of training rows. Against loop 2, it raised visible AUC by 0.003 for a loss increase of 0.0002, and the AUC gain held on the hidden suite (0.7927 to 0.7930).
• Speed. Loop 4 reached about the same loss as loop 3 in half the time: 0.021 s on the hidden suite, 7× faster than FastRiskScore and 220× faster than FasterRisk.
• Saturation. Loop 5's only kept change, an exact polish of the final score, lowered the visible loss by 0.0002. Almost all of that came from two near-separable problems (ionosphere at k = 7, mushroom at k = 10), and it changed 2 of the 135 hidden problems. On the problems where the polish runs, no single change of one feature's points, and no swap of one feature, lowers the loss of the final score. A sixth loop (30 attempts: new starting points, tabu walks, pair moves) kept nothing; running its whole search 4 times longer lowered the mean loss by 0.00001.

The fine-track versions are 3× faster than run 2's version, but on the decile suites they lost 0.0014 in visible loss, mostly on spambase and ionosphere. Two causes accounted for most of it: loop 3's minimum support per feature, and beam rules that kept one threshold per variable and one child per set of variables, which removed good feature sets on decile data. A loop on the decile visible suite started from loop 5's version with the minimum support switched off and kept four changes: candidates whose real-valued coefficient rounds to 0 are ranked last in the beam (spambase, whose large-scale columns collapsed 3 of 5 final feature sets to 2 features), the search kernels are compiled before the first fit (about 4 s had been spent compiling inside a few fits), the beam diversity rules are switched off, and the beam is trimmed back to 8 parents. A second loop from that version found nothing that held on the hidden suite. The package ships the second loop's final version (n17_lean, with the same hidden-suite results), and switches to the fine-track settings when n_thresholds > 19. Table 3 compares it with run 2's version.

Table 3. Mean training log loss, mean test AUC and geometric-mean fit time (s) of the shipped version and run 2's version on the suites of Table 1.

AUCs to four decimals, visible / hidden / full size: run 2's version 0.8424 / 0.7937 / 0.7926, shipped 0.8421 / 0.7939 / 0.7920. Times vary by about 10% between runs.

On the hidden suite the shipped version has a lower loss than run 2's version on 23 of the 135 problems and a higher one on 10, and it is 3 to 4× faster on the three suites. Test AUC is within 0.0006 of run 2's version on every suite.

A scoring system (risk score, point system) is a linear model score = sum_j lambda_j x_j with few nonzero coefficients, each a small integer (typically in {-5..5}). A prediction is a threshold on the score (SLIM) or a logistic map of it, P(y=1) = 1/(1+exp(-(lambda_0 + score)/m)) (RiskSLIM, FasterRisk). Work falls into three groups: (1) exact mixed-integer programming (MIP) formulations that can certify optimality but need a commercial solver and scale to tens of features; (2) heuristics that fit a sparse continuous model and then round it (FasterRisk, AutoScore, rounding methods); (3) generic L0 best-subset solvers (L0Learn, abess, fastSparse, OKRidge) that produce the sparse continuous model a rounding step needs. All entries below were checked against the arXiv/publisher page and the code repository in September 2026.

Methods

Both runs found this change first, and it lowered the loss more than any other: on the visible set it moved the mean excess loss from about +0.002 to −0.001 (measured against the original best known losses) in each run. FasterRisk judges a rounding by its loss at a fixed multiplier, but a risk score is used with its fitted curve, so a rounding that is worse at the fixed multiplier can be better once a and b are fit again. On spambase at k = 3, FasterRisk's rounding sets one of its three features to zero, because at that feature's scale the best integer point rounds to 0. The local search puts the feature back and lowers the loss by 0.095.

At first the local search added 60 to 80% to the run time, because every candidate move required fitting a and b again over all rows. Two changes made it cheap. The loss is computed from the counts of positive and negative rows at each value of the total score (with 5 features the total takes at most 51 values), and moves are ranked by an estimate so that only the best two are checked exactly. Run 2 found that its first estimate was too optimistic for moves that change the spread of the scores (0.58 estimated against 0.63 true on compas at k = 3). It fixed this by recomputing the cells a move touches after one Newton step in (a, b).

A profile of FasterRisk on bank at k = 10 (19.9 s) showed 82% of the time in coordinate descent, which reads every row for each coordinate step. Merging identical rows and fitting by Newton's method on the few chosen features made both runs 5× faster. Run 2 then grouped each beam candidate's rows by their values on its features (another 2.5×), moved the creation of each new candidate into one compiled function together with smaller changes (−34% time), and computed the gradient and swap screens as single matrix products (−14%, then −12%). After these changes each step reads the data once instead of once per candidate.

The integer search also made parts of FasterRisk's first stage unnecessary. Run 2 removed FasterRisk's diverse pool (swapping each feature for up to 50 others) and started the local search from the 10 beam solutions instead, which gave a lower loss. Run 1 kept the pool but reduced it from 50 to 10 solutions and the rounding multipliers from 20 to 5, with no change in loss.

• Local search from an empty score, without the beam search: much worse (excess loss +0.010 against −0.001), so the features chosen by the beam search matter.
• Cheaper ways to rank beam candidates before fitting them (a one-step score test, a one-coordinate fit): all raised the loss by 0.0002 to 0.0007.
• Smaller beams: 5 parents instead of 10 raised the loss by 0.003.
• Random restarts and perturbations of the local search: large gains on a few problems where it got stuck (mushroom at k = 7), but less than 0.0001 on average, for 2 to 3 times the time. 300 random restarts per problem on 15 problems beat the final version on only 4.
• More starts, deeper swaps, rescaling all points at once, refitting the best set of features and rounding again at several resolutions: gains of at most 0.0002, for 20 to 60% more time, so none passed the keep rule.

Visible set, as recorded by each run. Times vary by about 5% between runs of the same code. The runs' logs (runs/sep26-run1/notes.md, runs/sep26-run2/notes.md) list every attempt, kept or not, with its measured effect.