Risk-based Calibration for Probabilistic Classifiers
Aritz Pérez Carlos Echegoyen Guzmán Santafé





Risk-based Calibration for Probabilistic Classifiers
    Aritz Pérez, Carlos Echegoyen and Guzmán Santafé
Aritz Pérez is at the Basque Center for Applied Mathematics, 48009 Bilbao, Spain. Email: aperez@bcamath.org
Carlos Echegoyen and Guzmán Santafé are with the Spatial Statistics Group and INAMAT^2, Public University of Navarre, 31006 Pamplona, Spain. Email: {carlos.echegoyen, guzman.santafe}@unavarra.es

    September 9, 2024
====================================================================================================================================================================================================================================================================================================================================================================

§ ABSTRACT 
We introduce a general iterative procedure called risk-based calibration (RC) designed to minimize the empirical risk under the 0-1 loss (empirical error) for probabilistic classifiers. These classifiers are based on modeling probability distributions, including those constructed from the joint distribution (generative) and those based on the class conditional distribution (conditional). RC can be particularized to any probabilistic classifier provided a specific learning algorithm that computes the classifier's parameters in closed form using data statistics. RC reinforces the statistics aligned with the true class while penalizing those associated with other classes, guided by the 0-1 loss. The proposed method has been empirically tested on 30 datasets using naïve Bayes, quadratic discriminant analysis, and logistic regression classifiers. RC improves the empirical error of the original closed-form learning algorithms and, more notably, consistently outperforms the gradient descent approach with the three classifiers.



Supervised classification, probabilistic classifier, empirical risk minimization, iterative learning algorithm, gradient descent.





§ INTRODUCTION
 
Supervised classification is one of the most crucial problems in machine learning, entailing the acquisition of a classifier that minimizes the risk for the 0-1 loss (error). A classifier essentially maps input variables to a set of class labels. However, given the truth probability distribution is unknown, we can not compute the error and classifiers have to be learned from i.i.d. data according to the truth distribution. Traditionally, the learning problem is tackled by minimizing empirical surrogates of the error <cit.>. Support vector machines (SVM) <cit.> and logistic regression (LR) <cit.> exemplify such approaches, minimizing the average hinge and negative log loss on training samples, respectively.

Three main approaches to constructing classifiers are discriminative, conditional, and generative  <cit.>. Discriminative approaches directly model decision boundaries between classes, typified by SVMs. Conditional methods, on the other hand, construct class conditional distributions and determine a classifier by selecting the class label with the maximum probability, as seen in LR. LR achieves this by minimizing the average negative log loss through gradient descent (GD). Generative classifiers present a third alternative, rooted in a joint distribution framework and utilizing Bayes' rule to get class conditional distributions <cit.>. We call probabilistic classifiers conditional and generative approaches.

Unlike their counterparts, generative classifiers do not rely on minimizing an empirical error surrogate but instead on measures quantifying data fitness, such as the maximum likelihood principle. To this category belong classifiers from the exponential family, like Quadratic Discriminant Analysis (QDA), and those grounded in Bayesian networks under varied assumptions, including discrete Bayesian networks <cit.> and conditional Gaussian networks <cit.>. Notably, the direct impact of the likelihood function on classification performance can become negligible as the dimensionality increases <cit.>. Efforts have been made to learn generative classifiers by minimizing the average negative log loss for discrete Bayesian networks through GD <cit.>.

However, gradient descent could suffer several drawbacks, particularly those concerning constrained parameters. Violations of these constraints may yield invalid parameter values, numerical instability, and reduce model performance. Techniques like projection, parameter transformation <cit.>, or Lagrange multipliers <cit.> are indispensable for enforcing constraints but may introduce convergence challenges that hinder efficient optimization. Therefore, specialized methods become imperative to mitigate these challenges and ensure reliable learning of probabilistic classifiers with constrained parameters. Additionally, GD requires differentiable objective functions, typically addressing empirical error through surrogate losses like the negative log loss <cit.>. Furthermore, GD may entail high computational costs; for instance, computing the inverse of a covariance matrix represents a computationally intensive operation that can limit practical applications in high-dimensional domains.

One key advantage of discriminative classifiers concerning generative ones is that they directly model the decision boundary between classes, avoiding the more complex task of estimating the joint distribution, which often results in smaller errors. Additionally, discriminative classifiers tend to require fewer assumptions about the underlying data distribution, making them more flexible and generally more effective when the true data distribution is complex or unknown. 

On the other hand, generative classifiers offer several key advantages. First, they provide comprehensive data modeling by capturing the joint probability distribution of input features and class labels, offering a better understanding of the data-generating process <cit.>. Generative classifiers handle missing data effectively by marginalizing the missing values, enhancing robustness in practical scenarios <cit.>. Generative models can also generate synthetic data samples <cit.>, aiding in data augmentation and anomaly detection. Moreover, they facilitate the integration of prior knowledge and domain expertise through prior distributions in a Bayesian framework. They also fit well into a Bayesian decision theory framework, optimizing decisions under uncertainty <cit.>. Additionally, generative classifiers often perform better with smaller training set sizes, because they approach their best performance faster, possibly with sizes logarithmic in the number of parameters <cit.>. Finally,
When the class conditional distribution is accurately modeled within the joint distribution, generative classifiers can provide optimal predictions for a given cost-sensitive loss function without further adjustments <cit.>. Conditional classifiers represent an intermediate step between discriminative and generative approaches. They directly model the class conditional distribution while avoiding the need to model the marginal distribution of the input features, which is irrelevant for classification <cit.>.

