Robust XGBoosting for Regression

# Robust XGBoosting for Regression
Source: [https://arxiv.org/html/2608.13590](https://arxiv.org/html/2608.13590)
###### Abstract

XGBoost is a very popular and powerful method for prediction\. It iteratively fits simple decision trees to the residuals of the previous step\. An efficient and scalable implementation is available\. The standard loss function for XGBoost is the quadratic loss, but a Huber loss can also be used\. In this paper, we study the robustness of XGBoost and show that its performance can be affected by vertical outliers and leverage points\. To address this, we explore alternative loss functions, based on M\-, S\-, andτ\\tau\-estimators from robust regression\. Our results indicate that a two\-step procedure, referred to as MM\-XGBoost, provides the best trade\-off between robustness and prediction accuracy\.

## 1Introduction

Machine learning has become an essential tool for predictive modeling\. Among the many available methods, gradient boosting\(Friedman,[2001](https://arxiv.org/html/2608.13590#bib.bib4)\)is a particularly powerful approach, as it builds a strong predictor by iteratively combining simple base learners\. At every iteration, a base learner is constructed by minimizing a loss function that measures the discrepancy between the model predictions and the observed response\. Extreme gradient boosting, or XGBoost\(Chen and Guestrin,[2016](https://arxiv.org/html/2608.13590#bib.bib3)\), is an efficient and scalable implementation of gradient boosting\. It is widely used due to its strong predictive performance, its built\-in regularization, and its flexibility in handling different types of data\.

In this paper, we show that the classical XGBoost for regression can be sensitive to outliers\. Vertical outliers in the response variable and leverage points in the predictor variables can have a strong influence on the fitted model, reducing the accuracy and stability of the prediction\. This sensitivity is mainly observed when the commonly used quadratic loss function of XGBoost is used\. To address this issue, the XGBoost implementation provides a more robust alternative in the form of a smooth approximation of the Huber loss\. However, we show that even with this Huber loss, XGBoost can still be affected by leverage points\.

Motivated by these limitations, we aim to further improve the robustness of XGBoost for regression tasks by modifying the loss function used in the boosting procedure\. In particular, we study several alternative loss functions, inspired by robust regression estimators\. Implementations of XGBoost\(Chenet al\.,[2025](https://arxiv.org/html/2608.13590#bib.bib7)\)allow users to define custom loss functions without modifying the core algorithm, making it possible to incorporate such robust loss functions into the XGBoost algorithm without much difficulty\.

Our approach is based on ideas fromFriedman \([2001](https://arxiv.org/html/2608.13590#bib.bib4)\)andJu and Salibián\-Barrera \([2021](https://arxiv.org/html/2608.13590#bib.bib5)\), where the incorporation of robust loss functions into gradient boosting has been studied\. In particular,Friedman \([2001](https://arxiv.org/html/2608.13590#bib.bib4)\)andLutzet al\.\([2008](https://arxiv.org/html/2608.13590#bib.bib6)\)study loss functions based on M\-estimators of regression using a robust scale estimate of the residuals\.Ju and Salibián\-Barrera \([2021](https://arxiv.org/html/2608.13590#bib.bib5)\)use more sophisticated regression estimators and introduce loss functions based on S\- and MM\-estimators\. We adapt and extend their ideas to the extreme gradient boosting method\. We also include theτ\\tau\-loss into the boosting procedure, which to the best of our knowledge was not considered before\.

Through simulation studies, we evaluate which of our proposed alternative loss functions improve robustness against outlying observations while still maintaining good predictive performance\. For the remainder of the paper, Section[2](https://arxiv.org/html/2608.13590#S2)introduces the XGBoost method, Section[3](https://arxiv.org/html/2608.13590#S3)introduces our proposed robust extensions of XGBoost, and Section[4](https://arxiv.org/html/2608.13590#S4)presents the simulation studies and discusses the results\. Finally, Section[5](https://arxiv.org/html/2608.13590#S5)concludes the paper\.

## 2XGBoost for Regression

In this section, we first outline the principle of gradient boosting, followed by an introduction to extreme gradient boosting and the choice of loss functions\. We adopt a vectorized notation for the training data\. The number of predictor variables is denoted bypp, and the number of observations bynn\.

### 2\.1Gradient Boosting

Consider a training dataset with covariate matrix𝐗=\(𝐱1,…,𝐱n\)⊤∈ℝn×p\\mathbf\{X\}=\(\\mathbf\{x\}\_\{1\},\\dots,\\mathbf\{x\}\_\{n\}\)^\{\\top\}\\in\\mathbb\{R\}^\{n\\times p\}and response vector𝐲=\(y1,…,yn\)⊤∈ℝn\\mathbf\{y\}=\(y\_\{1\},\\dots,y\_\{n\}\)^\{\\top\}\\in\\mathbb\{R\}^\{n\}\. In gradient boosting with regression trees, the prediction model is constructed in a stage\-wise manner by sequentially adding regression trees\. Let𝐲^\(0\)∈ℝn\\hat\{\\mathbf\{y\}\}^\{\(0\)\}\\in\\mathbb\{R\}^\{n\}denote the initial prediction, defined as a constant vector

𝐲^\(0\)=y\(0\)​𝟏n,\\hat\{\\mathbf\{y\}\}^\{\(0\)\}=y^\{\(0\)\}\\mathbf\{1\}\_\{n\},where𝟏n=\(1​…​1\)⊤∈ℝ\\mathbf\{1\}\_\{n\}=\(1\\dots 1\)^\{\\top\}\\in\\mathbb\{R\}\. Here,y\(0\)∈ℝy^\{\(0\)\}\\in\\mathbb\{R\}is an initial estimate of the response, typically chosen as the mean of𝐲\\mathbf\{y\}\.

At each iteration, a regression treef:ℝp→ℝf:\\mathbb\{R\}^\{p\}\\rightarrow\\mathbb\{R\}is fitted and added to the current model\. A regression tree is defined by a tree structure

q:ℝp→\{1,…,T\},q:\\mathbb\{R\}^\{p\}\\rightarrow\\\{1,\\dots,T\\\},which assigns each observation to one of theTTleave nodes, and by an associated vector of leaf weights

𝐰=\(w1,…,wT\)⊤\.\\mathbf\{w\}=\(w\_\{1\},\\dots,w\_\{T\}\)^\{\\top\}\.Accordingly, for an observation𝐱∈ℝp\\mathbf\{x\}\\in\\mathbb\{R\}^\{p\}, the tree prediction is given by the weight of the corresponding leaf node, i\.e\.,

f​\(𝐱\)=wq​\(𝐱\)\.f\(\\mathbf\{x\}\)=w\_\{q\(\\mathbf\{x\}\)\}\.Applying the tree to all observations in the training data yields the prediction vector

F​\(𝐗\)=\(f​\(𝐱1\),…,f​\(𝐱n\)\)⊤∈ℝn\.F\(\\mathbf\{X\}\)=\\big\(f\(\\mathbf\{x\}\_\{1\}\),\\dots,f\(\\mathbf\{x\}\_\{n\}\)\\big\)^\{\\top\}\\in\\mathbb\{R\}^\{n\}\.
Letf\(k\)f^\{\(k\)\}be the tree obtained in iterationkkwith corresponding structureq\(k\)q^\{\(k\)\}and associated weights𝐰\(k\)\\mathbf\{w\}^\{\(k\)\}\. Let𝐲^\(k\)∈ℝn\\hat\{\\mathbf\{y\}\}^\{\(k\)\}\\in\\mathbb\{R\}^\{n\}denote the vector of predictions of thekk\-th boosting iteration\. At iteration\(k\+1\)\(k\+1\), the next regression tree is selected to improve the current model with respect to a specified loss functionℒ\\mathcal\{L\}, which measures the discrepancy between the observed responses and the model predictions\. More specifically, the treef\(k\+1\)f^\{\(k\+1\)\}is constructed such that, when added to the current fit, it reduces the loss

ℒ​\(𝐲,𝐲^\(k\)\+F\(k\+1\)​\(𝐗\)\)\.\\mathcal\{L\}\\big\(\\mathbf\{y\},\\hat\{\\mathbf\{y\}\}^\{\(k\)\}\+F^\{\(k\+1\)\}\(\\mathbf\{X\}\)\\big\)\.\(1\)The resulting tree is added to the existing model, yielding

𝐲^\(k\+1\)=𝐲^\(k\)\+η​F\(k\+1\)​\(𝐗\),\\hat\{\\mathbf\{y\}\}^\{\(k\+1\)\}=\\hat\{\\mathbf\{y\}\}^\{\(k\)\}\+\\eta F^\{\(k\+1\)\}\(\\mathbf\{X\}\),\(2\)whereη\>0\\eta\>0is a shrinkage parameter used to prevent overfitting\. The final predictions, afterKKiterations, are given by

𝐲^\(K\)=𝐲^\(0\)\+η​∑k=1KF\(k\)​\(𝐗\)\.\\hat\{\\mathbf\{y\}\}^\{\(K\)\}=\\hat\{\\mathbf\{y\}\}^\{\(0\)\}\+\\eta\\sum\_\{k=1\}^\{K\}F^\{\(k\)\}\(\\mathbf\{X\}\)\.Thus, the learned predictor for a new observation𝐱∈ℝp\\mathbf\{x\}\\in\\mathbb\{R\}^\{p\}can be expressed as the ensemble of trees

y^​\(𝐱\)=y\(0\)\+η​∑k=1Kf\(k\)​\(𝐱\),\\hat\{y\}\(\\mathbf\{x\}\)=y^\{\(0\)\}\+\\eta\\sum\_\{k=1\}^\{K\}f^\{\(k\)\}\(\\mathbf\{x\}\),wheref\(k\)​\(𝐱\)=wq\(k\)​\(𝐱\)\(k\)f^\{\(k\)\}\(\\mathbf\{x\}\)=w^\{\(k\)\}\_\{q^\{\(k\)\}\(\\mathbf\{x\}\)\}\. We refer to Chapter 10 ofHastieet al\.\([2009](https://arxiv.org/html/2608.13590#bib.bib15)\)for more details on gradient boosting\.

### 2\.2XGBoost

XGBoost is an extension of gradient boosting that aims to efficiently minimize the loss function \([1](https://arxiv.org/html/2608.13590#S2.E1)\) by using a second\-order Taylor approximation\. It also includes an explicit regularization term to control model complexity\. At each boosting iteration, the method builds a regression tree based on this local approximation of the loss around the current predictions\.

More specifically, at iterationkk, the loss is expanded around the prediction vector from the previous iteration,𝐲^\(k\)\\hat\{\\mathbf\{y\}\}^\{\(k\)\}, and a regularization termΩ​\(F\)\\Omega\(F\)is added to penalize the complexity of the newly added treeFF\. Let𝐠\\mathbf\{g\}and𝐇\\mathbf\{H\}denote the gradient and Hessian of the loss function with respect to𝐲^\\hat\{\\mathbf\{y\}\}, evaluated at𝐲^\(k\)\\hat\{\\mathbf\{y\}\}^\{\(k\)\}:

𝐠:=∇𝐲^ℒ​\(𝐲,𝐲^\)\|𝐲^=𝐲^\(k\)∈ℝn,𝐇:=∇𝐲^2ℒ​\(𝐲,𝐲^\)\|𝐲^=𝐲^\(k\)∈ℝn×n\.\\mathbf\{g\}:=\\nabla\_\{\\hat\{\\mathbf\{y\}\}\}\\mathcal\{L\}\(\\mathbf\{y\},\\hat\{\\mathbf\{y\}\}\)\\Big\|\_\{\\hat\{\\mathbf\{y\}\}=\\hat\{\\mathbf\{y\}\}^\{\(k\)\}\}\\in\\mathbb\{R\}^\{n\},\\qquad\\mathbf\{H\}:=\\nabla\_\{\\hat\{\\mathbf\{y\}\}\}^\{2\}\\mathcal\{L\}\(\\mathbf\{y\},\\hat\{\\mathbf\{y\}\}\)\\Big\|\_\{\\hat\{\\mathbf\{y\}\}=\\hat\{\\mathbf\{y\}\}^\{\(k\)\}\}\\in\\mathbb\{R\}^\{n\\times n\}\.\(3\)The approximated loss function at iterationkkis then given by

ℒ:=ℒ​\(𝐲,𝐲^\(k\)\+F​\(𝐗\)\)≈𝐠⊤​F​\(𝐗\)\+12​F​\(𝐗\)⊤​𝐇​F​\(𝐗\)\+Ω​\(F\),\\mathcal\{L\}:=\\mathcal\{L\}\\big\(\\mathbf\{y\},\\hat\{\\mathbf\{y\}\}^\{\(k\)\}\+F\(\\mathbf\{X\}\)\\big\)\\approx\\mathbf\{g\}^\{\\top\}F\(\\mathbf\{X\}\)\+\\frac\{1\}\{2\}F\(\\mathbf\{X\}\)^\{\\top\}\\mathbf\{H\}F\(\\mathbf\{X\}\)\+\\Omega\\big\(F\\big\),\(4\)where the constant termℒ​\(𝐲,𝐲^\(k\)\)\\mathcal\{L\}\\big\(\\mathbf\{y\},\\hat\{\\mathbf\{y\}\}^\{\(k\)\}\\big\)is omitted, as it does not influence the optimization\. The regularization termΩ​\(F\)\\Omega\\big\(F\\big\)is given by

Ω​\(F\)=γ​T\+12​λ​𝐰⊤​𝐰,\\Omega\\big\(F\\big\)=\\gamma T\+\\frac\{1\}\{2\}\\lambda\\mathbf\{w\}^\{\\top\}\\mathbf\{w\},whereTTdenotes the number of leaves in the tree and𝐰\\mathbf\{w\}the corresponding leaf weights\. The parametersγ\\gammaandλ\\lambdacontrol the penalization ofTTand𝐰\\mathbf\{w\}, respectively\.

For a fixed tree structureqq, letS∈\{0,1\}n×TS\\in\\\{0,1\\\}^\{n\\times T\}denote the leaf\-assignment matrix, defined by

Si​j=\{1,if​q​\(𝐱i\)=j,0,otherwise\.S\_\{ij\}=\\begin\{cases\}1,&\\text\{if \}q\(\\mathbf\{x\}\_\{i\}\)=j,\\\\ 0,&\\text\{otherwise\}\.\\end\{cases\}Each row ofSScontains exactly one nonzero entry, indicating the leaf to which the corresponding observation of the training data is assigned\. The vector of tree predictions can be now be written as

F​\(𝐗\)=S​𝐰\.F\(\\mathbf\{X\}\)=S\\mathbf\{w\}\.Thus, once the structureqqis fixed, determining the tree reduces to finding the optimal the leaf weights𝐰\\mathbf\{w\}\.SubstitutingF​\(𝐗\)=S​𝐰F\(\\mathbf\{X\}\)=S\\mathbf\{w\}into the approximation \([4](https://arxiv.org/html/2608.13590#S2.E4)\) yields

ℒ≈𝐠⊤​S​𝐰\+12​\(S​𝐰\)⊤​𝐇​\(S​𝐰\)\+12​λ​𝐰⊤​𝐰\+γ​T\.\\mathcal\{L\}\\approx\\mathbf\{g\}^\{\\top\}S\\mathbf\{w\}\+\\frac\{1\}\{2\}\(S\\mathbf\{w\}\)^\{\\top\}\\mathbf\{H\}\(S\\mathbf\{w\}\)\+\\frac\{1\}\{2\}\\lambda\\mathbf\{w\}^\{\\top\}\\mathbf\{w\}\+\\gamma T\.Differentiating with respect to𝐰\\mathbf\{w\}and setting the derivative equal to zero gives the optimal leaf weights

𝐰∗=−\(S⊤​𝐇​S\+λ​I\)−1​S⊤​𝐠\.\\mathbf\{w\}^\{\*\}=\-\\big\(S^\{\\top\}\\mathbf\{H\}S\+\\lambda I\\big\)^\{\-1\}S^\{\\top\}\\mathbf\{g\}\.\(5\)Hence, for a given tree structure, the optimal leaf weights can be obtained directly from the first\- and second\-order derivatives of the loss function\. The tree structureqqis constructed greedily in XGBoost by recursively evaluating candidate splits using a split\-gain criterion derived from the regularized objective\. The process continues until no further split provides sufficient improvement or a predefined stopping criterion is reached\.

#### Built\-in loss functions\.

The original XGBoost implementation ofChen and Guestrin \([2016](https://arxiv.org/html/2608.13590#bib.bib3)\)assumes the loss to be additive over individual observations\. The loss can be expressed in terms of the residualsri=yi−y^ir\_\{i\}=y\_\{i\}\-\\hat\{y\}\_\{i\}, fori=1,…,ni=1,\\dots,n, as

ℒ​\(𝐲,𝐲^\)=∑i=1nρ​\(ri\),\\mathcal\{L\}\(\\mathbf\{y\},\\hat\{\\mathbf\{y\}\}\)=\\sum\_\{i=1\}^\{n\}\\rho\(r\_\{i\}\),\(6\)whereρ:ℝ→ℝ\+\\rho:\\mathbb\{R\}\\rightarrow\\mathbb\{R\}^\{\+\}is an even positive function\. Thisρ\\rho\-function determines how prediction errors are penalized and therefore controls robustness\.

The default choice in the XGBoost implementation\(Chenet al\.,[2025](https://arxiv.org/html/2608.13590#bib.bib7)\)forρ\\rhois the quadraticρ\\rho\-function

ρ​\(u\)=12​u2\.\\rho\(u\)=\\tfrac\{1\}\{2\}u^\{2\}\.Thisρ\\rho\-function is optimal under Gaussian noise, but is highly sensitive to outliers as it penalizes large residuals quadratically\. We will refer to XGBoost using the quadratic loss as XGB Quadratic\. In Section[4](https://arxiv.org/html/2608.13590#S4)we show that XGB Quadratic is sensitive to both vertical outliers and leverage points in the training data\.

To improve robustness, the XGBoost implementation provides an alternative choice, the Pseudo\-Huberρ\\rho\-function\(Charbonnieret al\.,[1994](https://arxiv.org/html/2608.13590#bib.bib12); Zhang,[1997](https://arxiv.org/html/2608.13590#bib.bib13)\)

ρ​\(u\)=b2​\(1\+\(ub\)2−1\),\\rho\(u\)=b^\{2\}\\big\(\\sqrt\{1\+\\big\(\\frac\{u\}\{b\}\\big\)^\{2\}\}\-1\\big\),which is a smooth \(twice differentiable\) alternative for the Huberρ\\rho\-function\(Huber,[1992](https://arxiv.org/html/2608.13590#bib.bib2)\)

ρ​\(u\)=\{12​u2,\|u\|≤b,b​\(\|u\|−b2\),\|u\|\>b\.\\rho\(u\)=\\begin\{cases\}\\frac\{1\}\{2\}u^\{2\},&\|u\|\\leq b,\\\\ b\\big\(\|u\|\-\\frac\{b\}\{2\}\\big\),&\|u\|\>b\.\\end\{cases\}\(7\)The tuning parameterb\>0b\>0controls the transition between its quadratic and linear behavior\. Small residuals are penalized quadratically, whereas large residuals are penalized only linearly to reduce their influence\. In the XGBoost implementation\(Chenet al\.,[2025](https://arxiv.org/html/2608.13590#bib.bib7)\), the default tuning constant of the Pseudo\-Huberρ\\rho\-function isb=1b=1\. We refer to XGBoost using the Pseudo\-Huber loss as XGB Pseudo\-Huber\.

Because XGB Pseudo\-Huber penalizes large residuals linearly, unlike XGB Quadratic, it is expected to be more robust to outliers\. In the simulation studies presented in Section[4](https://arxiv.org/html/2608.13590#S4), we confirm this improvement in the presence of vertical outliers\. However, the results also show that its performance is still highly sensitive to leverage points\. This motivates the need for incorporating alternative loss functions within the boosting framework that are more robust than the loss functions provided by the standard XGBoost implementation\. An alternative loss function, very popular in robust statistics, is the Tukey’s bisquare function\(Tukey and others,[1977](https://arxiv.org/html/2608.13590#bib.bib14)\)

ρ​\(u\)=\{c26​\[1−\(1−\(uc\)2\)3\],\|u\|≤c,c26,\|u\|\>c\.\\rho\(u\)=\\begin\{cases\}\\tfrac\{c^\{2\}\}\{6\}\\big\[1\-\\big\(1\-\\big\(\\dfrac\{u\}\{c\}\\big\)^\{2\}\\big\)^\{3\}\\big\],&\|u\|\\leq c,\\\\ \\frac\{c^\{2\}\}\{6\},&\|u\|\>c\.\\end\{cases\}\(8\)Thisρ\\rho\-function flattens out beyond the cutoffc\>0c\>0, such that residuals with absolute value larger thancchave constant influence on the fit\. The derivativeψ\\psiof the Tukey’s bisquareρ\\rho\-function is redescending to zero\. The tuning constantcccontrols the trade\-off between robustness and efficiency: a smaller value ofccyields higher robustness by downweighting moderate residuals, while larger values allow more observations to contribute fully to the model fit\. The selected constantccensures high efficiency \(of95%95\\%\) under Gaussian errors while still providing high robustness to outliers\(Maronnaet al\.,[2006](https://arxiv.org/html/2608.13590#bib.bib1)\)\. Figure[1](https://arxiv.org/html/2608.13590#S2.F1)plots the consideredρ\\rho\-functions\.

![Refer to caption](https://arxiv.org/html/2608.13590v1/x1.png)Figure 1:Comparison of differentρ\\rho\-functions: Quadratic, Pseudo\-Huber \(withb=1b=1\), Huber \(withb=1\.345b=1\.345\), and Tukey’s bisquare \(withc=4\.865c=4\.865\)\.

## 3Robust XGBoost

We propose several variants of XGBoost, each defined by a loss function inspired by well\-known robust regression estimators\. Specifically, we introduceM\-XGBoost, based on M\-estimators\(Huber,[1992](https://arxiv.org/html/2608.13590#bib.bib2)\),S\-XGBoost, based on S\-estimators\(Rousseeuw and Yohai,[1984](https://arxiv.org/html/2608.13590#bib.bib8)\), andτ2\\tau^\{2\}\-XGBoost, based onτ\\tau\-estimators\(Yohai and Zamar,[1988](https://arxiv.org/html/2608.13590#bib.bib10)\)\. In addition, we consider a two\-stage procedure,MM\-XGBoost, similar in spirit to MM\-estimators\(Yohai,[1987](https://arxiv.org/html/2608.13590#bib.bib9)\)\. The latter combines the high robustness of S\-estimators with the more efficient M\-estimators\. We refer to Chapter44and55ofMaronnaet al\.\([2006](https://arxiv.org/html/2608.13590#bib.bib1)\)for an overview of robust regression methods\.

The XGBoost method requires the gradient𝐠\\mathbf\{g\}and Hessian𝐇\\mathbf\{H\}in the second\-order Taylor expansion of the loss function \([4](https://arxiv.org/html/2608.13590#S2.E4)\)\. Note that for an additive loss function of the form \([1](https://arxiv.org/html/2608.13590#S2.E1)\), all cross\-partial derivatives vanish, yielding a diagonal Hessian matrix𝐇\\mathbf\{H\}\. Since this Hessian matrix𝐇\\mathbf\{H\}becomes intractably large fornnbig, we approximate it by its diagonal\. We denote the diagonal vector of𝐇\\mathbf\{H\}by

𝐡=\(∂2ℒ∂y^12,…,∂2ℒ∂y^n2\)⊤\.\\mathbf\{h\}=\\big\(\\frac\{\\partial^\{2\}\\mathcal\{L\}\}\{\\partial\\hat\{y\}\_\{1\}^\{2\}\},\\dots,\\frac\{\\partial^\{2\}\\mathcal\{L\}\}\{\\partial\\hat\{y\}\_\{n\}^\{2\}\}\\big\)^\{\\top\}\.
So, to implement our custom loss functions, we only need to derive the corresponding first\- and second\-order derivatives

gj=∂ℒ​\(𝐲,𝐲^\)∂y^jandhj=∂2ℒ​\(𝐲,𝐲^\)∂y^j2,g\_\{j\}=\\frac\{\\partial\\mathcal\{L\}\(\\mathbf\{y\},\\hat\{\\mathbf\{y\}\}\)\}\{\\partial\\hat\{y\}\_\{j\}\}\\quad\\text\{and\}\\quad h\_\{j\}=\\frac\{\\partial^\{2\}\\mathcal\{L\}\(\\mathbf\{y\},\\hat\{\\mathbf\{y\}\}\)\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\},\(9\)forj=1,…,n\.j=1,\\dots,n\.In the next subsections we explicitly derive these quantities for each of the proposed XGBoost methods\.

### 3\.1M\-XGBoost

The first method we consider is the M\-XGBoost method\. Its loss function is similar to \([6](https://arxiv.org/html/2608.13590#S2.E6)\), but uses scaled residuals\. The resulting M\-loss is defined as

ℒM​\(𝐲,𝐲^\)=∑i=1nρ​\(riσ0\)​σ02,\\mathcal\{L\}\_\{M\}\(\\mathbf\{y\},\\hat\{\\mathbf\{y\}\}\)=\\sum\_\{i=1\}^\{n\}\\rho\\big\(\\frac\{r\_\{i\}\}\{\\sigma\_\{0\}\}\\big\)\\sigma\_\{0\}^\{2\},\(10\)whereri=yi−y^ir\_\{i\}=y\_\{i\}\-\\hat\{y\}\_\{i\}denotes the residual of theii\-th observation, fori=1,…,ni=1,\\dots,n\. The scaleσ0\\sigma\_\{0\}is treated as a fixed robust scale estimate, computed once from the response variable before the boosting iterations begin\. For M\-XGBoost, we use the median absolute deviation \(MAD\) of the responses

σ0=1\.48⋅medi⁡\|yi−medj⁡yj\|,\\sigma\_\{0\}=1\.48\\cdot\\operatorname\{med\}\_\{i\}\\,\|y\_\{i\}\-\\operatorname\{med\}\_\{j\}\\,y\_\{j\}\|,a well\-known robust alternative to the standard deviation\(Maronnaet al\.,[2006](https://arxiv.org/html/2608.13590#bib.bib1)\)\. Theρ\\rho\-function in \([10](https://arxiv.org/html/2608.13590#S3.E10)\) is the Huber function \([7](https://arxiv.org/html/2608.13590#S2.E7)\) with tuning constantb=1\.345b=1\.345, discussed in Section[2\.2](https://arxiv.org/html/2608.13590#S2.SS2)\. Note that settingσ0=1\\sigma\_\{0\}=1reduces the M\-loss in \([10](https://arxiv.org/html/2608.13590#S3.E10)\) to the standard XGBoost loss in \([6](https://arxiv.org/html/2608.13590#S2.E6)\)\. Furthermore, choosing the quadraticρ\\rho\-function reduces the M\-loss to the classical squared error loss\.

To incorporate the M\-loss \([10](https://arxiv.org/html/2608.13590#S3.E10)\) into the XGBoost framework, we derive the corresponding gradient and diagonal Hessian entries in Equation \([9](https://arxiv.org/html/2608.13590#S3.E9)\)\. Since the scaleσ0\\sigma\_\{0\}is treated as fixed and computed in advance from the response variable, the M\-loss depends on𝐲^\\hat\{\\mathbf\{y\}\}only through the residuals\. Throughout the paper, we denote byψ\\psiandψ′\\psi^\{\\prime\}the first and second derivatives of theρ\\rho\-function\. We obtain the following result \(the proof can be found in[Appendix](https://arxiv.org/html/2608.13590#A1.SSx1)\)\.

###### Result 1\.

The components of the gradient𝐠\\mathbf\{g\}and diagonal Hessian𝐡\\mathbf\{h\}of the M\-loss \([10](https://arxiv.org/html/2608.13590#S3.E10)\) are

gj=−ψ​\(r~j,0\)​σ0​and​hj=−ψ′​\(r~j,0\)g\_\{j\}=\-\\psi\(\\tilde\{r\}\_\{j,0\}\)\\,\\sigma\_\{0\}\\ \\mbox\{\\ and \\ \}\\ h\_\{j\}=\-\\psi^\{\\prime\}\(\\tilde\{r\}\_\{j,0\}\)\(11\)wherer~j,0=rj/σ0\\tilde\{r\}\_\{j,0\}=\{\{r\}\_\{j\}\}/\{\\sigma\_\{0\}\}forj=1,…,nj=1,\\dots,n\.

### 3\.2S\-XGBoost

In the M\-loss, the scaleσ0\\sigma\_\{0\}is treated as known and does not depend on the regression fit\. Although such loss functions can provide some robustness against outliers, it is better to estimate the scale jointly with the model predictions\. To address this, we consider the S\-XGBoost method, where the loss function is

ℒS=12​n​σ2,\\mathcal\{L\}\_\{S\}=\\tfrac\{1\}\{2\}n\\sigma^\{2\},\(12\)withσ\\sigmaa robust estimator of the residual scale\. For S\-XGBoost,σ\\sigmais taken as an M\-estimator of scale centered at zero, implicitly defined through the M\-scale equation\(Huber,[1992](https://arxiv.org/html/2608.13590#bib.bib2)\)

1n​∑i=1nρ1​\(riσ\)=κ\.\\frac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\\rho\_\{1\}\\big\(\\frac\{r\_\{i\}\}\{\\sigma\}\\big\)=\\kappa\.\(13\)The constantκ\\kappais taken asκ=𝔼​\[ρ1​\(ε\)\]\\kappa=\\mathbb\{E\}\[\\rho\_\{1\}\(\\varepsilon\)\]forε∼𝒩​\(0,1\)\\varepsilon\\sim\\mathcal\{N\}\(0,1\)\(Maronnaet al\.,[2006](https://arxiv.org/html/2608.13590#bib.bib1)\)\. We takeρ1\\rho\_\{1\}as Tukey’s bisquareρ\\rho\-function \([8](https://arxiv.org/html/2608.13590#S2.E8)\)\. Choosing the tuning constantc=1\.548c=1\.548yields a 50% breakdown point and corresponds toκ=0\.2\\kappa=0\.2\.

To incorporate the S\-loss \([12](https://arxiv.org/html/2608.13590#S3.E12)\) into the XGBoost framework, we need to derive the corresponding gradient and diagonal Hessian entries in Equation \([9](https://arxiv.org/html/2608.13590#S3.E9)\)\. Theσ\\sigmain \([12](https://arxiv.org/html/2608.13590#S3.E12)\) is not fixed, but is implicitly defined through the M\-scale \([13](https://arxiv.org/html/2608.13590#S3.E13)\)\. Since it depends on the residualsri=yi−y^ir\_\{i\}=y\_\{i\}\-\\hat\{y\}\_\{i\}, it therefore depends on the prediction vector𝐲^\\hat\{\\mathbf\{y\}\}\. We obtain the following result \(the proof can be found in[Appendix](https://arxiv.org/html/2608.13590#A1.SSx2)\)\.

###### Result 2\.

The components of the gradient𝐠\\mathbf\{g\}and diagonal Hessian𝐡\\mathbf\{h\}of the S\-loss \([12](https://arxiv.org/html/2608.13590#S3.E12)\), forj=1,…,nj=1,\\dots,n, are

gj=−n​σ​ψ1​\(r~j\)A1g\_\{j\}=\-\\,n\\sigma\\,\\frac\{\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\}\{A\_\{1\}\}\(14\)and

hj=n​\(ψ1′​\(r~j\)A1−2​ψ1′​\(r~j\)​ψ1​\(r~j\)​r~jA12\+ψ12​\(r~j\)A12​\(B1A1−1\)\+ψ12​\(r~j\)A12\),h\_\{j\}=n\\big\(\\frac\{\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\}\{A\_\{1\}\}\-2\\,\\frac\{\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\\tilde\{r\}\_\{j\}\}\{A\_\{1\}^\{2\}\}\+\\frac\{\\psi\_\{1\}^\{2\}\(\\tilde\{r\}\_\{j\}\)\}\{A\_\{1\}^\{2\}\}\\big\(\\frac\{B\_\{1\}\}\{A\_\{1\}\}\-1\\big\)\+\\frac\{\\psi\_\{1\}^\{2\}\(\\tilde\{r\}\_\{j\}\)\}\{A\_\{1\}^\{2\}\}\\big\),\(15\)whereA1=∑i=1nr~i​ψ1​\(r~i\)A\_\{1\}=\\sum\_\{i=1\}^\{n\}\\tilde\{r\}\_\{i\}\\,\\psi\_\{1\}\(\\tilde\{r\}\_\{i\}\)andB1=∑i=1nr~i​\(ψ1​\(r~i\)\+r~i​ψ1′​\(r~i\)\)B\_\{1\}=\\sum\_\{i=1\}^\{n\}\\tilde\{r\}\_\{i\}\\big\(\\psi\_\{1\}\(\\tilde\{r\}\_\{i\}\)\+\\tilde\{r\}\_\{i\}\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{i\}\)\\big\), withr~i=ri/σ\\tilde\{r\}\_\{i\}=r\_\{i\}/\\sigmafor1≤i≤n1\\leq i\\leq n\.

Note that for the quadraticρ\\rho\-function withκ=1/2\\kappa=\{1\}/\{2\}, the S\-loss reduces to the classical squared\-error loss, sinceσ2=∑i=1nri2/n\\sigma^\{2\}=\\sum\_\{i=1\}^\{n\}r\_\{i\}^\{2\}/n\. Consequently,A1=∑i=1nr~i2A\_\{1\}=\\sum\_\{i=1\}^\{n\}\\tilde\{r\}\_\{i\}^\{2\}andB1=2​A1B\_\{1\}=2A\_\{1\}, which impliesgj=−rjg\_\{j\}=\-r\_\{j\}andhj=1h\_\{j\}=1forj=1,…,nj=1,\\dots,n\. These derivatives coincide with Result[1](https://arxiv.org/html/2608.13590#Thmresult1)\.

S\-estimators of scale offer strong robustness, but have reduced efficiency under clean data\. Therefore, S\-XGBoost is expected to achieve high robustness, but possibly larger values for the prediction errors, in particular when no outliers are present\. To address this issue, we extend the ideas ofτ\\tau\- and MM\-estimators\(Maronnaet al\.,[2006](https://arxiv.org/html/2608.13590#bib.bib1)\)to XGBoost, resulting inτ2\\tau^\{2\}\-XGBoost and MM\-XGBoost\. The idea is that the improved efficiency of the regression estimators will translate into higher prediction accuracy under clean and contaminated data\.

### 3\.3τ2\\tau^\{2\}\-XGBoost

Like S\-XGBoost,τ2\\tau^\{2\}\-XGBoost minimizes a robust residual scale, but replaces the M\-scale by the more efficientτ\\tau\-scale\. Theτ2\\tau^\{2\}\-loss function is defined as

ℒτ2=n2​κ​τ2,\\mathcal\{L\}\_\{\\tau^\{2\}\}=\\frac\{n\}\{2\\kappa\}\\,\\tau^\{2\},\(16\)where theτ\\tau\-estimator of scale is defined as

τ2=1n​∑i=1nρ​\(riσ\)​σ2,\\tau^\{2\}=\\frac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\\rho\\big\(\\frac\{\{r\}\_\{i\}\}\{\\sigma\}\\big\)\\sigma^\{2\},\(17\)with residualsri=yi−y^ir\_\{i\}=y\_\{i\}\-\\hat\{y\}\_\{i\}\(Yohai and Zamar,[1988](https://arxiv.org/html/2608.13590#bib.bib10)\)\. The scaleσ\\sigmais a M\-scale, defined in Equation \([13](https://arxiv.org/html/2608.13590#S3.E13)\), where the functionρ1\\rho\_\{1\}is the Tukey’s bisquareρ\\rho\-function withc1=1\.548c\_\{1\}=1\.548\. For theρ\\rho\-function in \([17](https://arxiv.org/html/2608.13590#S3.E17)\), we take the Tukey’s bisquareρ\\rho\-function, now with a lager value for the tuning constant, namelyc=4\.865c=4\.865, to ensure high efficiency under Gaussian errors while maintaining strong robustness\.

To incorporate theτ2\\tau^\{2\}\-loss \([16](https://arxiv.org/html/2608.13590#S3.E16)\) into the XGBoost framework, we derive the corresponding gradient and diagonal Hessian entries in Equation \([9](https://arxiv.org/html/2608.13590#S3.E9)\)\. Note thatτ2\\tau^\{2\}, as defined in \([17](https://arxiv.org/html/2608.13590#S3.E17)\), depends on the residuals both directly and through the scaleσ\\sigma\. We obtain the following result \(the proof can be found in[Appendix](https://arxiv.org/html/2608.13590#A1.SSx3)\)\.

###### Result 3\.

The components of the gradient𝐠\\mathbf\{g\}and diagonal Hessian𝐡\\mathbf\{h\}of theτ2\\tau^\{2\}\-loss \([16](https://arxiv.org/html/2608.13590#S3.E16)\), forj=1,…,nj=1,\\dots,n, are

gj=−ψ1​\(r~j\)​σ2​κ​\(2​CA1−AA1\)−ψ​\(r~j\)​σ2​κg\_\{j\}=\-\\frac\{\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\\sigma\}\{2\\kappa\}\\big\(\\frac\{2C\}\{A\_\{1\}\}\-\\frac\{A\}\{A\_\{1\}\}\\big\)\-\\frac\{\\psi\(\\tilde\{r\}\_\{j\}\)\\sigma\}\{2\\kappa\}\(18\)and

hj=12​κ\[\\displaystyle h\_\{j\}=\\frac\{1\}\{2\\kappa\}\\Big\[\(∂σ∂y^j\)2​\(2​C\+B−3​A\)\+\(∂2σ∂y^j2\)​\(2​σ​C−σ​A\)\\displaystyle\\big\(\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)^\{2\}\(2C\+B\-3A\)\+\\big\(\\frac\{\\partial^\{2\}\\sigma\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}\\big\)\\big\(2\\sigma C\-\\sigma A\\big\)\+\\displaystyle\+\(∂σ∂y^j\)\(2ψ′\(r~j\)r~j−2ψ\(r~j\)\)\+ψ′\(r~j\)\],\\displaystyle\\big\(\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)\\big\(2\\psi^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\tilde\{r\}\_\{j\}\-2\\psi\(\\tilde\{r\}\_\{j\}\)\\big\)\+\\psi^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\Big\],\(19\)whereA=∑i=1nr~i​ψ​\(r~i\),C=∑i=1nρ​\(r~i\)A=\\sum\_\{i=1\}^\{n\}\\tilde\{r\}\_\{i\}\\psi\(\\tilde\{r\}\_\{i\}\),C=\\sum\_\{i=1\}^\{n\}\\rho\(\\tilde\{r\}\_\{i\}\),A1=∑i=1nr~i​ψ1​\(r~i\)A\_\{1\}=\\sum\_\{i=1\}^\{n\}\\tilde\{r\}\_\{i\}\\psi\_\{1\}\(\\tilde\{r\}\_\{i\}\)andB=∑i=1nr~i​\(ψ​\(r~i\)\+r~i​ψ′​\(r~i\)\)B=\\sum\_\{i=1\}^\{n\}\\tilde\{r\}\_\{i\}\\big\(\\psi\(\\tilde\{r\}\_\{i\}\)\+\\tilde\{r\}\_\{i\}\\psi^\{\\prime\}\(\\tilde\{r\}\_\{i\}\)\\big\), withr~i=ri/σ\\tilde\{r\}\_\{i\}=r\_\{i\}/\\sigmafor1≤i≤n1\\leq i\\leq n\.

Forρ=ρ1\\rho=\\rho\_\{1\}with equal tuning constants, theτ2\\tau^\{2\}\-loss coincides with the S\-loss\. Thenψ=ψ1\\psi=\\psi\_\{1\},A=A1A=A\_\{1\}andC=n​κC=n\\kappa\. Thejj\-th gradient entry of theτ2\\tau^\{2\}\-loss, forj=1,…,nj=1,\\dots,n, simplifies togj=−n​σ​ψ​\(r~j\)/Ag\_\{j\}=\-n\\sigma\\,\{\\psi\(\\tilde\{r\}\_\{j\}\)\}/\{A\}, which is exactly the S\-loss gradient entry \([14](https://arxiv.org/html/2608.13590#S3.E14)\)\. An analogous simplification holds for the diagonal Hessian entry\.

### 3\.4MM\-XGBoost

In addition toτ\\tau\-estimators, MM\-estimators provide another approach to improve efficiency while maintaining strong robustness properties\(Yohai,[1987](https://arxiv.org/html/2608.13590#bib.bib9)\)\. We adapt this idea to the XGBoost framework and refer to the resulting method as MM\-XGBoost\. The method follows a two\-stage procedure that combines a highly robust initial fit with a refinement step\. In the first stage, we obtain robust predictions using S\-XGBoost\. In the second stage, a robust scale estimate is computed from the residuals of the initial fit\. The predictions from the first stage are used as the starting vector in the second stage\. In summary,

#### Stage 1: Robust Initial Fit\.

Apply S\-XGBoost \(Section[3\.2](https://arxiv.org/html/2608.13590#S3.SS2)\) withK1K\_\{1\}boosting iterations to obtain a highly robust initial predictor\. Let

rS,i=yi−y^S,i,i=1,…,n,r\_\{S,i\}=y\_\{i\}\-\\hat\{y\}\_\{S,i\},\\qquad i=1,\\dots,n,denote the residuals from the final S\-XGBoost prediction\. Define the scale parameterσ0\\sigma\_\{0\}as the M\-estimator of scale \([13](https://arxiv.org/html/2608.13590#S3.E13)\) computed from the residualsrS,ir\_\{S,i\}, withρ1\\rho\_\{1\}the Tukey’s bisquare function \([8](https://arxiv.org/html/2608.13590#S2.E8)\) with tuning constantc=1\.345c=1\.345\.

#### Stage 2: Refinement\.

Apply XGBoost with the M\-loss defined in Equation \([10](https://arxiv.org/html/2608.13590#S3.E10)\) withK2K\_\{2\}boosting iterations, using as initial predictions the final output from stage 1:

𝐲^\(0\)=\(y^S,1,…,y^S,n\)⊤\.\\hat\{\\mathbf\{y\}\}^\{\(0\)\}=\(\\hat\{y\}\_\{S,1\},\\dots,\\hat\{y\}\_\{S,n\}\)^\{\\top\}\.In this second stage, the M\-loss employs the Tukey’s bisquareρ\\rho\-function \([8](https://arxiv.org/html/2608.13590#S2.E8)\) with tuning constantc=4\.685c=4\.685, which ensures high efficiency under Gaussian errors while maintaining robustness\. The scaleσ0\\sigma\_\{0\}is the M\-estimator of scale obtained in stage 1\. We take the number of boosting iterationsK1K\_\{1\}in the first step equal to the number of boosting iterationsK2K\_\{2\}in the second step\.

## 4Simulation Study

To evaluate the robustness and predictive performance of the proposed XGBoost\-based methods introduced in Section[3](https://arxiv.org/html/2608.13590#S3), we conducted a simulation study\. Robustness is assessed by comparing prediction performance when models are trained under various contamination scenarios\. In particular, we compare the results obtained from the contaminated training data with those in the clean data case\. We consider the baseline methods XGB Quadratic and XGB Pseudo\-Huber together with the proposed XGBoost variants: M\-XGBoost, S\-XGBoost,τ2\\tau^\{2\}\-XGBoost, and MM\-XGBoost\. After presenting the simulation settings, we discuss the details of the model implementation\. Next, we present the performance measure used and we discuss the simulation results at length\.

### 4\.1Data Generation

Consider a training set\{\(𝐱i,yi\),i=1,…,n\}\\\{\(\\mathbf\{x\}\_\{i\},y\_\{i\}\),i=1,\\dots,n\\\}consisting of predictor variables𝐱i∈ℝp\\mathbf\{x\}\_\{i\}\\in\\mathbb\{R\}^\{p\}and responsesyi∈ℝy\_\{i\}\\in\\mathbb\{R\}\. The responses are generated from the regression model

yi=g​\(𝐱i\)\+εi,i=1,…,n,y\_\{i\}=g\(\\mathbf\{x\}\_\{i\}\)\+\\varepsilon\_\{i\},\\qquad i=1,\\dots,n,where the errorsεi\\varepsilon\_\{i\}arei\.i\.d\.​𝒩​\(0,1\)\\text\{i\.i\.d\. \}\\mathcal\{N\}\(0,1\)and independent of the predictor variables𝐱i\\mathbf\{x\}\_\{i\}\. We consider several regression functions forg​\(⋅\)g\(\\cdot\)\. Letp≥5p\\geq 5and write𝐱=\(x1,…,xp\)⊤\\mathbf\{x\}=\(x\_\{1\},\\dots,x\_\{p\}\)^\{\\top\}\. We consider a low\-dimensional setting withp=10p=10covariates, and a high\-dimensional setting withp=100p=100\. The first regression function is the linear functiong1​\(⋅\)g\_\{1\}\(\\cdot\):

g1​\(𝐱\)=∑j=15xj\.g\_\{1\}\(\\mathbf\{x\}\)=\\sum\_\{j=1\}^\{5\}x\_\{j\}\.In addition, we include two non\-linear regression functionsg2​\(⋅\)g\_\{2\}\(\\cdot\)andg3​\(⋅\)g\_\{3\}\(\\cdot\)adapted from the simulation settings inJu and Salibián\-Barrera \([2021](https://arxiv.org/html/2608.13590#bib.bib5)\):

g2​\(𝐱\)\\displaystyle g\_\{2\}\(\\mathbf\{x\}\)=5​x12\+\(x2​x3−1x2​x4\)2,\\displaystyle=5\\sqrt\{\{x\}\_\{1\}^\{2\}\+\\big\(\{x\}\_\{2\}\{x\}\_\{3\}\-\\frac\{1\}\{\{x\}\_\{2\}\{x\}\_\{4\}\}\\big\)^\{2\}\},g3​\(𝐱\)\\displaystyle g\_\{3\}\(\\mathbf\{x\}\)=\[∑j=15\(1\+\(−1\)j​0\.8​xj\+sin⁡\(6​xj\)\)\]​\[∑j=13\(1\+xj3\)\]\.\\displaystyle=\\big\[\\sum\_\{j=1\}^\{5\}\\big\(1\+\(\-1\)^\{j\}0\.8\{x\}\_\{j\}\+\\sin\(6\{x\}\_\{j\}\)\\big\)\\big\]\\big\[\\sum\_\{j=1\}^\{3\}\\big\(1\+\\frac\{\{x\}\_\{j\}\}\{3\}\\big\)\\big\]\.
The components of the predictor variables are generated from a uniform distribution𝒰​\(0,1\)\\mathcal\{U\}\(0,1\), except forg2g\_\{2\}, wherex2x\_\{2\}andx4x\_\{4\}follow a𝒰​\(1,2\)\\mathcal\{U\}\(1,2\)distribution due to constraints ong2g\_\{2\}\.

For the correlation structures of the predictor variables, we consider three scenario’s, again adapted from the simulation settings ofJu and Salibián\-Barrera \([2021](https://arxiv.org/html/2608.13590#bib.bib5)\)\. The first structure,S1S\_\{1\}, assumes independent components\. So the correlation matrix of𝐱\\mathbf\{x\},Cor​\(𝐱\)\\mathrm\{Cor\}\(\\mathbf\{x\}\), is diagonal\. The second structure,S2S\_\{2\}, introduces a decreasing correlation pattern

Cor​\(𝐱\)j​k=0\.8\|j−k\|\.\\mathrm\{Cor\}\(\\mathbf\{x\}\)\_\{jk\}=0\.8^\{\|j\-k\|\}\.This structure induces stronger dependence between predictors with nearby indices and weaker dependence as their indices are farther apart\. The third structure,S3S\_\{3\}, follows a block correlation pattern

Cor​\(𝐱\)j​k=\{0\.8if​xj​and​xk​belong to the same block,0otherwise\.\\mathrm\{Cor\}\(\\mathbf\{x\}\)\_\{jk\}=\\begin\{cases\}0\.8&\\text\{if \}x\_\{j\}\\text\{ and \}x\_\{k\}\\text\{ belong to the same block\},\\\\\[2\.84526pt\] 0&\\text\{otherwise\}\.\\end\{cases\}The composition of the blocks depends on the underlying regression function\. Forg1g\_\{1\}andg3g\_\{3\}, block 1 consists of\(x1,x2,x3\)\(x\_\{1\},x\_\{2\},x\_\{3\}\)and block 2 of\(x4,x5\)\(x\_\{4\},x\_\{5\}\)\. Forg2g\_\{2\}, the blocks are\(x1,x2\)\(x\_\{1\},x\_\{2\}\)and\(x3,x4\)\(x\_\{3\},x\_\{4\}\)\. This structure represents settings in which predictors are strongly correlated within groups but independent between groups\.

We consider two contamination rates in the training data: low contamination \(5%5\\%\) and high contamination \(20%20\\%\)\. Letℐout⊆\{1,…,n\}\\mathcal\{I\}\_\{\\text\{out\}\}\\subseteq\\\{1,\\dots,n\\\}denote a randomly selected index set of contaminated observations\. We consider the following three contamination scenarios along with the clean data case\.

#### \(a\) Vertical Outliers\.

The responses are shifted:

yi\(vertical\)=yi\+δiy,fori∈ℐout,y\_\{i\}^\{\(\\text\{vertical\}\)\}=y\_\{i\}\+\\delta\_\{i\}^\{y\},\\qquad\\text\{for\}\\quad i\\in\\mathcal\{I\}\_\{\\text\{out\}\},whereδiy∼t1\\delta\_\{i\}^\{y\}\\sim t\_\{1\}, a Student\-ttdistribution with 1 degree of freedom\. The predictor variables remain unchanged\.

#### \(b\) Good Leverage Points\.

The predictor variables are shifted:

xi​j\(good\)=xi​j\+δi​jx,fori∈ℐout,\{x\}\_\{ij\}^\{\(\\text\{good\}\)\}=\{x\}\_\{ij\}\+\\delta\_\{ij\}^\{x\},\\qquad\\text\{for\}\\quad i\\in\\mathcal\{I\}\_\{\\text\{out\}\},whereδi​jx∼𝒩​\(μi​j,10\)\\delta\_\{ij\}^\{x\}\\sim\\mathcal\{N\}\(\\mu\_\{ij\},10\)andμi​j∼𝒩​\(10,10\)\\mu\_\{ij\}\\sim\\mathcal\{N\}\(10,10\), forj=1,…,pj=1,\\dots,p\. The response values are then regenerated according to the model:

yi\(good\)=g​\(𝐱i\(good\)\)\+εi,fori∈ℐout\.y\_\{i\}^\{\(\\text\{good\}\)\}=g\\big\(\\mathbf\{x\}\_\{i\}^\{\(\\text\{good\}\)\}\\big\)\+\\varepsilon\_\{i\},\\qquad\\text\{for\}\\quad i\\in\\mathcal\{I\}\_\{\\text\{out\}\}\.

#### \(c\) Bad Leverage Points\.

First, predictor variables are shifted:

xi​j\(bad\)=xi​j\+δi​jx,fori∈ℐout,\{x\}\_\{ij\}^\{\(\\text\{bad\}\)\}=\{x\}\_\{ij\}\+\\delta\_\{ij\}^\{x\},\\qquad\\text\{for\}\\quad i\\in\\mathcal\{I\}\_\{\\text\{out\}\},whereδi​jx∼𝒩​\(μi​j,10\)\\delta\_\{ij\}^\{x\}\\sim\\mathcal\{N\}\(\\mu\_\{ij\},10\)andμi​j∼𝒩​\(10,10\)\\mu\_\{ij\}\\sim\\mathcal\{N\}\(10,10\), forj=1,…,pj=1,\\dots,p\. Then, the responses are regenerated and additionally shifted:

yi\(bad\)=g​\(𝐱i\(bad\)\)\+εi\+δiy,fori∈ℐout,y\_\{i\}^\{\(\\text\{bad\}\)\}=g\\big\(\\mathbf\{x\}\_\{i\}^\{\(\\text\{bad\}\)\}\\big\)\+\\varepsilon\_\{i\}\+\\delta\_\{i\}^\{y\},\\qquad\\text\{for\}\\quad i\\in\\mathcal\{I\}\_\{\\text\{out\}\},whereδiy∼t1\\delta\_\{i\}^\{y\}\\sim t\_\{1\}\.

### 4\.2Model Implementation

The loss functions for the baseline methods \(XGB Quadratic and XGB Pseudo\-Huber\) are built\-in options ofxgboostin R/Python\. In contrast, the other methods \(M\-XGBoost, S\-XGBoost,τ2\\tau^\{2\}\-XGBoost, and MM\-XGBoost\) use custom loss functions, which are supported byxgb\.train\. Such a custom loss function requires the vectors𝐲^\\hat\{\\mathbf\{y\}\}and𝐲\{\\mathbf\{y\}\}as inputs, and should return the gradient vector𝐠\\mathbf\{g\}and the diagonal Hessian vector𝐡\\mathbf\{h\}as output\. Hence, it is not necessary to compute the loss functionℒ\\mathcal\{L\}itself, as one only needs𝐠\\mathbf\{g\}and𝐡\\mathbf\{h\}for the optimal leaf weights in \([5](https://arxiv.org/html/2608.13590#S2.E5)\)\. Implementing the custom loss functions therefore becomes immediate, as one only needs to code the explicit expressions for𝐠\\mathbf\{g\}and𝐡\\mathbf\{h\}given in Results 1\-3 of Section[3](https://arxiv.org/html/2608.13590#S3)\.

The initial prediction used as a start for the boosting iterations is chosen as the median of the response values

y\(0\)=median⁡\(y1,…,yn\)\.\{y\}^\{\(0\)\}=\\operatorname\{median\}\(y\_\{1\},\\dots,y\_\{n\}\)\.This provides a more robust starting point than the commonly used mean, and we use it for all methods in this study\. As will be shown, only taking a robust starting value is not sufficient to obtain a fully robust version of XGBoost\. The number of boosting iterations is fixed atK=100K=100\. For larger values ofKKno further improvement in prediction accuracy was observed, for the simulation settings considered\. We preferred to work with a fixed choice for the number of iterations, to focus solely on the comparison of the different loss functions\. An alternative strategy is to selectKKusing a validation set, combined with an early stopping rule\. For MM\-XGBoost we setK1=50K\_\{1\}=50iterations in stage 1 andK2=50K\_\{2\}=50iterations in stage 2, resulting in the same total number of iterations as for the other XGBoost methods\. Numerical experiments indicate that choosingK2≥K1K\_\{2\}\\geq K\_\{1\}yields the best performance\. The first stage is designed to ensure robustness and should therefore not involve an excessive number of iterations, whereas the second stage focuses on improving accuracy and thus benefits from a larger number of iterations\.

A set of tuning parameters must be specified when applying XGBoost\. The maximum tree depth is fixed atd=1d=1, thereby restricting the base learners to decision stumps\.

The tree construction method is set toexact, which employs an exact greedy search for optimal split finding rather than a histogram based method\(Chen and Guestrin,[2016](https://arxiv.org/html/2608.13590#bib.bib3)\)\. All other hyperparameters of XGBoost in the R/Python implementation\(Chenet al\.,[2025](https://arxiv.org/html/2608.13590#bib.bib7)\)are kept at their default\. These include the learning rateη\\eta\(default value0\.30\.3\), the regularization parametersγ\\gammaandλ\\lambda\(default values0and11, respectively\), and the minimum sum of instance weights required in a leaf node \(default value11\)\.

### 4\.3Performance Measure

For each simulation run, we randomly partition the data into a training set, used to fit the model, and a test set, used exclusively to evaluate predictive performance\. Predictive performance is evaluated using the Root Mean Squared Error \(RMSE\)\. The test set is kept clean, which means that no vertical outliers nor leverage points are introduced into it\. Letℐtest\\mathcal\{I\}\_\{\\text\{test\}\}denote the index set of test observations andy^i\\hat\{y\}\_\{i\}the prediction for observationi∈ℐtesti\\in\\mathcal\{I\}\_\{\\text\{test\}\}\. The RMSE is defined as

RMSE=1\|ℐtest\|​∑i∈ℐtest\(yi−y^i\)2\.\\text\{RMSE\}=\\sqrt\{\\frac\{1\}\{\|\\mathcal\{I\}\_\{\\text\{test\}\}\|\}\\sum\_\{i\\in\\mathcal\{I\}\_\{\\text\{test\}\}\}\\big\(y\_\{i\}\-\\hat\{y\}\_\{i\}\\big\)^\{2\}\}\.
Each simulation scenario is replicated100100times to ensure stable performance estimates\. For every scenario, we report the mean RMSE across the replicates and its corresponding standard error \(SE\)\. An alternative performance measure is the Mean Absolute Error, yielding results that are fully in line with the ones presented\.

### 4\.4Results

For each combination of regression function, dimensionality, correlation structure, and contamination scenario, we evaluate the predictive performance of all proposed XGBoost variants as well as the baseline methods XGB Quadratic and XGB Pseudo\-Huber\. The training sample size is fixed atntrain=700n\_\{\\text\{train\}\}=700and the test sample size atntest=300n\_\{\\text\{test\}\}=300\. We discuss two settings in detail\. The tables of the RMSE for all other simulation settings can be found in the supplementary material\.

#### Setting 1: Non\-linear function, low\-dimensional, correlated predictors, low contamination rate\.

The first setting considers the non\-linear regression functiong2g\_\{2\}in a scenario withp=10p=10predictors following correlation structureS3S\_\{3\}, and a contamination rate of5%5\\%\. Figure[2](https://arxiv.org/html/2608.13590#S4.F2)presents the boxplots of thelog⁡\(RMSE\)\\log\(\\text\{RMSE\}\)over100100simulation replicates\. Table[1](https://arxiv.org/html/2608.13590#S4.T1)reports the mean RMSE values with standard errors\. As expected, XGB Quadratic is not robust\. It is affected by both vertical outliers and leverage points, resulting in increased prediction error and variability compared to the clean data case\. Applying XGBoost with the more robust Pseudo\-Huber loss function improves performance in the presence of vertical outliers, but remains sensitive to leverage points\. A similar pattern is observed for M\-XGBoost: Although it provides robustness against vertical outliers, it remains highly sensitive to leverage points\. In contrast, S\-XGBoost shows strong robustness across all contamination scenarios, with stable prediction errors\. However, this robustness comes at the cost of slightly reduced efficiency in the clean data case, as it yields higher RMSE values than the baseline methods\. Theτ2\\tau^\{2\}\-XGBoost method improves efficiency relative to S\-XGBoost, achieving lower RMSE in the clean data case while maintaining good robustness\. Finally, MM\-XGBoost balances robustness and efficiency most effectively\. It substantially improves prediction error compared to S\-XGBoost, more so thanτ2\\tau^\{2\}\-XGBoost, while preserving a high level of robustness\. Across all contamination scenarios, it achieves overall the lowest RMSE with low variability, making it the overall best\-performing method in this setting\.

Table 1:Mean RMSE \(SE\) forg2g\_\{2\}withp=10p=10under clean and 5% contamination scenarios withS3S\_\{3\}forM=100M=100independent runs![Refer to caption](https://arxiv.org/html/2608.13590v1/x2.png)Figure 2:Boxplots oflog⁡\(R​M​S​E\)\\log\(\{RMSE\}\)across 100 independent runs for simulation setting11, and for different contamination scenarios\.
#### Setting 2: Nonlinear function, high\-dimensional, correlated predictors, high contamination rate\.

The second setting represents a more challenging scenario, combining the nonlinear regression functiong3g\_\{3\}, strong predictor correlationS2S\_\{2\}, high dimensionality \(p=100p=100\), and a high contamination rate of20%20\\%\. Figure[3](https://arxiv.org/html/2608.13590#S4.F3)shows the correspondinglog⁡\(RMSE\)\\log\(\\text\{RMSE\}\)boxplots over100100simulation replicates, and Table[2](https://arxiv.org/html/2608.13590#S4.T2)reports the mean RMSE values with standard errors\. XGB Quadratic is unstable, with large increases in both prediction error and variability under vertical outliers as well as leverage points\. XGB Pseudo\-Huber improves robustness to vertical outliers, but the method remains sensitive to leverage points\. A similar lack of robustness is observed for M\-XGBoost, which is affected by leverage points\. In contrast, S\-XGBoost remains highly robust across all contamination types, maintaining stable RMSE values\. However, its robustness comes at the cost of efficiency: under clean data it yields noticeably higher RMSE than the more efficient XGB Quadratic\. Theτ2\\tau^\{2\}\-XGBoost method improves prediction error relative to S\-XGBoost, achieving lower RMSE under clean data\. However, this comes with a slight reduction in robustness, as its prediction error increases more noticeably under contamination compared to S\-XGBoost\. Moreover, it remains less efficient than XGB Quadratic in the clean setting\. Finally, MM\-XGBoost provides the best overall performance\. It substantially improves efficiency compared to S\-XGBoost while maintaining strong robustness across all contamination scenarios\. In particular, it reaches low prediction error under clean data and maintains stable performance in the presence of vertical outliers, leverage points, strong predictor correlation, high dimensionality, and substantial amounts of outliers\. These results demonstrate the ability of MM\-XGBoost to combine high efficiency with strong robustness in a more challenging setting\.

Table 2:Mean RMSE \(SE\) forg3g\_\{3\}withp=100p=100under clean and 20% contamination scenarios withS2S\_\{2\}forM=100M=100independent runs![Refer to caption](https://arxiv.org/html/2608.13590v1/x3.png)Figure 3:Boxplots oflog⁡\(R​M​S​E\)\\log\(\{RMSE\}\)across 100 independent runs for simulation setting22and for different contamination scenarios\.
#### Summary of Results Across All Configurations\.

The general patterns observed in the two representative settings extend across all simulation configurations\. From the tables in the supplementary file it is confirmed that \(i\) XGB Quadratic is sensitive to all types of contamination, while XGB Pseudo\-Huber improves robustness to vertical outliers but remains sensitive to leverage points\. M\-XGBoost shows a similar lack of robustness to leverage points\. \(ii\) S\-XGBoost provides strong robustness in all settings, with stable performance under both vertical outliers and leverage points\. However, this robustness is achieved at the cost of increased prediction error\. \(iii\) Theτ2\\tau^\{2\}\-XGBoost method improves efficiency relative to S\-XGBoost, but this gain becomes modest in presence of outliers\. \(iv\) Across all combinations of regression functions, contamination levels, and correlation structures, MM\-XGBoost achieves the best overall trade\-off between robustness and predictive accuracy\. In presence of contamination, it consistently outperforms all competing methods\. In absence of outliers, the loss in predictive performance compared to the standard XGBoost is limited\.

## 5Conclusion

This paper proposes highly robust versions of the widely used XGBoost prediction method for regression\. An important practical advantage of the proposed methods is that they are straightforward to implement within the XGBoost framework, which supports customized loss functions\. The user only needs to specify the corresponding gradient and diagonal Hessian of the loss function, which have been mathematically derived in Section[3](https://arxiv.org/html/2608.13590#S3)\.

The proposed methods were evaluated through extensive simulation studies encompassing a wide range of outlier scenarios\. In Section[4](https://arxiv.org/html/2608.13590#S4)we discussed two representative simulation settings in detail\. Results for other settings are included in the supplementary file, including a situation where the number of predictorsp=400p=400exceeds the size of the training setn=300n=300\. The results demonstrate that standard XGBoost is not robust to outliers, even when employing its built\-in Pseudo\-Huber loss function\. In contrast, the proposed MM\-XGBoost method achieves both a high degree of robustness and strong predictive accuracy, making it the most attractive approach overall\.

Our work is closely related to the robust gradient boosting methods proposed byJu and Salibián\-Barrera \([2021](https://arxiv.org/html/2608.13590#bib.bib5)\)\. Specifically, their methods, SBoost and RRBoost, can be viewed as counterparts of S\-XGBoost and MM\-XGBoost within the standard gradient boosting framework\. Unlike XGBoost, standard gradient boosting relies on a first\-order approximation of the loss function and does not require second\-order derivative information\. A key advantage of XGBoost is its computational efficiency and scalability, which largely explain its widespread popularity\. The robust methods proposed in this paper retain these advantages, as they modify only the loss function while preserving the core XGBoost algorithm\. Consequently, in our simulation settings, SBoost and RRBoost required up to one hundred times more computation time than S\-XGBoost and MM\-XGBoost, respectively\. In terms of predictive performance, however, the two families of methods yielded broadly comparable results, with MM\-XGBoost tending to slightly outperform RRBoost\. Hence, MM\-XGBoost offers a clear practical advantage by combining the robustness and predictive performance of RRBoost with the computational efficiency and scalability of the XGBoost framework\.

The robust loss functions considered in this paper can also be incorporated into the LightGBM framework\(Keet al\.,[2017](https://arxiv.org/html/2608.13590#bib.bib11)\)\. LightGBM is closely related to XGBoost and is often more computationally efficient, particularly on large datasets\. Our exploratory analyses indicate, however, that LightGBM is considerably more vulnerable to outliers than XGBoost, even when robust loss functions are employed\. This suggests that robustness may require modifications beyond the loss function alone, potentially involving the histogram\-based splitting mechanism\. We leave the development of robust LightGBM methods for future research\.

In conclusion, we proposed a robust extension of the widely used XGBoost regression method\. By providing an appropriate loss function, the proposed methods retain the computational efficiency and scalability of XGBoost while substantially improving robustness to outliers\. We hope that this work will stimulate further research on robust machine learning methods\.

## Appendix AppendixProofs

### Proof of Result[1](https://arxiv.org/html/2608.13590#Thmresult1)

The partial derivative of the M\-loss defined in Equation \([10](https://arxiv.org/html/2608.13590#S3.E10)\) with respect toy^j\\hat\{y\}\_\{j\}is given by

gj=∑i=1nψ​\(riσ0\)​\(−δi​jσ0\)​σ02,g\_\{j\}=\\sum\_\{i=1\}^\{n\}\\psi\\big\(\\frac\{r\_\{i\}\}\{\\sigma\_\{0\}\}\\big\)\\big\(\\frac\{\-\\delta\_\{ij\}\}\{\\sigma\_\{0\}\}\\big\)\\sigma\_\{0\}^\{2\},whereδi​j\\delta\_\{ij\}is the Kronecker delta\. Simplifying the expression gives the gradient entrygjg\_\{j\}of Equation \([11](https://arxiv.org/html/2608.13590#S3.E11)\) in Result[1](https://arxiv.org/html/2608.13590#Thmresult1)\. Differentiating the gradient entry with respect toy^j\\hat\{y\}\_\{j\}once more, yields the diagonal Hessian entryhjh\_\{j\}of Equation \([11](https://arxiv.org/html/2608.13590#S3.E11)\) in Result[1](https://arxiv.org/html/2608.13590#Thmresult1)\.

### Proof of Result[2](https://arxiv.org/html/2608.13590#Thmresult2)

The first derivative of the S\-loss defined in Equation \([12](https://arxiv.org/html/2608.13590#S3.E12)\) with respect toy^j\\hat\{y\}\_\{j\}, forj=1,…,nj=1,\\dots,n, is

gj=n​σ​∂σ∂y^j\.g\_\{j\}=n\\sigma\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\.\(20\)The M\-scale estimatorσ\\sigmais given by the solution of the implicit equation

1n​∑i=1nρ1​\(riσ\)=κ,\\frac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\\rho\_\{1\}\\big\(\\frac\{r\_\{i\}\}\{\\sigma\}\\big\)=\\kappa,\(21\)Defineψ1=ρ1′\\psi\_\{1\}=\\rho\_\{1\}^\{\\prime\}andr~i=ri/σ\\tilde\{r\}\_\{i\}=\{r\_\{i\}\}/\{\\sigma\}\. Differentiating \([21](https://arxiv.org/html/2608.13590#A1.E21)\) with respect toy^j\\hat\{y\}\_\{j\}, forj=1,…,nj=1,\\dots,n, and setting it equal to zero gives

0=∂∂y^j​\[1n​∑i=1nρ1​\(r~i\)\]=1n​∑i=1nψ1​\(r~i\)​∂r~i∂y^j\.0=\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\[\\frac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\\rho\_\{1\}\(\\tilde\{r\}\_\{i\}\)\\big\]=\\frac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\\psi\_\{1\}\(\\tilde\{r\}\_\{i\}\)\\frac\{\\partial\\tilde\{r\}\_\{i\}\}\{\\partial\\hat\{y\}\_\{j\}\}\.\(22\)The derivative ofr~i\\tilde\{r\}\_\{i\}with respect toy^j\\hat\{y\}\_\{j\}is

∂r~i∂y^j=∂∂y^j​\(riσ\)=−δi​j​σ−ri​∂σ∂y^jσ2=−δi​jσ−r~iσ​∂σ∂y^j\.\\frac\{\\partial\\tilde\{r\}\_\{i\}\}\{\\partial\\hat\{y\}\_\{j\}\}=\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\(\\frac\{r\_\{i\}\}\{\\sigma\}\\big\)=\\frac\{\-\\delta\_\{ij\}\\sigma\-r\_\{i\}\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\}\{\\sigma^\{2\}\}=\-\\frac\{\\delta\_\{ij\}\}\{\\sigma\}\-\\frac\{\\tilde\{r\}\_\{i\}\}\{\\sigma\}\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\.\(23\)Substituting this derivative into \([22](https://arxiv.org/html/2608.13590#A1.E22)\) and multiplying it byn​σn\\sigmaresults in

0=−ψ1​\(r~j\)−∂σ∂y^j​∑i=1nr~i​ψ1​\(r~i\)\.0=\-\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\-\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\sum\_\{i=1\}^\{n\}\\tilde\{r\}\_\{i\}\\psi\_\{1\}\(\\tilde\{r\}\_\{i\}\)\.By defining

A1=∑i=1nr~i​ψ1​\(r~i\),A\_\{1\}\\;=\\;\\sum\_\{i=1\}^\{n\}\\tilde\{r\}\_\{i\}\\psi\_\{1\}\(\\tilde\{r\}\_\{i\}\),\(24\)the first derivative ofσ\\sigmawith respect toy^j\\hat\{y\}\_\{j\}, forj=1,…,nj=1,\\dots,n, is now given by,

∂σ∂y^j=−ψ1​\(r~j\)A1\.\{\\quad\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\;=\\;\-\\frac\{\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\}\{A\_\{1\}\}\.\\quad\}\(25\)Substituting \([25](https://arxiv.org/html/2608.13590#A1.E25)\) into \([20](https://arxiv.org/html/2608.13590#A1.E20)\) yields the gradient entry \([14](https://arxiv.org/html/2608.13590#S3.E14)\) in Result[2](https://arxiv.org/html/2608.13590#Thmresult2)\. The diagonal Hessian entry is obtained by differentiating the gradient entry with respect toy^j\\hat\{y\}\_\{j\}once more:

hj=n​\(σ​∂2σ∂y^j2\+\(∂σ∂y^j\)2\)\.h\_\{j\}=n\\big\(\\sigma\\,\\frac\{\\partial^\{2\}\\sigma\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}\+\\Big\(\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\Big\)^\{2\}\\big\)\.\(26\)The second derivative of the M\-scale estimatorσ\\sigma, is obtained by differentiating the identity of the first derivative \([22](https://arxiv.org/html/2608.13590#A1.E22)\) again with respect toy^j\\hat\{y\}\_\{j\}forj=1,…,nj=1,\\dots,n:

∂∂y^j​ψ1​\(r~j\)\+∂2σ∂y^j2​A1\+∂σ∂y^j​∂A1∂y^j=0\.\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\+\\frac\{\\partial^\{2\}\\sigma\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}\\,A\_\{1\}\+\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\frac\{\\partial A\_\{1\}\}\{\\partial\\hat\{y\}\_\{j\}\}=0\.\(27\)We evaluate the first and third term of above equation\. For the first term, using the chain rule we obtain

∂∂y^j​ψ1​\(r~j\)=ψ1′​\(r~j\)​∂r~j∂y^j=ψ1′​\(r~j\)​\(−1σ−r~jσ​∂σ∂y^j\)\.\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)=\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\frac\{\\partial\\tilde\{r\}\_\{j\}\}\{\\partial\\hat\{y\}\_\{j\}\}=\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\big\(\-\\frac\{1\}\{\\sigma\}\-\\frac\{\\tilde\{r\}\_\{j\}\}\{\\sigma\}\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)\.For the third term, we compute∂A1/∂y^j\{\\partial A\_\{1\}\}/\{\\partial\\hat\{y\}\_\{j\}\},

∂A1∂y^j=∑i=1n∂∂y^j​\(r~i​ψ1​\(r~i\)\)\\displaystyle\\frac\{\\partial A\_\{1\}\}\{\\partial\\hat\{y\}\_\{j\}\}=\\sum\_\{i=1\}^\{n\}\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\(\\tilde\{r\}\_\{i\}\\psi\_\{1\}\(\\tilde\{r\}\_\{i\}\)\\big\)=∑i=1n\{∂r~i∂y^j​ψ1​\(r~i\)\+r~i​ψ1′​\(r~i\)​∂r~i∂y^j\}\\displaystyle=\\sum\_\{i=1\}^\{n\}\\big\\\{\\frac\{\\partial\\tilde\{r\}\_\{i\}\}\{\\partial\\hat\{y\}\_\{j\}\}\\psi\_\{1\}\(\\tilde\{r\}\_\{i\}\)\+\\tilde\{r\}\_\{i\}\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{i\}\)\\frac\{\\partial\\tilde\{r\}\_\{i\}\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\\\}=∑i=1n∂r~i∂y^j​\(ψ1​\(r~i\)\+r~i​ψ1′​\(r~i\)\)\.\\displaystyle=\\sum\_\{i=1\}^\{n\}\\frac\{\\partial\\tilde\{r\}\_\{i\}\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\(\\psi\_\{1\}\(\\tilde\{r\}\_\{i\}\)\+\\tilde\{r\}\_\{i\}\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{i\}\)\\big\)\.Define

B1=∑i=1nr~i​\(ψ1​\(r~i\)\+r~i​ψ1′​\(r~i\)\)\.B\_\{1\}\\;=\\;\\sum\_\{i=1\}^\{n\}\\tilde\{r\}\_\{i\}\\big\(\\psi\_\{1\}\(\\tilde\{r\}\_\{i\}\)\+\\tilde\{r\}\_\{i\}\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{i\}\)\\big\)\.\(28\)Substituting the derivative ofr~i\\tilde\{r\}\_\{i\}\([23](https://arxiv.org/html/2608.13590#A1.E23)\) into above equation yields

∂A1∂y^j=−1σ​\(ψ1​\(r~j\)\+r~j​ψ1′​\(r~j\)\)−∂σ∂y^j​B1σ\.\\frac\{\\partial A\_\{1\}\}\{\\partial\\hat\{y\}\_\{j\}\}=\-\\frac\{1\}\{\\sigma\}\\big\(\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\+\\tilde\{r\}\_\{j\}\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\big\)\-\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\frac\{B\_\{1\}\}\{\\sigma\}\.Inserting the evaluated first and third terms into \([27](https://arxiv.org/html/2608.13590#A1.E27)\) leads to

ψ1′​\(r~j\)​\(−1σ−r~jσ​∂σ∂y^j\)\+∂2σ∂y^j2​A1\+∂σ∂y^j​\[−1σ​\(ψ1​\(r~j\)\+r~j​ψ1′​\(r~j\)\)−∂σ∂y^j​B1σ\]=0\.\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\Big\(\-\\frac\{1\}\{\\sigma\}\-\\frac\{\\tilde\{r\}\_\{j\}\}\{\\sigma\}\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\Big\)\+\\frac\{\\partial^\{2\}\\sigma\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}\\,A\_\{1\}\+\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\Big\[\-\\frac\{1\}\{\\sigma\}\\big\(\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\+\\tilde\{r\}\_\{j\}\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\big\)\-\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\frac\{B\_\{1\}\}\{\\sigma\}\\Big\]=0\.Simplifying the above equation and writing in terms of∂2σ/∂y^j2\\partial^\{2\}\\sigma/\\partial\\hat\{y\}\_\{j\}^\{2\}gives

σ​A1​∂2σ∂y^j2=ψ1′​\(r~j\)\+ψ1′​\(r~j\)​r~j​∂σ∂y^j\+∂σ∂y^j​\(ψ1​\(r~j\)\+r~j​ψ1′​\(r~j\)\)\+B1σ​\(∂σ∂y^j\)2\.\\sigma A\_\{1\}\\frac\{\\partial^\{2\}\\sigma\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}=\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\+\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\tilde\{r\}\_\{j\}\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\+\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\(\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\+\\tilde\{r\}\_\{j\}\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\big\)\+\\frac\{B\_\{1\}\}\{\\sigma\}\\big\(\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)^\{2\}\.Finally, we substitute∂σ/∂y^j=−ψ1​\(r~j\)/A1\\partial\\sigma/\\partial\\hat\{y\}\_\{j\}=\-\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)/A\_\{1\}, obtained in \([25](https://arxiv.org/html/2608.13590#A1.E25)\), in the above equation\. After reordering the terms, the second derivative ofσ\\sigmawith respect toy^j\\hat\{y\}\_\{j\}, forj=1,…,nj=1,\\dots,n, is

∂2σ∂y^j2=ψ1′​\(r~j\)σ​A1−2​ψ1′​\(r~j\)​ψ1​\(r~j\)​r~jσ​A12\+ψ12​\(r~j\)σ​A12​\(B1A1−1\)\{\\quad\\frac\{\\partial^\{2\}\\sigma\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}=\\frac\{\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\}\{\\sigma A\_\{1\}\}\-2\\,\\frac\{\\psi\_\{1\}^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\\tilde\{r\}\_\{j\}\}\{\\sigma A\_\{1\}^\{2\}\}\+\\frac\{\\psi\_\{1\}^\{2\}\(\\tilde\{r\}\_\{j\}\)\}\{\\sigma A\_\{1\}^\{2\}\}\\big\(\\frac\{B\_\{1\}\}\{A\_\{1\}\}\-1\\big\)\\quad\}\(29\)whereA1A\_\{1\}andB1B\_\{1\}are defined in \([24](https://arxiv.org/html/2608.13590#A1.E24)\) and \([28](https://arxiv.org/html/2608.13590#A1.E28)\) respectively\. Substituting \([25](https://arxiv.org/html/2608.13590#A1.E25)\) and \([29](https://arxiv.org/html/2608.13590#A1.E29)\) into \([26](https://arxiv.org/html/2608.13590#A1.E26)\) yields the diagonal Hessian entry \([15](https://arxiv.org/html/2608.13590#S3.E15)\) in Result[2](https://arxiv.org/html/2608.13590#Thmresult2)\.

### Proof of Result[3](https://arxiv.org/html/2608.13590#Thmresult3)

The first derivative of theτ2\\tau^\{2\}\-loss derived in Equation \([16](https://arxiv.org/html/2608.13590#S3.E16)\) with respect toy^j\\hat\{y\}\_\{j\}, forj=1,…,nj=1,\\dots,n, follows directly from the chain rule\. Since theτ2\\tau^\{2\}\-loss depends ony^j\\hat\{y\}\_\{j\}throughτ2\\tau^\{2\}, we obtain

gj=n2​κ​∂τ2∂y^j\.g\_\{j\}=\\frac\{n\}\{2\\kappa\}\\,\\frac\{\\partial\\tau^\{2\}\}\{\\partial\\hat\{y\}\_\{j\}\}\.\(30\)Theτ\\tau\-scale estimator is given by

τ2=1n​∑i=1nρ​\(riσ\)​σ2,\\tau^\{2\}=\\frac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\\rho\\big\(\\frac\{r\_\{i\}\}\{\\sigma\}\\big\)\\sigma^\{2\},\(31\)withσ\\sigmaan M\-scale estimator which is the solution to the equation

1n​∑i=1nρ1​\(riσ\)=κ\.\\frac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\\rho\_\{1\}\\big\(\\frac\{r\_\{i\}\}\{\\sigma\}\\big\)=\\kappa\.To obtain the derivative of theτ\\tau\-scale estimator we differentiate \([31](https://arxiv.org/html/2608.13590#A1.E31)\) with respect toy^j\\hat\{y\}\_\{j\}forj=1,…​nj=1,\\dots n\. Note thatτ2\\tau^\{2\}depends ony^j\\hat\{y\}\_\{j\}both directly and throughσ\\sigma\. This results in

∂∂y^j​\[1n​∑i=1nρ​\(r~i\)​σ2\]=1n​∑i=1nψ​\(r~i\)​∂r~i∂y^j​σ2\+1n​∑i=1nρ​\(r~i\)​2​σ​∂σ∂y^j,\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\[\\frac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\\rho\(\\tilde\{r\}\_\{i\}\)\\sigma^\{2\}\\big\]=\\frac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\\psi\(\\tilde\{r\}\_\{i\}\)\\frac\{\\partial\\tilde\{r\}\_\{i\}\}\{\\partial\\hat\{y\}\_\{j\}\}\\sigma^\{2\}\+\\frac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\\rho\(\\tilde\{r\}\_\{i\}\)2\\sigma\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\},\(32\)whereψ=ρ′\\psi=\\rho^\{\\prime\}andr~i=ri/σ\\tilde\{r\}\_\{i\}=\{r\_\{i\}\}/\{\\sigma\}\. Define

C=∑i=1nρ​\(r~i\),andA=∑i=1nr~i​ψ​\(r~i\)\.C=\\sum\_\{i=1\}^\{n\}\\rho\(\\tilde\{r\}\_\{i\}\),\\qquad\\text\{and\}\\qquad A=\\sum\_\{i=1\}^\{n\}\\tilde\{r\}\_\{i\}\\psi\(\\tilde\{r\}\_\{i\}\)\.\(33\)Substituting the derivative ofr~i\\tilde\{r\}\_\{i\}\([23](https://arxiv.org/html/2608.13590#A1.E23)\) and the derivative ofσ\\sigma\([25](https://arxiv.org/html/2608.13590#A1.E25)\) into \([32](https://arxiv.org/html/2608.13590#A1.E32)\) gives

∂τ2∂y^j=−ψ1​\(r~j\)A1​\(−1n​A​σ\+2​σ​1n​C\)−ψ​\(r~j\)​σn,\\frac\{\\partial\\tau^\{2\}\}\{\\partial\\hat\{y\}\_\{j\}\}=\-\\frac\{\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\}\{A\_\{1\}\}\\big\(\\frac\{\-1\}\{n\}A\\sigma\+2\\sigma\\frac\{1\}\{n\}C\\big\)\-\\frac\{\\psi\(\\tilde\{r\}\_\{j\}\)\\sigma\}\{n\},\(34\)withA1A\_\{1\}defined in \([24](https://arxiv.org/html/2608.13590#A1.E24)\)\. After reordering Equation \([34](https://arxiv.org/html/2608.13590#A1.E34)\) we get the first derivative of theτ\\tau\-scale estimator with respect toy^j\\hat\{y\}\_\{j\}, forj=1,…,nj=1,\\dots,n,

∂τ2∂y^j=−ψ1​\(r~j\)​σn​\(2​CA1−AA1\)−ψ​\(r~j\)​σn\.\{\\quad\\frac\{\\partial\\tau^\{2\}\}\{\\partial\\hat\{y\}\_\{j\}\}=\-\\frac\{\\psi\_\{1\}\(\\tilde\{r\}\_\{j\}\)\\sigma\}\{n\}\\big\(\\frac\{2C\}\{A\_\{1\}\}\-\\frac\{A\}\{A\_\{1\}\}\\big\)\-\\frac\{\\psi\(\\tilde\{r\}\_\{j\}\)\\sigma\}\{n\}\.\}\(35\)Substituting \([35](https://arxiv.org/html/2608.13590#A1.E35)\) into \([30](https://arxiv.org/html/2608.13590#A1.E30)\) yields the gradient entry \([18](https://arxiv.org/html/2608.13590#S3.E18)\) in Result[3](https://arxiv.org/html/2608.13590#Thmresult3)\. The diagonal Hessian entry is obtained by differentiating the gradient entry with respect toy^j\\hat\{y\}\_\{j\}once more,

hj=n2​κ​∂2τ2∂y^j2\.h\_\{j\}=\\frac\{n\}\{2\\kappa\}\\frac\{\\partial^\{2\}\\tau^\{2\}\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}\.To obtain the second derivative of theτ\\tau\-scale estimatorτ2\\tau^\{2\}, we differentiate the identity of the first derivative \([32](https://arxiv.org/html/2608.13590#A1.E32)\) once again with respect toy^j\\hat\{y\}\_\{j\}forj=1,…,nj=1,\\dots,n,

∂2τ2∂y^j2=1n​\[∂∂y^j​\(−σ​A​∂σ∂y^j\)\+∂∂y^j​\(2​σ​C​∂σ∂y^j\)−∂∂y^j​\(σ​ψ​\(r~j\)\)\]\.\\displaystyle\\frac\{\\partial^\{2\}\\tau^\{2\}\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}=\\frac\{1\}\{n\}\\big\[\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\(\-\\sigma A\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)\+\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\(2\\sigma C\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)\-\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\(\\sigma\\psi\(\\tilde\{r\}\_\{j\}\)\)\\big\]\.\(36\)We evaluate the terms of the above equation separately\. Define

B=∑i=1nr~i​\(ψ​\(r~i\)\+r~i​ψ′​\(r~i\)\)\.B\\;=\\;\\sum\_\{i=1\}^\{n\}\\tilde\{r\}\_\{i\}\\big\(\\psi\(\\tilde\{r\}\_\{i\}\)\+\\tilde\{r\}\_\{i\}\\psi^\{\\prime\}\(\\tilde\{r\}\_\{i\}\)\\big\)\.Expanding and simplifying the first term in \([36](https://arxiv.org/html/2608.13590#A1.E36)\) gives

∂∂y^j​\(−σ​A​∂σ∂y^j\)\\displaystyle\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\(\-\\sigma A\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)=−A​\(∂σ∂y^j\)2−σ​\(−v′​\(r~j\)σ−1σ​∂σ∂y^j​B\)​∂σ∂y^j−σ​A​∂2σ∂y^j2\\displaystyle=\-A\\big\(\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)^\{2\}\-\\sigma\\big\(\\frac\{\-v^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\}\{\\sigma\}\-\\frac\{1\}\{\\sigma\}\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}B\\big\)\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\-\\sigma A\\frac\{\\partial^\{2\}\\sigma\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}=\(∂σ∂y^j\)2​\(B−A\)−σ​A​∂2σ∂y^j2\+v′​\(r~j\)​∂σ∂y^j\.\\displaystyle=\\big\(\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)^\{2\}\(B\-A\)\-\\sigma A\\frac\{\\partial^\{2\}\\sigma\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}\+v^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\.Expanding the second term in \([36](https://arxiv.org/html/2608.13590#A1.E36)\) gives the equation

∂∂y^j​\(2​σ​C​∂σ∂y^j\)=2​C​\(∂σ∂y^j\)2\+2​σ​\(∂C∂y^j​∂σ∂y^j\+C​∂2σ∂y^j2\)\.\\displaystyle\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\(2\\sigma C\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)=2C\\big\(\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)^\{2\}\+2\\sigma\\big\(\\frac\{\\partial C\}\{\\partial\\hat\{y\}\_\{j\}\}\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\+C\\frac\{\\partial^\{2\}\\sigma\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}\\big\)\.\(37\)To evaluate this second term we compute∂C/∂y^j\{\\partial C\}/\{\\partial\\hat\{y\}\_\{j\}\}, withCCdefined in \([33](https://arxiv.org/html/2608.13590#A1.E33)\)\. Using the derivative ofr~i\\tilde\{r\}\_\{i\}obtained in \([23](https://arxiv.org/html/2608.13590#A1.E23)\) we get

∂C∂y^j=∂∂y^j​∑i=1nρ​\(r~i\)=∑i=1nψ​\(r~i\)​∂r~i∂y^j=−Aσ​∂σ∂y^j−ψ​\(r~j\)σ\.\\displaystyle\\frac\{\\partial C\}\{\\partial\\hat\{y\}\_\{j\}\}=\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\sum\_\{i=1\}^\{n\}\\rho\(\\tilde\{r\}\_\{i\}\)=\\sum\_\{i=1\}^\{n\}\\psi\(\\tilde\{r\}\_\{i\}\)\\frac\{\\partial\\tilde\{r\}\_\{i\}\}\{\\partial\\hat\{y\}\_\{j\}\}=\\frac\{\-A\}\{\\sigma\}\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\-\\frac\{\\psi\(\\tilde\{r\}\_\{j\}\)\}\{\\sigma\}\.Plugging the derivative ofCCinto the second term \([37](https://arxiv.org/html/2608.13590#A1.E37)\) and rearranging the terms gives now the second term

∂∂y^j​\(2​σ​C​∂σ∂y^j\)=2​\(C−A\)​\(∂σ∂y^j\)2\+2​σ​C​∂2σ∂y^j2−2​ψ​\(r~j\)​∂σ∂y^j\.\\displaystyle\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\(2\\sigma C\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)=2\(C\-A\)\\big\(\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)^\{2\}\+2\\sigma C\\frac\{\\partial^\{2\}\\sigma\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}\-2\\psi\(\\tilde\{r\}\_\{j\}\)\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\.The third term of the second derivative \([36](https://arxiv.org/html/2608.13590#A1.E36)\) is

∂∂y^j​\(−σ​ψ​\(r~j\)\)=−ψ​\(r~j\)​∂σ∂y^j−σ​ψ′​\(r~j\)​∂r~j∂y^j\.\\displaystyle\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\(\-\\sigma\\psi\(\\tilde\{r\}\_\{j\}\)\\big\)=\-\\psi\(\\tilde\{r\}\_\{j\}\)\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\-\\sigma\\psi^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\frac\{\\partial\\tilde\{r\}\_\{j\}\}\{\\partial\\hat\{y\}\_\{j\}\}\.Substituting the derivative ofr~i\\tilde\{r\}\_\{i\}obtained in \([23](https://arxiv.org/html/2608.13590#A1.E23)\) gives

∂∂y^j​\(−σ​ψ​\(r~j\)\)=∂σ∂y^j​\(ψ′​\(r~j\)​r~j−ψ​\(r~j\)\)−ψ′​\(r~j\)\.\\displaystyle\\frac\{\\partial\}\{\\partial\\hat\{y\}\_\{j\}\}\(\-\\sigma\\psi\(\\tilde\{r\}\_\{j\}\)\)=\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\(\\psi^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\tilde\{r\}\_\{j\}\-\\psi\(\\tilde\{r\}\_\{j\}\)\\big\)\-\\psi^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\.Finally, we insert all evaluated terms into the second derivative equation of theτ\\tau\-scale estimator \([36](https://arxiv.org/html/2608.13590#A1.E36)\) w\.r\.ty^j\\hat\{y\}\_\{j\}forj=1,…,nj=1,\\dots,n, which results in

∂2τ2∂y^j2=1n\[\\displaystyle\\frac\{\\partial^\{2\}\\tau^\{2\}\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}=\\frac\{1\}\{n\}\\Big\[\(∂σ∂y^j\)2​\(2​C\+B−3​A\)\\displaystyle\\big\(\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)^\{2\}\(2C\+B\-3A\)\+\\displaystyle\+\(∂2σ∂y^j2\)​\(2​σ​C−σ​A\)\\displaystyle\\big\(\\frac\{\\partial^\{2\}\\sigma\}\{\\partial\\hat\{y\}\_\{j\}^\{2\}\}\\big\)\\big\(2\\sigma C\-\\sigma A\\big\)\+\\displaystyle\+\(∂σ∂y^j\)\(2ψ′\(r~j\)r~j−2ψ\(r~j\)\)\+ψ′\(r~j\)\]\.\\displaystyle\\big\(\\frac\{\\partial\\sigma\}\{\\partial\\hat\{y\}\_\{j\}\}\\big\)\\big\(2\\psi^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\tilde\{r\}\_\{j\}\-2\\psi\(\\tilde\{r\}\_\{j\}\)\\big\)\+\\psi^\{\\prime\}\(\\tilde\{r\}\_\{j\}\)\\Big\]\.\(38\)Substituting \([25](https://arxiv.org/html/2608.13590#A1.E25)\) and \([29](https://arxiv.org/html/2608.13590#A1.E29)\) into the second derivative ofτ2\\tau^\{2\}yields the diagonal Hessian entry \([3](https://arxiv.org/html/2608.13590#S3.Ex17)\) in Result[3](https://arxiv.org/html/2608.13590#Thmresult3)\.

## References

- P\. Charbonnier, L\. Blanc\-Feraud, G\. Aubert, and M\. Barlaud \(1994\)Two deterministic half\-quadratic regularization algorithms for computed imaging\.InProceedings of 1st international conference on image processing,Vol\.2,pp\. 168–172\.Cited by:[§2\.2](https://arxiv.org/html/2608.13590#S2.SS2.SSS0.Px1.p3.1)\.
- T\. Chen and C\. Guestrin \(2016\)Xgboost: a scalable tree boosting system\.InProceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining,pp\. 785–794\.Cited by:[§1](https://arxiv.org/html/2608.13590#S1.p1.1),[§2\.2](https://arxiv.org/html/2608.13590#S2.SS2.SSS0.Px1.p1.2),[§4\.2](https://arxiv.org/html/2608.13590#S4.SS2.p4.7)\.
- T\. Chen, T\. He, M\. Benesty, V\. Khotilovich, Y\. Tang, H\. Cho, K\. Chen, R\. Mitchell, I\. Cano, T\. Zhou, M\. Li, J\. Xie, M\. Lin, Y\. Geng, Y\. Li, J\. Yuan, and X\. contributors \(2025\)Xgboost: extreme gradient boosting\.Note:https://CRAN\.R\-project\.org/package=xgboostR package version 3\.1\.2\.1External Links:[Document](https://dx.doi.org/10.32614/CRAN.package.xgboost)Cited by:[§1](https://arxiv.org/html/2608.13590#S1.p3.1),[§2\.2](https://arxiv.org/html/2608.13590#S2.SS2.SSS0.Px1.p2.2),[§2\.2](https://arxiv.org/html/2608.13590#S2.SS2.SSS0.Px1.p3.5),[§4\.2](https://arxiv.org/html/2608.13590#S4.SS2.p4.7)\.
- J\. H\. Friedman \(2001\)Greedy function approximation: a gradient boosting machine\.Annals of statistics,pp\. 1189–1232\.Cited by:[§1](https://arxiv.org/html/2608.13590#S1.p1.1),[§1](https://arxiv.org/html/2608.13590#S1.p4.1)\.
- T\. Hastie, R\. Tibshirani, and J\. Friedman \(2009\)The elements of statistical learning: data mining, inference, and prediction\.2nd edition,Springer,New York, NY, USA\.Cited by:[§2\.1](https://arxiv.org/html/2608.13590#S2.SS1.p3.13)\.
- P\. J\. Huber \(1992\)Robust estimation of a location parameter\.InBreakthroughs in Statistics,S\. Kotz and N\. L\. Johnson \(Eds\.\),Springer Series in Statistics,pp\. 492–518\.External Links:[Document](https://dx.doi.org/10.1007/978-1-4612-4380-9%5F35)Cited by:[§2\.2](https://arxiv.org/html/2608.13590#S2.SS2.SSS0.Px1.p3.2),[§3\.2](https://arxiv.org/html/2608.13590#S3.SS2.p1.3),[§3](https://arxiv.org/html/2608.13590#S3.p1.4)\.
- X\. Ju and M\. Salibián\-Barrera \(2021\)Robust boosting for regression problems\.Computational Statistics & Data Analysis153,pp\. 107065\.Cited by:[§1](https://arxiv.org/html/2608.13590#S1.p4.1),[§4\.1](https://arxiv.org/html/2608.13590#S4.SS1.p1.14),[§4\.1](https://arxiv.org/html/2608.13590#S4.SS1.p4.4),[§5](https://arxiv.org/html/2608.13590#S5.p3.1)\.
- G\. Ke, Q\. Meng, T\. Finley, T\. Wang, W\. Chen, W\. Ma, Q\. Ye, and T\. Liu \(2017\)Lightgbm: a highly efficient gradient boosting decision tree\.Advances in neural information processing systems30\.Cited by:[§5](https://arxiv.org/html/2608.13590#S5.p4.1)\.
- R\. W\. Lutz, M\. Kalisch, and P\. Bühlmann \(2008\)Robustified l2 boosting\.Computational Statistics & Data Analysis52\(7\),pp\. 3331–3341\.Cited by:[§1](https://arxiv.org/html/2608.13590#S1.p4.1)\.
- R\. A\. Maronna, R\. D\. Martin, and V\. J\. Yohai \(2006\)Robust statistics: theory and methods\.John Wiley & Sons\.External Links:ISBN 978\-0\-470\-01092\-1Cited by:[§2\.2](https://arxiv.org/html/2608.13590#S2.SS2.SSS0.Px1.p4.10),[§3\.1](https://arxiv.org/html/2608.13590#S3.SS1.p1.8),[§3\.2](https://arxiv.org/html/2608.13590#S3.SS2.p1.10),[§3\.2](https://arxiv.org/html/2608.13590#S3.SS2.p4.2),[§3](https://arxiv.org/html/2608.13590#S3.p1.4)\.
- P\. Rousseeuw and V\. Yohai \(1984\)Robust regression by means of s\-estimators\.InRobust and nonlinear time series analysis: Proceedings of a Workshop Organized by the Sonderforschungsbereich 123 “Stochastische Mathematische Modelle”, Heidelberg 1983,pp\. 256–272\.Cited by:[§3](https://arxiv.org/html/2608.13590#S3.p1.4)\.
- J\. W\. Tukeyet al\.\(1977\)Exploratory data analysis\.Vol\.2,Springer\.Cited by:[§2\.2](https://arxiv.org/html/2608.13590#S2.SS2.SSS0.Px1.p4.11)\.
- V\. J\. Yohai and R\. H\. Zamar \(1988\)High breakdown\-point estimates of regression by means of the minimization of an efficient scale\.Journal of the American Statistical Association83\(402\),pp\. 406–413\.External Links:[Document](https://dx.doi.org/10.1080/01621459.1988.10478611),[Link](https://www.tandfonline.com/doi/abs/10.1080/01621459.1988.10478611),https://www\.tandfonline\.com/doi/pdf/10\.1080/01621459\.1988\.10478611Cited by:[§3\.3](https://arxiv.org/html/2608.13590#S3.SS3.p1.12),[§3](https://arxiv.org/html/2608.13590#S3.p1.4)\.
- V\. J\. Yohai \(1987\)High breakdown\-point and high efficiency robust estimates for regression\.The Annals of statistics,pp\. 642–656\.Cited by:[§3\.4](https://arxiv.org/html/2608.13590#S3.SS4.p1.1),[§3](https://arxiv.org/html/2608.13590#S3.p1.4)\.
- Z\. Zhang \(1997\)Parameter estimation techniques: a tutorial with application to conic fitting\.Image and vision Computing15\(1\),pp\. 59–76\.Cited by:[§2\.2](https://arxiv.org/html/2608.13590#S2.SS2.SSS0.Px1.p3.1)\.