The part of the code that you are pointing out is the computation of the intercept. The estimator of a model minimizing the mean absolute error is the median. So the computation of the intercept looks logical to me.

From my understanding, optimizing Mean Absolute Error and Median Absolute Error should give different results, but I believe that either I have not been able to communicate my question clearly, or have I not provided sufficient backing for my question, or I have confused the code for the intercept above.

Any resources would be much appreciated.

I would appreciate if there was more explanation regarding on using custom loss functions Custom Loss function? #21614. This way, we could optimize for Mean Absolute Error, Median Absolute Error, Log Cosh, etc. as per the requirement.

Aside from that, I will try my best to research further and formulate a better question in the future. You may close this issue. Thank you!

The median is the minimizer of the mean absolute error. Let's sample same data distributed in such a way that the mean and the median are different:

import numpy as np

data = np.random.lognormal(mean=0, sigma=1, size=1000)
np.median(data).round(4), np.mean(data).round(4)

(np.float64(0.9989), np.float64(1.6216))

Let's find the minimizer of MAE using scipy.optimize:

from scipy.optimize import fmin


def mae(x, data):
    return np.mean(np.abs(data - x))

x = np.zeros(shape=1, dtype=np.float64)
fmin(mae, x, args=(data,), disp=False)

array([0.99875])

Using HistGradientBoostingRegressor(loss="absolute_error") optimizes for the median of absolute errors

No, HistGradientBoostingRegressor(loss="absolute_error") minimizes the sum or mean absolute error of the conditional estimates (as does the quantile loss with q=0.5). Both choices lead to an estimate of the conditional median of Y|X as a result.

Please point us to the specific part of the documentation that was not clear and suggest a way to improve it. Otherwise I think we can close.

The median is then minimizer of the mean absolute error.

Noted. You are correct; I replicated your code, thanks!

No, HistGradientBoostingRegressor(loss="absolute_error") minimizes the (conditional) mean absolute error (as does the quantile loss with q=0.5) and both estimate the conditional median of Y|X as a result.

But code for quantile = 0.5 (ie, median) gives me a different result from that of mean; this is my point:

def quantileae(x, data):
    return np.quantile(np.abs(data - x), 0.5)

x = np.zeros(shape=1, dtype=np.float64)
fmin(quantileae, x, args=(data,), disp=False)

array([0.58675])

def medianae(x, data):
    return np.median(np.abs(data - x))

x = np.zeros(shape=1, dtype=np.float64)
fmin(medianae, x, args=(data,), disp=False)

array([0.58675])

Please point us to the specific part of the documentation that was not clear and suggest a way to improve it. Otherwise I think we can close.

Problem
The term "squared_error" and "absolute_error" is not very clear: ie, is it Mean Squared Error or Median Squared Error; Mean Absolute Error or Median Absolute Error. The documentation should make it clear that it is Mean Squared Error and Mean Absolute Error.

Specific part of the documentation
https://scikit-learn.org/1.5/modules/generated/sklearn.ensemble.HistGradientBoostingRegressor.html

loss{‘squared_error’, ‘absolute_error’, ‘gamma’, ‘poisson’, ‘quantile’}, default=’squared_error’
The loss function to use in the boosting process. Note that the “squared error”, “gamma” and “poisson” losses actually implement “half least squares loss”, “half gamma deviance” and “half poisson deviance” to simplify the computation of the gradient. Furthermore, “gamma” and “poisson” losses internally use a log-link, “gamma” requires y > 0 and “poisson” requires y >= 0. “quantile” uses the pinball loss.

Suggestion

The loss function to use in the boosting process. Note that the “squared error”, “gamma” and “poisson” losses actually implement “half least squares loss”, “half gamma deviance” and “half poisson deviance” to simplify the computation of the gradient. Furthermore, “gamma” and “poisson” losses internally use a log-link, “gamma” requires y > 0 and “poisson” requires y >= 0. “quantile” uses the pinball loss. “absolute_error” refers to Mean Absolute Error; for other variants of absolute error, you need to write your own loss function, “squared_error” refers to Mean Squared Error; for other variants of squared error, you need to write your own loss function

But code for quantile = 0.5 (ie, median) gives me a different result from that of mean; this is my point

The quantile loss of HistGradientBoostingRegressor is not equivalent to the code you wrote above. The quantile loss is also known as the pinball loss, and it is documented here:

https://scikit-learn.org/stable/modules/model_evaluation.html#pinball-loss

Here is the equivalent Python code:

def quantile_loss(x, data, q=0.5):
    return np.mean(np.maximum(q * (data - x), (q - 1) * (data - x)))

x = np.zeros(shape=1, dtype=np.float64)
fmin(quantile_loss, x, args=(data,), disp=False)

array([0.99875])

What we could do is expand the following section:

https://scikit-learn.org/stable/modules/ensemble.html#loss-functions

to link to the matching metrics function from https://scikit-learn.org/stable/modules/model_evaluation.html.

To make that more explicit and reference this section from the docstring.

“absolute_error” refers to Mean Absolute Error; for other variants of absolute error, you need to write your own loss function, “squared_error” refers to Mean Squared Error; for other variants of squared error, you need to write your own loss function

I would not write that: for all loss functions, it's always computed as the sum (or equivalently the mean) of the individual loss function computed on each data point in the training set: this is the case for absolute and squared error but also for quantile, Poisson, gamma... Furthermore, we should not advertise for custom loss functions in the doc as there is no official public API to do so in scikit-learn.

we should not advertise for custom loss functions in the doc as there is no official public API to do so in scikit-learn.

Noted. That is right.

for all loss functions, it's always computed as the sum (or equivalently the mean) of the individual loss function computed on each data point in the training set

This is the point I'm not able to understand. Why is that so: I'm assuming this is implied or an unspoken rule in Machine Learning that mean is the standard aggregate function for loss functions, is it?

So far, I've always specified mean/median whenever I mention a loss function as I was not aware of such. If you could share any resource that I could read, I would appreciate it.

This is the point I'm not able to understand. Why is that so: I'm assuming this is implied or an unspoken rule in Machine Learning that mean is the standard aggregate function for loss functions, is it?

You probably want to refer to:
https://en.wikipedia.org/wiki/Empirical_risk_minimization

You are interested in minimizing the risk that is defined as the expectation of the loss function and thus the sum aggregate.