Add partial_fit Functionality to LinearDiscriminantAnalysis Classifier

### Describe the workflow you want to enable

Currently, Scikit-learn's LinearDiscriminantAnalysis (LDA) classifier does not support incremental learning through the partial_fit method.   This poses challenges when processing large scale classification problems for which the full training set might not fit in memory. 

### Describe your proposed solution

Implementing partial_fit would allow users to train the LDA model incrementally, updating the model with batches of data as they become available.  This is consistent with the existing scikit-learn API, which currently supports this functionality for various other models such as the Naive Bayes classifiers.

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

It may prove to be difficult to support all possible solvers and functionality (e.g.., shrinkage).  As an alternative, an IncrementalLinearDiscriminantAnalysis classifier could be introduced that doesn't support all of the LInearDiscriminantAnalysis parameter options.

### Additional context

I've implemented a basic extension of the current LinearDiscriminantAnalysis classifier that accomplishes this, and am willing to do the development.

COMMENT:
Here's a version I've done some limited testing with.  I can create a github repo with code for the various test cases I've tested so far.   I've only experimented with incremental SVD up to this point.  (Edit:  I did a couple of quick tests with lsqr and eigen solvers, and those seem to work ok as well as long as shrinkage is "None")

```python
import numpy as np
import unittest
import logging
from sklearn.utils.validation import check_array
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.model_selection import train_test_split
from sklearn.datasets import (make_classification, load_wine,
                              make_multilabel_classification)
from sklearn.metrics import accuracy_score, classification_report
from sklearn.utils import shuffle
from scipy.sparse import csr_matrix

logging.basicConfig(level=logging.INFO, format='%(message)s')

class LinearDiscriminantAnalysis2(LinearDiscriminantAnalysis):
    """
    Extension of LinearDiscriminantAnalysis to add partial_fit functionality.
    Allows for incremental learning and the addition of new classes during partial_fit calls.

    **Note:** Currently, only the 'svd' solver is supported for incremental learning.
    Other solvers ('lsqr', 'eigen') are not supported and will raise an error if used.

    This class utilizes Welford's online algorithm for updating the mean and covariance incrementally.
    Reference: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm
    """
    def __init__(self, solver='svd', shrinkage=None, priors=None, n_components=None,
                 store_covariance=False, tol=1.0e-4):
        """
        Initialize the LinearDiscriminantAnalysis2.

        Parameters
        ----------
        solver : {'svd', 'lsqr', 'eigen'}, default='svd'
            Solver to use. Currently, only 'svd' is supported for incremental learning.
            'lsqr' and 'eigen' solvers are not supported and will raise an error if used.

        shrinkage : string, float or None, default=None
            Shrinkage parameter. **Currently unused**, as only 'svd' solver is supported.
            If using 'lsqr' or 'eigen', shrinkage would determine the regularization applied
            to the covariance estimates.

        priors : array-like of shape (n_classes,), default=None
            Prior probabilities of the classes. If specified, the priors are not
            adjusted according to the data.

        n_components : int or None, default=None
            Number of components (< n_classes - 1) for dimensionality reduction.
            **Currently unused**, as only 'svd' solver is supported.

        store_covariance : bool, default=False
            If True, the covariance matrix is stored in the `covariance_` attribute.
            **Currently unused**, as only 'svd' solver is supported.

        tol : float, default=1e-4
            Tolerance for rank estimation in 'eigen' and 'lsqr' solver.
            **Currently unused**, as only 'svd' solver is supported.
        """
        # Initialize the superclass with supported parameters
        super().__init__(solver=solver, shrinkage=shrinkage, priors=priors,
                         n_components=n_components, store_covariance=store_covariance,
                         tol=tol)
        
        # Error handling: Ensure only 'svd' solver is used
        if self.solver != 'svd':
            raise NotImplementedError(
                f"Solver '{self.solver}' is not supported in LinearDiscriminantAnalysis2. "
                "Only 'svd' solver is supported for incremental learning."
            )
        
        # Reset internal state
        self._reset()

    def _reset(self):
        """Reset the internal data-dependent attributes."""
        self.classes_ = None
        self._means = None
        self._n_samples_class = None
        self.coef_ = None
        self.intercept_ = None
        self._covariance_accumulator = None

    def partial_fit(self, X, y):
        """
        Incrementally fit the LDA model with new data. Allows for new classes to be added 
        even if they weren't present in the initial fit.

        This method updates the mean and covariance using Welford's online algorithm.
        Reference: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            Training data.

        y : array-like of shape (n_samples,)
            Target values.

        Returns
        -------
        self : object
            Fitted estimator.
        """
        # Validate and convert input data
        X, y = check_array(X), np.asarray(y)

        if self.classes_ is None:
            # First partial_fit call: initialize classes and statistics
            self.classes_ = np.unique(y)
            n_classes = len(self.classes_)
            n_features = X.shape[1]

            # Initialize means, covariance accumulators, and sample counts
            self._means = np.zeros((n_classes, n_features))
            self._covariance_accumulator = np.zeros((n_classes, n_features, n_features))
            self._n_samples_class = np.zeros(n_classes)
        else:
            current_classes = self.classes_

            # Identify unique classes in the current batch of data
            unique_y = np.unique(y)

            # Determine new classes that are not already in self.classes_
            new_classes = np.setdiff1d(unique_y, current_classes)

            if new_classes.size > 0:
                # Append new classes to self.classes_
                self.classes_ = np.concatenate([self.classes_, new_classes])
                n_new = len(new_classes)
                n_features = X.shape[1]

                # Expand _means to accommodate new classes
                self._means = np.vstack([self._means, np.zeros((n_new, n_features))])

                # Expand _covariance_accumulator for new classes
                new_cov_acc = np.zeros((n_new, n_features, n_features))
                self._covariance_accumulator = np.concatenate([self._covariance_accumulator, new_cov_acc], axis=0)

                # Expand _n_samples_class for new classes
                self._n_samples_class = np.concatenate([self._n_samples_class, np.zeros(n_new)])

        # Compute unique classes and their respective sample counts in the current batch
        unique_y, counts_y = np.unique(y, return_counts=True)

        # Iterate over each class present in the current batch
        for cls, count in zip(unique_y, counts_y):
            # Find the index of the current class in self.classes_
            idx_class = np.where(self.classes_ == cls)[0][0]

            # Extract all samples belonging to the current class
            X_c = X[y == cls]

            # Compute the mean of the new samples for the current class
            x_mean_new = X_c.mean(axis=0)

            # Center the new samples by subtracting the new mean
            centered_X_new = X_c - x_mean_new

            # Compute the covariance accumulator for the new samples
            S_new = np.dot(centered_X_new.T, centered_X_new)

            # Retrieve the existing number of samples, mean, and covariance accumulator for the class
            n_c = self._n_samples_class[idx_class]
            mu_c = self._means[idx_class]
            S_c = self._covariance_accumulator[idx_class]

            if n_c == 0:
                # First time seeing this class: initialize mean and covariance accumulator
                mu_c_new = x_mean_new
                S_c_new = S_new
            else:
                # Update mean and covariance accumulator using Welford's online algorithm
                delta = x_mean_new - mu_c
                mu_c_new = (n_c * mu_c + count * x_mean_new) / (n_c + count)
                S_c_new = S_c + S_new + (n_c * count) / (n_c + count) * np.outer(delta, delta)

            # Update the class-specific mean, covariance accumulator, and sample count
            self._means[idx_class] = mu_c_new
            self._covariance_accumulator[idx_class] = S_c_new
            self._n_samples_class[idx_class] = n_c + count

        # Total number of samples across all classes
        total_samples = np.sum(self._n_samples_class)

        # Update class prior probabilities based on the accumulated sample counts
        self.priors_ = self._n_samples_class / total_samples

        # Calculate the degrees of freedom for covariance computation
        n_total_samples_minus_classes = total_samples - len(self.classes_)
        if n_total_samples_minus_classes <= 0:
            # Avoid division by zero; defer covariance computation until more data is available
            return self

        # Compute the pooled covariance matrix by summing accumulators and normalizing
        self.covariance_ = np.sum(self._covariance_accumulator, axis=0) / n_total_samples_minus_classes

        # Update the inverse of the covariance matrix, handling potential singularity
        inv_covariance = np.linalg.pinv(self.covariance_)

        # Update the coefficients (weights) for each class
        self.coef_ = np.dot(inv_covariance, self._means.T).T

        # Compute the intercepts for each class
        with np.errstate(divide='ignore'):
            self.intercept_ = -0.5 * np.sum(self._means * self.coef_, axis=1) + np.log(self.priors_)

        return self

    def fit(self, X, y):
        """
        Complete fit using accumulated data.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            Training data.

        y : array-like of shape (n_samples,)
            Target values.

        Returns
        -------
        self : object
            Fitted estimator.
        """
        self._reset()
        self.partial_fit(X, y)
        return self

```

Labels: New Feature