In this work, we present a method that combines the strengths of both probabilistic and discriminative approaches. The proposed method, called risk-based calibration (RC), is designed for learning the parameters of probabilistic classifiers, ensuring their performance is comparable to discriminative classifiers while preserving the advantages derived from modeling probability distributions. focuses on minimizing the empirical risk for 0-1 loss (empirical error) by using learning algorithms that compute the parameters in closed-form from data statistics, such as maximum log-likelihood (ML) or maximum a posterior (MAP) learning procedures.



The rest of the paper is organized as follows. Section <ref> provides the background of the proposal. Section <ref> presents a detailed and formal explanation of the proposed method. Section <ref> introduces the datasets used in the experiments as well as some common aspects of all the experiments. Sections <ref> and <ref> show empirical results for generative and conditional classifiers respectively. Section <ref> summarizes the main conclusions of the current work. Finally, in the appendices, the reader can find further details on the connection of with previous methods (Appendix <ref>), maximum a posteriori estimation of the parameters (Appendix <ref>), the implementation of GD (Appendix <ref>), and additional experimental results (Appendix <ref>).




§ PRELIMINARIES
 



 §.§ Supervised classification

The supervised classification task involves learning a classifier from a training data set that minimizes the expected loss (risk). However, this is often unreliable since we do not know the underlying probability distribution, and the problem is reformulated using surrogate functions of the available training data, such as the average loss in training (empirical risk) <cit.>. 

Let X⊂^n and Y={1, ⋯, r} be the input space and the set of class labels, respectively. A classifier h is a function from instances to labels, h: X→Y, and the set of classifiers is denoted by H. Classifier families can be defined in terms of a particular functional form of their parameters.  We denote by H_Θ={h(;): ∈Θ} the family of classifiers with the parameter space Θ⊂^d. One of the simplest classifiers is the linear discriminant model <cit.>

    h()=max_y_y^T · (1,),

for every input ∈X,  where _y=(θ_0,y,θ_1,y,⋯,θ_n,y) ∈^n+1 and n is the dimension of . In a linear discriminant model, the decision boundaries that separate the classes are linear functions of the inputs  given in terms of hyperplanes in the input space X. A more general form of linear discriminant function is given by

    h()=max_y^T ·ϕ(,y),

where ∈^d are the parameters and ϕ(,y):X,Y↦^d is the feature mapping. Intuitively, the feature mapping defines what is relevant for classifying . The linear discriminant function of Eq. <ref> corresponds to the parameters =(_1, ⋯, _r), with _y ∈^n+1 for y ∈Y, and the one-hot class encoding of the linear function ψ()=(1,), 

    ϕ(,y)= (y=1·ψ(), ⋯, y=r·ψ()),

where · is the indicator function that takes value one when its argument is true, and zero otherwise. We call Eq. <ref> with the linear function ψ()=(1,) linear feature mapping. In general, these classifiers are linear models in the feature space defined by the mapping ϕ(·). Other alternatives to feature mappings that are not linear include polynomial feature mapping <cit.>, radial basis functions <cit.>, Fourier random features <cit.>, and embeddings based on deep learning models.

The loss function measures the discrepancy between the predicted class labels and the true class labels. Formally, the loss of a classifier h evaluated at instance (,y), is a function l: H,(X,Y) → [0,∞). The natural loss in classification is the 0-1 loss, also known as the misclassification loss, l_01(h,(,y))= y≠ h(). Formally, the goal of supervised learning can be defined as selecting the classifier h ∈H that minimizes the risk under 0-1 loss (expected 0-1 loss or error):

    min_h ∈H E_p^* [l_01(h,(,y))],
 
where p^* ∈Δ(X,Y) is the underlying (unknown) distribution of the data. The classifier that minimizes the error is named the Bayes classifier.

In practice, the supervised classification problem is adapted to be tractable. Generally, following a divide-and-conquer approach, the supervised classification problem is addressed for a specific parametric family of classifiers H_Θ, where Θ is the support of the parameters. This approach enables the development of efficient learning algorithms that leverage the functional form of the chosen classifier family. Besides, the adaptations of the supervised classification problem typically involve minimizing a surrogate for the error. In the standard supervised classification settings, the underlying distribution of the data p^* ∈Δ(X,Y) is unknown, and we have access to a supervised training set, (X,Y) ∈X^m ×Y^m, with i.i.d. instances according to p^*,  (X, Y)={(^i,y^i)}_i=1^m. Often, the empirical risk under the 0-1 loss (the average of the 0-1 loss in the training data or empirical error) is used as a surrogate for the true error, and the learning reduces to the empirical error minimization:

    min_h ∈H_Θ1/m∑_,y ∈ X,Y l_01(h,(,y)).

This work is focused on the minimization of the empirical error. For alternatives to the empirical error minimization, see robust risk minimization approaches <cit.>.

Even when the supervised classification problem is restricted to a particular parametric family, optimizing the empirical error can be challenging, and therefore is often replaced by an alternative loss with suitable properties for its minimization. An example of such modifications is learning logistic regression by minimizing the empirical negative log loss, l_log(h,(,y))=- log h(y|), facilitated by its differentiability.



 §.§ Generative classifiers and closed-form learning algorithms


Conditional classifiers are constructed upon a class conditional distribution. The usual approach to learning conditional classifiers is focused on obtaining an accurate model of the class conditional distribution, h(· | ) ∈Δ(Y) for each ∈X, and the classification corresponds to the class label with higher probability h():= max_y h(y|). A typical example of conditional classifiers is LR,

    h(y|)∝^T ·ϕ(,y),

for y ∈Y and for each ∈X.


