Implement temperature scaling for (multi-class) calibration

### Describe the workflow you want to enable

It would be great to have temperature scaling available as a post-hoc calibration method for binary and multi-class classifiers, for example in `CalibratedClassifierCV`.

### Describe your proposed solution

Temperature scaling is a simple, efficient, and very popular post-hoc calibration method that also naturally supports the multi-class classification setting. It has been proposed in Guo et al. (2017) with >5000 citations, so it meets the inclusion criterion: http://proceedings.mlr.press/v70/guo17a.html
It also does not affect rank-based metrics (if the temperature is restricted to positive values) unlike isotonic regression (https://github.com/scikit-learn/scikit-learn/issues/16321). Moreover, it avoids the infinite-log-loss problems of isotonic regression.
Temperature scaling has been discussed in https://github.com/scikit-learn/scikit-learn/discussions/21785
I experimented with different post-hoc calibration methods on 71 medium-sized (2K-50K samples) tabular classification data sets. For NNs and XGBoost, temperature scaling is competitive with isotonic regression and considerably better than Platt scaling (if Platt scaling is applied to probabilities, as implemented in scikit-learn, and not logits). For AUC, it is considerably better than isotonic regression. 

Here is a simple implementation using PyTorch (can be adapted to numpy). It is derived from the popular but no longer maintained implementation at https://github.com/gpleiss/temperature_scaling/blob/master/temperature_scaling.py
with the following changes:
- using inverse temperatures to prevent division by zero errors
- using 50 optimizer steps instead of a single one (seemingly an error in the mentioned repo). (The original paper mentions that 10 CG iterations should be enough, here it is 50 L-BFGS iterations.)
- accepting probabilities as provided by many scikit-learn estimators using `predict_proba()`. The code converts probabilities to logits using `log(probs + 1e-10)`. While the logits are only determined up to a constant shift, the choice of the constant does not affect the result of temperature scaling.
```python
import torch
import torch.nn as nn
import numpy as np
from sklearn.base import BaseEstimator

class InverseTemperatureScalingCalibrator(BaseEstimator):
    # following https://github.com/gpleiss/temperature_scaling/blob/master/temperature_scaling.py
    def _get_logits(self, X):
        X = X + 1e-10
        X /= np.sum(X, axis=-1, keepdims=True)
        return torch.as_tensor(np.log(X), dtype=torch.float32)

    def fit(self, X, y):
        # X should be the probabilities as output by predict_proba()
        logits = self._get_logits(X)
        labels = torch.as_tensor(y)
        self.inv_temperature_ = nn.Parameter(torch.ones(1) / 1.5)
        criterion = nn.CrossEntropyLoss()

        optimizer = torch.optim.LBFGS([self.inv_temperature_], lr=0.01, max_iter=50)

        def eval():
            optimizer.zero_grad()
            y_pred = logits * self.inv_temperature_[:, None]
            loss = criterion(y_pred, labels)
            loss.backward()
            return loss

        for i in range(50):
            optimizer.step(eval)

        print(f'Optimal temperature: {(1./self.inv_temperature_).item():g}')
        return self

    def predict_proba(self, X):
        # X should be the probabilities as output by predict_proba()
        logits = self._get_logits(X)
        with torch.no_grad():
            y_pred = logits * self.inv_temperature_[:, None]
            return torch.softmax(y_pred, dim=-1).detach().numpy()
```

### Describe alternatives you've considered, if relevant

Centered isotonic regression (https://github.com/scikit-learn/scikit-learn/pull/21454) is less popular and does not fully solve the problem of affecting rank-based metrics.
Beta-calibration (https://github.com/scikit-learn/scikit-learn/issues/25552) seems very similar or even partially identical but is less well-cited, and only formulated for binary classification.

### Additional context

_No response_

COMMENT:
Hello,

I've implemented a simple version of temperature scaling with `numpy` and `scipy` using the `_temperature_scaling` function.

The function takes the following arguments:

- `predictions`: `(n_samples, n_classes)` array, the output of the `predict_proba` method.
- `labels`: `(n_samples, )` array, the correct label in ordinal fashion.
- `initial_temperature`: float, the initial temperature to start with (as seen in the OP, it's set to 1.5).

I've provided my proof of work below:
```python
import numpy as np
from scipy.optimize import minimize
from scipy.special import softmax

def _row_max_normalization(data: np.ndarray) -> np.ndarray:
    '''Normalise the output by subtracting
       the per-row maximum element.
    '''
    row_max: np.ndarray = np.max(data, 
                                 axis = 1, 
                                 keepdims = True
                                )
    
    return data - row_max

def _softmax_T(predictions: np.ndarray, 
               temperature: float,
              ) -> np.ndarray:
    '''Softmax function scaled by the
       inverse temperature.
    '''
    
    softmax_T_output: np.ndarray = predictions
    softmax_T_output = _row_max_normalization(softmax_T_output)  
    softmax_T_output /= temperature  
    softmax_T_output = softmax(softmax_T_output, 
                               axis = 1
                              )
    softmax_T_output = softmax_T_output.astype(dtype = predictions.dtype)
    
    return softmax_T_output

def _exp_T(predictions: np.ndarray, 
           temperature: float
          ) -> np.ndarray:
    '''Scale by inverse temperature,
       and then apply the nature
       exponential function
    '''
    
    exp_T_output: np.ndarray = predictions
    exp_T_output = _row_max_normalization(exp_T_output)
    exp_T_output /= temperature
    exp_T_output = np.exp(exp_T_output)
    
    return exp_T_output 

def _temperature_scaling(predictions: np.ndarray, 
                         labels: np.ndarray, 
                         initial_temperature: float
                        ) -> float:
    
    def negative_log_likelihood(temperature: float):
        '''Negative Log Likelihood Loss with respect
           to Temperature
        '''
        
        # Losses
        losses: np.ndarray = _softmax_T(predictions, 
                                          temperature
                                         )
            
        # Select the probability of the correct class
        losses = losses[np.arange(losses.shape[0]), 
                        labels
                       ]
        
        losses = np.log(losses)
        
        # Derivates with respect to Temperature
        exp_T: np.ndarray = _exp_T(predictions, temperature)
        exp_T_sum = exp_T.sum(axis = 1)
        
        term_1: np.ndarray = _row_max_normalization(predictions)
        term_1 /= temperature ** 2
        term_1 = - term_1[np.arange(term_1.shape[0]), 
                          labels
                         ]
        term_1 *= exp_T_sum
        
        
        
        term_2: np.ndarray = _row_max_normalization(predictions)
        term_2 /= temperature ** 2
        term_2 = _row_max_normalization(term_2)
        term_2 *= exp_T
        term_2 = term_2.sum(axis = 1)
        
        dL_dts: np.ndarray = (term_1 + term_2) / exp_T_sum
            
        # print(f"{-losses.sum() = },  {-dL_dts.sum() = }")
            
        return -losses.sum(),  -dL_dts.sum()
    
    temperature_minimizer: minimize = minimize(negative_log_likelihood, 
                                               initial_temperature, 
                                               method = "L-BFGS-B",
                                               jac = True,
                                               options = {"gtol": 1e-6,
                                                           "ftol": 64 * np.finfo(float).eps,
                                                         }
                                              )
        
    return temperature_minimizer.x[0]
```

I tested the function with the iris dataset, employing both the support vector classifier and the logistic regressor. The initial temperature of 1.5 converged to 0.15 and 0.12, respectively:
```python
>>> import numpy as np
>>> from sklearn import datasets
>>> from sklearn.model_selection import train_test_split
>>> from sklearn.svm import SVC
>>> from sklearn.linear_model import LogisticRegression
>>> X, y = datasets.load_iris(return_X_y=True)
>>> X_train, X_test, y_train, y_test = train_test_split(X, y)
>>> SVC_classifier = SVC(probability=True)
>>> SVC_classifier.fit(X_train, y_train)
>>> svc_predictions = SVC_classifier.predict_proba(X_train)
>>> Logistic_Regression = LogisticRegression()
>>> Logistic_Regression.fit(X_train, y_train)
>>> logistic_predictions = Logistic_Regression.predict_proba(X_train)
>>> _temperature_scaling(svc_predictions, y_train, 1.5)
0.1491147513643915
>>> _temperature_scaling(logistic_predictions, y_train, 1.5)
0.1197697499802383
```

Any comments are welcome! Hope this can bring something new to the table.

Labels: New Feature, Moderate, help wanted, module:calibration