The generative classifiers are constructed upon a joint probability distribution h(,y) ∈Δ(X,Y), which by the Bayes rule obtains the class conditional distribution h(y|)=h(,y)/∑_y' ∈Yh(,y'). Generative classifiers are fundamentally motivated by their ability to represent the Bayes classifier, assuming they accurately capture the conditional distribution. Generative classifiers are usually constructed upon a joint distribution from a parametric family, such as the exponential family. Common examples of generative classifiers are the quadratic discriminant analysis (QDA) and the classifiers based on Bayesian networks <cit.>. Generative classifiers focus the learning on obtaining a good estimate of the joint distribution, and thus, they use surrogates for the empirical risk indirectly related to classification.
Hence, generative classifiers are typically learned by maximizing the log-likelihood of the joint distribution. In contrast, this work proposes a learning procedure guided by classification performance. 

The current work is mainly devoted to generative classifiers with a closed-form learning algorithm based on statistics obtained from data. We say that a learning algorithm, a, has closed-form when it is a function composition of a statistics mapping function s: X^m,Y^m →^k and a parameter mapping function θ: ^k →Θ, a:= θ∘ s. Statistics mapping summarizes the relevant information in the training data (X,Y) into k statistics, which are used to compute analytically the d parameters of classifier. Generative classifiers based on the exponential family have closed-form algorithms that maximize the likelihood of the training data. In this family, the feature mapping corresponds to statistics mapping, ϕ(·)=s(·), and therefore k=d. Usually, statistics mapping s(·) and feature mapping ϕ(·) are closely related, however, the statistics mapping can include more terms than the feature mapping (see for instance Section <ref>).

From here on, we consider that the statistics mapping involved in the closed-form learning algorithm is additively decomposable, i.e., for (X,Y) ∈X^m,×Y^m, we have that s =s(X,Y)=∑_,y ∈ X,Y s(,y), where with a slight abuse in the notation s(X,Y) and s(,y) denotes the statistics mapping over a training set (X,Y) and over an instance (,y) respectively.




 §.§ Examples of generative classifiers and learning algorithms


Next, we illustrate how the statistics mapping, s(·), manages the statistics calculated from the training sets with two well-known generative classifiers: naïve Bayes (NB) for discrete variables, and quadratic discriminant analysis (QDA) for continuous variables. Most of the classifiers use one-hot encoding for the class-related statistics, and this is the case for both NB and QDA. The statistics mapping for both classifiers can be given by the class one-hot encoding
s(,y)= (y=1·ψ(),⋯ ,y=r·ψ()),
where ψ() extracts the statistics from the features that typically correspond to those required to compute the zeroth, first, and second moments.

In the case of NB with discrete variables, the i-th input feature, x_i, has support X_i={1,...,r_i} for i=1,...,n. NB assumes that the input features are independent given the class variable, which leads to the classification rule for :

    h()= max_y p(y)·∏_i=1^n p(x_i|y),


where p(y) is the marginal probability of the class label y and p(x_i|y) is the probability of i-th input variable taken the value x_i given the class label y. These conditional distributions are assumed to be categorical, and their parameters can be given by the maximum likelihood estimates obtained from the counting statistics. In NB ψ(x)=(1,1=x_1,...,r_1=x_1,...,1=x_n,...,r_n=x_n). For NB the statistics obtained from data are s(X,Y)=(s_0,1,s_1,1,1,...,s_1,r_1,1,...,s_n,r_n,1,...,s_0,r,s_1,1,r,...,s_1,r_1,r,...,
s_n,r_n,r), where s_0,y is the statistics associated to the label y of the class variable and s_i,x_i,y is the statistics associated to the value x_i of input i given the label y. The maximum likelihood parameter mapping is given by p(y)=s_0,y/∑_y' ∈Ys_0,y' and p(x_i|y)=s_i,x_i,y/∑_x_i' ∈X_is_i,x_i',y for y ∈Y, i=1,...n and x_i=1,...,r_i.

On the other hand, the QDA's classification rule is given by: 

    h()= max_y p(y)· |Σ_y|^-1/2· exp{(-μ_y)^T Σ_y^-1 (-μ_y)},


where μ_y and Σ_y are the mean vector and covariance matrix of  given the class y. This classifier is learned by estimating the maximum likelihood mean vector and covariance matrix using the statistics mapping given by ψ()=(1,, ^2), with ^2=·^T being a n × n matrix. The statistics obtained from data are s(X,Y)= (s_0,1,s_1,s^2_1,...,s_0,r,s_r,s^2_r) with s_y'=∑_,y ∈ X,Yy'=y·x, and s^2_y'=∑_,y ∈ X,Yy'=y·^2 for y' ∈Y. The components 1, , and ^2 are used to get the maximum likelihood estimation for p(y), μ_y, and Σ_y using the parameter mapping p(y)=s_0,y/∑_y' ∈Ys_0,y', μ_y= s_y/s_0,y and Σ_y= s^2_y/s_0,y - μ_y ·μ_y^T, for each class label y ∈Y.

Using the statistics mapping described here, we can also adapt the parameter mapping for NB and QDA to create a closed-form algorithm that maximizes the a posteriori distribution using conjugate priors (see Appendix <ref> for further details).



§ RISK-BASED CALIBRATION
 

The risk-based calibration algorithm (RC) is an iterative heuristic method to improve the empirical risk under 0-1 loss of a probabilistic classifier using a closed-form learning algorithm. This method is founded on the very basic intuition of modifying the statistics used by the closed-form learning algorithm, guided by the stochastic 0-1 loss for each instance. The stochastic 0-1 loss of the probabilistic classifier h in (,y) is given by 

    l_s01 (h,(,y)) = 1-h(y|)= ∑_y'≠ y h(y'|),

that corresponds to the expected loss of a randomized classifier that selects
label y with probability h(y|) for y ∈Y.

We aim at finding the parameters ^* that minimize the empirical risk of the stochastic 0-1 loss (empirical stochastic error). Given a training set (X,Y), the empirical stochastic error is zero when h(y|)=1 for all ,y ∈ X,Y. Intuitively, it is possible to lead the statistic s, and thus the model parameters, θ(s), towards the optimal classifier by strengthening the statistics s(,y) and weakening s(,y') for ,y ∈ X,Y and every y' ∈Y with y' ≠ y. The strengthening-weakening update is given by the classifier's stochastic 0-1 loss at point (,y).  

We propose to shrink the empirical stochastic error by raising h(y|) and by dropping h(y'|) for (,y) according to the stochastic 0-1 loss by calibrating the set of statistics, s, used to obtain the model parameter. The calibration of s is performed by adding s(,y) with a weight 1 - h(y|) (strengthening h(y|)), and by subtracting s(,y') with a weight h(y'|) for all y'≠ y (weakening h(y'|)). The strengthening-weakening calibration is directly give by Ep. <ref> Combining both terms, given the classifier h, we have the following updating rule of the statistics s given (,y):

    s= s + s(,y) - ∑_y' ∈Y h(y'|) · s(,y').

Given the data (X,Y) and the classifier h, and due to the additive nature of the statistics s, the updating rule is simply given by:

    s= s + s(X,Y) - s(X,h),
 
where s(X,h)= ∑_∈ X∑_y' ∈Y h(y'|)· s(,y') is the probabilistic estimate of the statistics given the probabilistic classifier h. 

The procedure is described in Algorithm <ref>, where, at each iteration, h corresponds to the probabilistic classifier with parameters θ(s), and lr > 0 is the learning rate.

(X,Y)


  * s← s(X,Y)

  * 


  * s←s - lr ·(s(X,h)- s(X,Y)), for h with θ(s)


  * Stop criterion is met



The initialization of the statistics (step 1 in Algorithm <ref>) can be arbitrary, as long as the statistics remain consistent and produce valid parameters. However, we recommend initializing them using s(X,Y). The statistics obtained from training data provide a more competitive starting point from the empirical error point of view. However, it is possible to try different runs using bootstrap samples from the training data to avoid the convergence to a poor local optima.

The computational cost of each iteration is given by the statistics mapping s(·) and parameter mapping θ(·). The computational complexity of 
the statistics mapping is linear in the number of training samples 
m and the dimension of the statistics k, O(d · k), while the computation of the parameter mapping is independent of m
and is usually linear in the number of statistics and parameters O(d + ||). Due to the additively decomposable assumption for the statistics mapping, it is possible to speed up the by using stochastic and minibatch versions, which process subsets of the data in each iteration.

A relevant property of is that the sample size of the updated statistics remains invariant to iterations. This is because the strengthening-weakening strategy comes from the stochastic 0-1 loss (Eq. <ref>) and satisfies ∑_y' ∈Yy=y'-h(y'|)= 0 for any ,y ∈X,Y. However, discrepancies between the class conditional probabilities given by the classifier and the class labels in the training set can sometimes lead to invalid parameters. To address this, we simply allow the learning process to continue by updating only those statistics that would result in valid parameters, freezing the rest.

The is a general-purpose algorithm for learning probabilistic classifiers that presents some connections to other existing methods proposed in the literature. The analysis of the connections with three iterative learning algorithms is presented in Appendix <ref>: Discriminative frequency estimate for classifiers based on discrete Bayesian networks <cit.>, GD for LR, and the TM algorithm for generative classifiers from the exponential family <cit.>.



§ DATASETS AND EXPERIMENTAL SETTING
 

In this section, we will introduce the datasets used in the experiments and some common elements to all the experiments in the current paper. All the experiments are focused on the minimization of the empirical error. 

The implementations of the classifiers, learning algorithms, and experiments are available online in the public Python repository at <https://gitlab.bcamath.org/aperez/risk-based_calibration>.

We utilize 30 public available datasets <cit.>, each characterized by different numbers of instances (m) and variables (n).
Table <ref> provides an overview of the datasets used for the experiments, where the "Index" column serves as the identifier referenced along the empirical analysis. This comprehensive collection covers a wide range of domains and complexities, ensuring a robust evaluation of the proposed method across diverse real-world scenarios. Note that the 512 features of the datasets {3, 4, 8, 20, 30} are provided by ResNet18 <cit.> pre-trained deep neural networks for image classification problems. These specific datasets contain a larger number of instances and input variables compared to the others.

All algorithms have been run 64 iterations with a fixed learning rate of lr= 0.1. The performance is measured in terms of the empirical error. We use two closed-form learning algorithms for the probabilistic classifiers: Maximum Likelihood (ML) and maximum a posteriori (MAP) (see Section <ref> and <ref>, and Appendix <ref>, respectively). The parameters of the classifiers are also initialized according to ML and MAP, respectively.

The proposed method, , is compared with GD using the same parameter initialization. The main results are summarized in tables indexed by dataset with the following description of the columns: The "ML" (or "MAP") column shows empirical error given by ML (or MAP) parameters and constitutes the initialization of and GD; the "" and "GD" columns provide the minimum empirical error reached by and GD algorithm in 64 iterations, respectively, and the best result for each dataset is highlighted in bold; the "Iter" columns contain the iteration at which and GD reach the minimum empirical error in 64 iterations. The "Reach" column shows the number of iterations required by RC to achieve an error that is less than or equal to the lowest error obtained by GD in 64 iterations. As the reach is smaller, it indicates that the error reduction of RC is steeper and faster than that of GD. The "avg." row presents the average errors and number of iterations across the 30 datasets, and the average reach is computed with those datasets in which obtains an error that is less than or equal to the lowest error of GD.



































§ EXPERIMENTS WITH NB AND QDA
 

NB and QDA were introduced in section <ref> as illustrative examples of generative classifiers alongside the closed-form algorithm that utilizes to perform the calibration of these classifiers. In this section, we summarize the empirical results obtained with NB and QDA. 



 §.§ for NB


NB is a classification model that deals with discrete variables. To adapt the datasets from Table <ref> for use with NB models, each continuous variable is discretized into 5 categories according to a k-means strategy <cit.>, where the values in each bin belong to the same cluster. The results are presented in Tables <ref> and <ref>. 


According to Table <ref>, 
improves the ML initialization in 27 out of 30 datasets. In all these cases, obtains lower errors than GD. The average reach of is 5 indicating that it produces a steeper and faster empirical error reduction than GD.



In Table  <ref>, when MAP learning is used, reduces the initial error across 21 datasets, with proving superior to GD in 18 of them. The average reach is 5. The behavior of is clearly superior to GD for learning NB.










 §.§ for QDA


The experimental results obtained with QDA are summarized in Tables <ref> and <ref> for ML and MAP, respectively. 
According to Table <ref>, 
improves the ML initialization in 29 out of 30 datasets. In these datasets, obtains lower errors than GD except in (5) where they tie. The average reach in this scenario is 6. In Table  <ref>, when MAP learning is used, also reduces the initial error in 29 datasets. In all these datasets, achieves lower errors than GD except in three cases, obtaining equal results in two of them. The average reach is also 6. Again, the behavior of is clearly superior to GD for learning QDA.







 









§ FOR LOGISTIC REGRESSION


One of the most popular conditional classifiers is LR (see Eq. <ref>), which is typically learned using GD to minimize the empirical risk under the negative log loss. In this section, we present a closed-form algorithm for LR, enabling its learning using . The proposed closed-form learning algorithm is based on a generative formulation of LR under parametric assumptions.





 §.§ Closed-form algorithm for LR
 

Let h(,y)= h(y)∏_i=1^n h(x_i|y), where h(y) is a categorical distribution with parameters p_1,...,p_r, h(x_i|y) is a Gaussian distribution with mean μ_i,y and variance σ_i^2 for i=1,...,n and y ∈Y. This generative classifier corresponds to the naïve Bayes classifier based on conditional Gaussian networks <cit.> under homoscedasticity assumption, i.e., the variance of each continuous variable does not depend on the class. The connection between both models is explained and analyzed in <cit.>. By the Bayes rule h(x,y) leads to the conditional class distribution:

    h(y|) ∝ exp{ϕ(,y)^T· - c_y},

where ϕ(,y) is the linear feature mapping, and =(_1,...,_r) are the parameters with _y=(θ_0,y, θ_1,y,⋯θ_n,y) for y ∈Y; and θ_0,y= log p_y, θ_i,y= logμ_i,y/σ_i^2, and c_y= ∑_i=1^n μ_i,y^2/σ_i^2 + 1/2 logσ_i^2 for y=1,...,r and i=1,...,n. By grouping terms, the conditional class distribution corresponds to Eq. <ref> with =(_1,...,_r), θ_y,0= c_y + log p_y for y=1,...,r.


Now, we are in a position to define a closed-form learning algorithm for LR that maximizes the likelihood of the data (ML). The statistics mapping is simply given by the concatenation of the linear feature mapping, and the squares of the input features (x_1^2,...,x_n^2), s(,y)= (ϕ(,y),x_1^2,...,x_n^2), with s(X,Y)= (s_0,1,s_1, ⋯, s_0,r,s_r, s^2) with s^2=(s_1^2,...,s_n^2)=∑_∈ X(x_1^2,...,x_n^2). The parameter mapping corresponds to p(y)= s_0,y/∑_y' s_0,y', the mean vector μ_y= (μ_1,y,...,μ_n,y)= s_y/s_0,y for y ∈Y, and the variance vector σ^2= (σ_1^2,...,σ_n^2)= s^2/∑_y'∈Ys_0,y' - ∑_y' ∈Yp(y)·μ_y^2 with μ_y^2=(μ_1,y^2,...,μ_n,y^2).[In this work, we are considering ψ(·)=(1,) to be linear on , but the proposal can be easily extended to arbitrary ψ(·) by using its corresponding class one-hot feature mapping (Eq. <ref>), and a statistics mapping corresponding to ϕ(,y) concatenated with the squares of ψ(·), s(,y)=(ϕ(,y), ψ()^2).]
Using the same statistics mapping, we can adapt the parameter mapping to create a closed-form algorithm that maximizes a posteriori distribution using the conjugate priors (see Appendix <ref> for further details).



 §.§ Experiments using for LR


Next, we provide the set of experiments on the minimization of the empirical error for LR. The experiments have been performed with using the closed-form algorithms given in <ref> (ML and MAP) and GD of the empirical risk under the negative log loss. and GD starts from the same initialization corresponding to ML and MAP parameters obtained from data. Additionally, in Appendix <ref> we include results for random initialization that highlights the robustness and efficiency of .
 
Tables <ref> and <ref> show the performance of for ML and MAP, respectively. reduces the initial error in 29 datasets for both initializations, obtaining better results than GD in 24 out of 29 cases with ML and in 25 out of 29 cases with MAP. The average reach of is 13 and 16, for ML and MAP respectively. 
Once again, the results show the ability of to improve the closed-form algorithms and its superiority concerning GD for learning LR. 









§ CONCLUSIONS
 

This work proposes an iterative learning algorithm called risk-based calibration () to minimize the empirical error of probabilistic classifiers. 

can be used to learn any probabilistic classifier, whether generative or conditional, as long as the classifier has a closed-form learning algorithm that involves two steps: statistics collection from data and analytical parameter computation from the obtained statistics. Typical examples include classifiers from the exponential family, such as naïve Bayes and quadratic discriminant analysis, which can benefit from the algorithm. Additionally, we show how to use the proposed procedure to learn the logistic regression classifier. The main difference between and other alternatives is that focuses on the calibration of the statistics rather than the parameters. This calibration is performed using information from the stochastic 0-1 loss function, ensuring that the learning process directly targets the minimization of the empirical error.

In the experiments, consistently achieved lower empirical errors than the gradient descent approach when learning naïve Bayes, quadratic discriminant analysis, and logistic regression classifiers. This demonstrates the effectiveness and robustness of to learn probabilistic classifiers by minimizing the empirical error. The ability of to reduce the empirical error reveals its potential as a preferred choice when learning probabilistic classifiers from data. 

The implementations of the classifiers, learning algorithms, and experiments are available online in the public Python repository at <https://gitlab.bcamath.org/aperez/risk-based_calibration>.

IEEEtran





§ CONNECTIONS BETWEEN AND OTHER METHODS
 

This section shows the connection of with three iterative learning procedures: Discriminative frequency estimate for learning classifiers based on discrete Bayesian networks, GD for LR, and the TM algorithm for generative classifiers from the exponential family.



 §.§ Discriminative frequency estimate


Discriminative Frequency Estimate (DFE) <cit.> is used for learning the parameters of classifiers based on discrete Bayesian networks <cit.>, such as naïve Bayes (NB). The main motivation for DFE is that learning discrete Bayesian networks based on gradient descent over the empirical risk under the negative log loss <cit.> is computationally demanding. DFE is an iterative procedure where, at each iteration, the statistics used to derive the parameters of a discrete Bayesian network classifier are updated. Using our notation, the updating rule is given by:

    s= s + ∑_,y ∈ X,Y(1-h(y|))· s(,y)

This corresponds to a heuristic that strengthens the statistics s(,y) according to the stochastic 0-1 loss. The main difference from in the context of discrete Bayesian network classifiers is that DFE does not weaken the incorrect statistics associated with the instance (,y), s(,y') for y' ≠ y. An important negative consequence is that DFE increases the equivalent sample size of the statistics 
s, at each iteration by 
O(ϵ_s01· m), where ϵ_s01 is the empirical risk under the stochastic 0-1 loss of the classifier obtained in the previous iteration. This can have dramatic effects with large training sets or when the number of iterations for convergence is large. The main limitation of DFE is that it only applies to discrete Bayesian network classifiers.

In Table <ref> of Appendix <ref> we have summarized an experimental comparison between using ML and DFE for learning NB. These results clearly show that outperforms DFE.




 §.§ The connection between and gradient descent for LR


In the case of LR, is equivalent to GD for the empirical risk under the negative log loss when uses a particular closed-form learning algorithm. This algorithm is related to the empirical average of the feature mapping used by LR. In Section <ref>, we present LR with linear feature mapping, while here we consider Eq. <ref> an arbitrary function ψ(·).


Let h be a classifier given by Eq. <ref> for the feature mapping given in Eq. <ref> with arbitrary ψ(·), a= θ∘ s be a learning algorithm for h, and (X,Y) be a training set of size m. The update rule is equivalent to GD of the empirical risk under the log loss if the learning algorithm a is given by the statistics mapping s(,y)= ϕ(,y) and the parameter mapping θ(s)=s/m.


On the one hand, the updating rule with learning rate lr>0 is given by:

    s_y'= s_y' - lr·∑_,yϕ(x)· (h(y'|)- y'=y),

for s_y'∈Y. Then, using the parameter mapping, we have that

    =  - lr/m·∑_,y ∈ X,Y∑_y' ∈Yϕ(,y) · (h(y'|) - y'=y)


On the other hand, the gradient descent of the average negative log loss of h, R(h)=-1/m ∑_(,y)log h(y|), with respect to  is

    δ R(h)/δ=1/m∑_,y ∈ X,Y∑_y'∈Yϕ(,y)· (h(y'|)-y'=y),

which leads to the same parameter updating rule as .






 §.§ The connection between and TM


The TM <cit.> is an iterative algorithm for maximizing the conditional likelihood using maximum likelihood learning procedures. Conditional log likelihood is proportional to the empirical risk under the negative log loss. TM is a general purpose algorithm that can be used for both regression (Y⊂ continuous) and classification (Y={1,⋯ , r} categorical).

At each iteration t, the TM solves two steps:

    
  * T-step: Compute the gradient of the marginal log-likelihood with respect to the parameters: 
    δ LL(X;)/δ,

    with LL(X,)= ∑_∈ Xlog∑_y ∈Y h(,y;).
    
  * M-step: solve the maximization problem.
    
    ^(t+1)= max_ LL(X,Y;) - δ LL(X;^(t))/δ·,
     

where δ LL(X;^(t))/δ denotes the derivative of the marginal log-likelihood with respect to the parameters  evaluated at ^(t).
The T-step is based on an approximation to the conditional log-likelihood function, obtained by linearizing the marginal log-likelihood. 

In classification, Y={1,⋯ , r}, the M-step is equivalent to finding  so that the next equality holds:

    δ LL(X,Y;)/δ = δ LL(X,h(· | X;^(t)))/δ,

with LL(X,h(· | X;))=∑_∈ X∑_y ∈Y h(y|;) ·log h(,y; ). The TM differentiates the minimal sufficient statistics that depend on the class variable, u, from those independent of the class v, s(X,Y)=(u(X,Y),v(X)). For the particular case of generative classifiers from the exponential family, the TM iterates by updating the minimal sufficient statistics that depend on the class. Following our notation, it can be shown that TM reduces to:

    u^(t+1)= u^(t) + u(X,Y) - u(X,h^(t)),

where the parameters of the model at iteration t+1 are given by the maximum likelihood sufficient statistics s=(u^(t+1),v^(0)), with v^(0) = v(X). In summary, for the particular case of generative classifiers from the exponential family using the maximum likelihood learning algorithm, and a learning rate of lr=1, TM and are equivalent. 


We believe the main reason for the limited use of TM is the difficulty of its implementation. This is mainly due to the necessity of understanding the intricate details of the exponential family model, such as minimal sufficient statistics and those dependent or independent of the class variable. For instance, one of the few examples of TM usage is outlined in <cit.>, where the authors used it to learn classifiers based on Bayesian networks with categorical variables. The primary methodological challenge in that work involves updating the minimal sufficient statistics of conditional distributions over categorical variables while maintaining their consistency. This unnecessarily complicates the implementation of the learning method, particularly since its complexity strongly depends on the number of states of the categorical variables. This complexity is apparent when compared to , which merely requires the iterative application of maximum likelihood learning with probabilistically labeled data.



§ MAXIMUM A POSTERIORI PARAMETER MAPPING


An alternative to the maximum likelihood parameter estimation is a Bayesian estimation of the parameters. In Bayesian estimation, we assume a prior for the distribution of the parameters. Then, given the data, we obtain the posterior distribution of the parameters and select the parameters of its mode (maximum a posteriori parameters, MAP). For certain, parametric distributions, there are prior distributions over their parameters that allow obtaining the MAP parameters in closed form.

NB is based on the categorical distribution and the Bayesian conjugate of the categorical distribution is the Dirichlet distribution. Let's take the categorical distribution p(y) for y ∈Y={1,...,r} and the prior Dirichlet distribution for its parameters

    =(θ_1,...,θ_r) ∼ Dir(α),

with hyperparameters α=(α_1,...,α_r). The posterior distribution of the parameters of p(y) after observing Y={y_i}_i=1^m is given by:

    ∼ Dir(α + (m_1,...,m_r)),

where m_y'= ∑_y ∈ Yy'=y for y' ∈Y. The MAP corresponds to:

    '_y= m_y + α_y - 1/m + α - r,

with α= ∑_y=1^r α_y. By taking α_y= m_0/r +1 for y ∈ 1 we have the more intuitive MAP

    '_y= m_y + m_0/r/m + m_0,

where m_0 can be interpreted as the equivalent sample size of the prior. In the experiments with MAP we have taken m_0= r. The same analysis follows for all the conditional distributions that are involved in NB, p(x_i|y) for y ∈Y, and i=1,...,n with x_i ∈{1,...,r_i}.

QDA is based on a categorical distribution p(y) and n-dimensional Gaussian density functions p(|y) for y ∈Y. The Bayesian conjugate of the parameters of a multivariate Gaussian distribution is the normal distribution for the mean and the inverse-Wishart distribution for the covariance:

    μ|Σ∼   N(μ_0,1/κ_0Σ)
    Σ∼   W^-1(T_0, ν_0)

The posterior distribution for the parameters given the observations X={_1,...,_m} is given by

    μ|Σ∼N(μ', 1/κ'Σ)
    Σ∼W^-1(Σ',ν')

Then, the MAP parameters for the mean and the covariance matrix are

    μ'   =   κ_0μ_0 + mμ̂/κ_0 + m
    Σ'   =   T_0 + mΣ̂ + κ_0 · m/κ_0 + m(μ̂-μ_0)· (μ̂-μ_0)^T/ν_0+m+n+1,

where μ̂ and Σ̂ are the sample mean and covariance matrix. We propose the following re-parametrization in terms of m_1, m_2 ≥ 0, κ_0= m_1, T_0= m_2Σ_0 and ν_0=(m_2- n -1), for the sake of interpretability. Under this re-parametrization and neglecting the term (κ_0 · m)(μ̂-μ_0)· (μ̂-μ_0)^T/((ν_0+m+n+1)·(κ_0 + m)) because usually κ_0<<m, we have the next intuitive MAP parameters in terms of prior mean vector μ_0 and covariance matrix Σ_0 with weights m_1 and m_2, respectively:

    μ'= m_1·μ_0 + mμ̂/m_1 + m
    Σ'= m_2·Σ_0 + mΣ̂/m_2 + m

Here, m_1 and m_2 can be interpreted as the equivalent sample size of the priors for the mean vector and covariance matrix, respectively. In the experiments with MAP we have taken m_1=m_2=10.
This MAP estimate is used for the parameters of p(|y) for y ∈Y, while the MAP parameters of p(y) are obtained using the same procedure of NB. We have followed a similar approach for the parameter mapping used in the closed-form algorithm of LR.



§ GRADIENT DESCENT


In this section, we give the details of the gradient descent updating rule for the empirical risk under the negative log loss for the three classifiers considered in the experiments of this paper. In addition, we explain the transformations used to satisfy the parameter constraints. All the classifiers have been expressed into the exponential family form, 
    h(y|) ∝η^T· s(,y) + A_y(η^T).
 This simplifies obtaining a clear expression for the gradient of the parameters and puts all the considered models into a comparable form.

The exponential family form of the conditional probability modeled by NB is given by feature mapping corresponding to the one-hot encoding of each feature and the class
ϕ(,y)=(1=y,1=y1=x_1,...,1=yr_1=x_1,...,,1=y,1=y1=x_n,...,r=yr_n=x_n,...,r=y1=x_n,...,r=yr_n=x_n) and the parameters
η=(η_0,η_1,0,...,η_n,0,... ,η_1,r,...,η_n,r) with η_0=(η_0,y= log p(y))_y=1^r and η_i,y=(η_i,x_i,y = log p(x_i|y))_x_i=1^r_i for i=1,...,n and y ∈Y. For NB A_y(η^T)=0. The gradient descent updating rule is given by:

    η_0,y'=   η_0,y' - lr/m∑_,y ∈ X,Y  (h(y'|)- y'=y)
    η_i,x_i',y'=   η_i,x_i',y'
       - lr/m∑_,y ∈ X,Yx_i'= x_i·(h(y'|) - y'=y)

for y'=1,...,r, i=1,...,n, and x'_i=1,...,r_i.

Then, after applying the gradient descent updating rule, the natural parameters are transformed into probabilities by exponentiation and by projecting them into the simplex. Alternatives to the projection to the simplex include using softmax to obtain proper probability distributions. Unfortunately, by using these transformations the descent in the average negative log loss is no longer guaranteed.

The exponential-family distribution form of the conditional probability modeled by QDA is given by the feature mapping corresponding to the class one-hot coding ϕ(,y)=(1=y,1=y·,1=y·^T,..., r=y,r=y·,r=y··^T); the parameters η=(η_0,η_1,η_2) with η_0= (η_0,y=log p(y))_y=1^r, η_1= (η_1,y= Σ_y^-1·μ_y)_y=1^r and η_2=(η_2,y= -1/2 ·Σ_y^-1)_y=1^r, being μ_y and Σ_y the mean vector and covariance matrix conditioned to y; and A_y(η)= 1/4·η_2,y^-1·η_1,y·η_1,y^T ·η_2,y^-1 - 1/2·η_2,y^-1. The gradient descent updating rules are given by:

    η_1,y' =   η_1,y' - lr/m·∑_x,y (h(y'|) - y=y')
       ·(x + 1/2 ·η_2,y'^-1·η_1,y')
    η_2,y' =   η_2,y' - lr/m·∑_,y ∈ X,Y (h(y'|) - y=y')
       ·(^2 - 1/4 ·η_2,y'^-1·η_1,y'·η_1,y'^T ·η_2,y'^-1 + 1/2· tr(η_2,y'^-1))

for y'=1,...,r, and being tr(·) the trace of a matrix. Then, after every gradient descent updating the natural parameters η_0 are transformed into probabilities by exponentiation and by projecting them into the simplex; η_2,y is transformed into the covariance matrix Σ_y=-1/2·η_2,y^-1 and is ensured to be a positive semi-definite matrix by: i) obtaining the singular value decomposition, ii) guaranteeing that all the eigenvalues are no smaller than ϵ= 10^-2, and iii) reconstructing the matrix using the constrained eigenvalues. Again, by transforming the obtained parameters to fulfill their associated constraints (probabilities and covariance matrices) we can not ensure that the average negative log loss descents.

LR is directly given in the exponential-family distribution form, η=. In this work, we simply consider the feature mapping corresponding to the one-hot class encoding ϕ(,y)=(1=y, 1=y·,..., r=y, r=y·) and the parameters =(_0,_1) with _0=(θ_0,y∈)_y=1^r _1=(_1,y∈^n)_y=1^r. In this model A_y()=0. The gradient descent updating rule for LR is:

    =  + lr/m∑_,y ∈ X,Yϕ(,y).

In this model, the parameters have no constraint, and thus it is guaranteed a monotonic descent of the average negative log loss, but not in terms of the empirical error.



§ ADDITIONAL RESULTS




 §.§ Convergence curves in mnist


Next, we show the curves of the evolution of the empirical error with respect to the number of iterations of and GD in mnist dataset (20) for NB, QDA, and LR. The error curves represent 128 iterations.



In Figure <ref>, we show the error curves of using ML and MAP for learning NB. can reduce the error of the initial model using ML and MAP parameters. with MAP converges to a local optimum in 60 iterations, while with ML does not reach convergence in 128 iterations. However, the behavior using ML is better than MAP for mnist.





Figure <ref> shows the error curves of QDA with ML and MAP. shows the same behavior with ML and MAP. RC shows a significant reduction in error in the first iterations and achieves an error close to the minimum in less than 30 iterations.




Figure <ref> shows the error curves of LR with ML and MAP. with ML and MAP shows a similar behavior. In 128 iterations, reduces the error from 0.14 to 0.05, and the curve's downward slope suggests that the error will keep decreasing if more iterations are used. The error curve of decreases faster than that obtained with GD.








 §.§ Random initialization of parameters

Table <ref> shows the behavior of for LR with random initialization and updating the parameters using ML. In this case, the column "Random" contains the empirical error when the model is initialized at random. obtain lower errors than GD in 27 datasets. In two of the three remaining datasets, both algorithms tie. The results from Table <ref> suggest that reaches the best GD result, on average, in 19 iterations.
The minimum errors and number of iterations required to obtain them with using ML and random initialization are similar (see Table <ref>). Figure <ref> and Figure <ref> show the evolution of the error with random and ML initialization in mnist, respectively. When departing from random initialization produces a steeper descent in the first iterations and obtains a similar error to the ML initialization in only 2 iterations. These results suggest that exhibits a robust behavior independently of the initial parameters.
  













 §.§ Comparison with DFE





Table <ref> summarizes the comparison between with ML as closed-form learning algorithm and the discriminative frequency estimate (DFE) for discrete NB (see Appendix <ref>). achieves lower errors in 25 datasets and equal in 5. The average reach of is 12, indicating that it produces a steeper and faster empirical error reduction than DFE. These results clearly show that using ML is better than DFE for learning NB.



