Repository: scikit-learn/scikit-learn
Pull request: FEA Allow passing custom covariance estimators to LDA
Fixes #14436 
Plug-in any covariance estimator into Fisher discriminants
This PR allows the user to plug-in any sklearn covariance estimator into LDA and QDA. This option is only possible with the "lsqr" and "eigen" solver of LDA and for the "lsqr" solver of QDA (initially only an "svd" solver was available for QDA so I made this explicit and added the "lsqr").

Problem statement:
Linear and Quadratic discriminant analysis could be more modular to allow for user specific estimation of covariance
By default sample estimate is used to estimate the covariance in Fisher discriminant such as LDA and QDA. However it could be useful to make the API a bit more general so that any estimator of covariance can be plugged in (I know that Ledoit-Wolf is already implemented). 

What I propose is to have two additional parameters in the API:
covariance_estimator and standardize
- covariance estimator would be a sklearn-like estimator of covariance such as OAS or LedoitWolf
- standardize tells you if the data should be normalized for estimation of covariance (current behaviour when shrinkage="auto")

I have a pull request ready for this. Tell me if you think this is of any use.
Hugo

Hints from issue discussion:
@hugorichard if you have a ready PR, I'd say submit and we can see how we feel about it then.

Commit message: correct docs issues

Patch to review:
diff --git a/sklearn/discriminant_analysis.py b/sklearn/discriminant_analysis.py
index 7a0bdcf398aa2..3438c6ce9b2bb 100644
--- a/sklearn/discriminant_analysis.py
+++ b/sklearn/discriminant_analysis.py
@@ -25,47 +25,63 @@
 from .utils.extmath import softmax
 from .preprocessing import StandardScaler
 
-
 __all__ = ['LinearDiscriminantAnalysis', 'QuadraticDiscriminantAnalysis']
 
 
-def _cov(X, shrinkage=None):
-    """Estimate covariance matrix (using optional shrinkage).
-
+def _cov(X, shrinkage=None, covariance_estimator=None, standardize=False):
+    """Estimate covariance matrix (using covariance_estimator).
     Parameters
     ----------
     X : array-like, shape (n_samples, n_features)
         Input data.
-
+    covariance_estimator : BaseEstimator
+        A covariance estimator, possible values
+            - None: shrinkage parameter drives the estimate
+            - BaseEstimator with a fit method and
+            a .covariance_ attribute like the estimators
+            in sklearn.covariance
+    standardize : bool
+        if True standardize the data before applying the covariance estimator
+        Note: only used when covariance_estimator is not None
     shrinkage : string or float, optional
         Shrinkage parameter, possible values:
           - None or 'empirical': no shrinkage (default).
           - 'auto': automatic shrinkage using the Ledoit-Wolf lemma.
           - float between 0 and 1: fixed shrinkage parameter.
-
     Returns
     -------
     s : array, shape (n_features, n_features)
         Estimated covariance matrix.
     """
-    shrinkage = "empirical" if shrinkage is None else shrinkage
-    if isinstance(shrinkage, str):
-        if shrinkage == 'auto':
-            sc = StandardScaler()  # standardize features
+    if covariance_estimator is None:
+        shrinkage = "empirical" if shrinkage is None else shrinkage
+        if isinstance(shrinkage, str):
+            if shrinkage == 'auto':
+                sc = StandardScaler()  # standardize features
+                X = sc.fit_transform(X)
+                s = ledoit_wolf(X)[0]
+                # rescale
+                s = sc.scale_[:, np.newaxis] * s * sc.scale_[np.newaxis, :]
+            elif shrinkage == 'empirical':
+                s = empirical_covariance(X)
+            else:
+                raise ValueError('unknown shrinkage parameter')
+        elif isinstance(shrinkage, float) or isinstance(shrinkage, int):
+            if shrinkage < 0 or shrinkage > 1:
+                raise ValueError('shrinkage parameter must be between 0 and 1')
+            s = shrunk_covariance(empirical_covariance(X), shrinkage)
+        else:
+            raise TypeError('shrinkage must be of string or int type')
+    else:
+        if standardize:
+            sc = StandardScaler()
             X = sc.fit_transform(X)
-            s = ledoit_wolf(X)[0]
-            # rescale
+            covariance_estimator.fit(X)
+            s = covariance_estimator.covariance_
             s = sc.scale_[:, np.newaxis] * s * sc.scale_[np.newaxis, :]
-        elif shrinkage == 'empirical':
-            s = empirical_covariance(X)
         else:
-            raise ValueError('unknown shrinkage parameter')
-    elif isinstance(shrinkage, float) or isinstance(shrinkage, int):
-        if shrinkage < 0 or shrinkage > 1:
-            raise ValueError('shrinkage parameter must be between 0 and 1')
-        s = shrunk_covariance(empirical_covariance(X), shrinkage)
-    else:
-        raise TypeError('shrinkage must be of string or int type')
+            covariance_estimator.fit(X)
+            s = covariance_estimator.covariance_
     return s
 
 
@@ -93,26 +109,37 @@ def _class_means(X, y):
     return means
 
 
-def _class_cov(X, y, priors, shrinkage=None):
+def _class_cov(
+        X,
+        y,
+        priors=None,
+        shrinkage=None,
+        covariance_estimator=None,
+        standardize=False
+):
     """Compute class covariance matrix.
-
     Parameters
     ----------
     X : array-like, shape (n_samples, n_features)
         Input data.
-
     y : array-like, shape (n_samples,) or (n_samples, n_targets)
         Target values.
-
     priors : array-like, shape (n_classes,)
         Class priors.
-
+    covariance_estimator : BaseEstimator
+        A covariance estimator, possible values
+            - None: shrinkage parameter drives the estimate
+            - BaseEstimator with a fit method and
+            a .covariance_ attribute like the estimators
+            in sklearn.covariance
+    standardize : bool
+        if True standardize the data before applying the covariance estimator
+        Note: only used when covariance_estimator is not None
     shrinkage : string or float, optional
         Shrinkage parameter, possible values:
           - None: no shrinkage (default).
           - 'auto': automatic shrinkage using the Ledoit-Wolf lemma.
           - float between 0 and 1: fixed shrinkage parameter.
-
     Returns
     -------
     cov : array-like, shape (n_features, n_features)
@@ -122,10 +149,64 @@ def _class_cov(X, y, priors, shrinkage=None):
     cov = np.zeros(shape=(X.shape[1], X.shape[1]))
     for idx, group in enumerate(classes):
         Xg = X[y == group, :]
-        cov += priors[idx] * np.atleast_2d(_cov(Xg, shrinkage))
+        cov += priors[idx] * np.atleast_2d(_cov(
+            Xg,
+            shrinkage,
+            covariance_estimator,
+            standardize))
     return cov
 
 
+def _classes_cov(
+        X,
+        y,
+        shrinkage=None,
+        covariance_estimator=None,
+        standardize=False
+):
+    """Compute covariance matrix of each class.
+    Parameters
+    ----------
+    X : array-like, shape (n_samples, n_features)
+        Input data.
+    y : array-like, shape (n_samples,) or (n_samples, n_targets)
+        Target values.
+    covariance_estimator : BaseEstimator
+        A covariance estimator, possible values
+            - None: shrinkage parameter drives the estimate
+            - BaseEstimator with a fit method and
+            a .covariance_ attribute like the estimators
+            in sklearn.covariance
+    standardize : bool
+        if True standardize the data before applying the covariance estimator
+        Note: only used when covariance_estimator is not None
+    shrinkage : string or float, optional
+        Shrinkage parameter, possible values:
+          - None: no shrinkage (default).
+          - 'auto': automatic shrinkage using the Ledoit-Wolf lemma.
+          - float between 0 and 1: fixed shrinkage parameter.
+    Returns
+    -------
+    cov : n_classes list of array-like shape (n_features, n_features)
+        Class covariance matrix of each class.
+    """
+    classes = np.unique(y)
+    covs = []
+    for idx, group in enumerate(classes):
+        Xg = X[y == group, :]
+        if len(Xg) == 1:
+            raise ValueError(
+                'y contains a class with'
+                'only 1 sample. Covariance '
+                'is ill defined.')
+        covs.append(np.atleast_2d(_cov(
+            Xg,
+            shrinkage,
+            covariance_estimator,
+            standardize)))
+    return np.array(covs)
+
+
 class LinearDiscriminantAnalysis(BaseEstimator, LinearClassifierMixin,
                                  TransformerMixin):
     """Linear Discriminant Analysis
@@ -162,6 +243,17 @@ class LinearDiscriminantAnalysis(BaseEstimator, LinearClassifierMixin,
 
         Note that shrinkage works only with 'lsqr' and 'eigen' solvers.
 
+    covariance_estimator : BaseEstimator
+        A covariance estimator, possible values
+            - None: shrinkage parameter drives the estimate
+            - BaseEstimator with a fit method and
+            a .covariance_ attribute like the estimators
+            in sklearn.covariance
+
+    standardize : bool
+        if True standardize the data before applying the covariance estimator
+        Note: only used when covariance_estimator is not None
+
     priors : array, optional, shape (n_classes,)
         Class priors.
 
@@ -247,16 +339,27 @@ class LinearDiscriminantAnalysis(BaseEstimator, LinearClassifierMixin,
     [1]
     """
 
-    def __init__(self, solver='svd', shrinkage=None, priors=None,
-                 n_components=None, store_covariance=False, tol=1e-4):
+    def __init__(
+            self,
+            solver='svd',
+            shrinkage=None,
+            covariance_estimator=None,
+            standardize=False,
+            priors=None,
+            n_components=None,
+            store_covariance=False,
+            tol=1e-4
+    ):
         self.solver = solver
         self.shrinkage = shrinkage
         self.priors = priors
         self.n_components = n_components
         self.store_covariance = store_covariance  # used only in svd solver
         self.tol = tol  # used only in svd solver
+        self.covariance_estimator = covariance_estimator
+        self.standardize = standardize
 
-    def _solve_lsqr(self, X, y, shrinkage):
+    def _solve_lsqr(self, X, y, shrinkage, covariance_estimator, standardize):
         """Least squares solver.
 
         The least squares solver computes a straightforward solution of the
@@ -279,7 +382,19 @@ def _solve_lsqr(self, X, y, shrinkage):
               - 'auto': automatic shrinkage using the Ledoit-Wolf lemma.
               - float between 0 and 1: fixed shrinkage parameter.
 
-        Notes
+        covariance_estimator : BaseEstimator
+            A covariance estimator, possible values
+                - None: shrinkage parameter drives the estimate
+                - BaseEstimator with a fit method and
+                a .covariance_ attribute like the estimators
+                in sklearn.covariance
+
+        standardize : bool
+            if True standardize the data before applying the
+            covariance estimator
+            Note: only used when covariance_estimator is not None
+
+       Notes
         -----
         This solver is based on [1]_, section 2.6.2, pp. 39-41.
 
@@ -290,12 +405,20 @@ def _solve_lsqr(self, X, y, shrinkage):
            0-471-05669-3.
         """
         self.means_ = _class_means(X, y)
-        self.covariance_ = _class_cov(X, y, self.priors_, shrinkage)
+        self.covariance_ = _class_cov(
+                X,
+                y,
+                self.priors_,
+                shrinkage,
+                covariance_estimator,
+                standardize)
         self.coef_ = linalg.lstsq(self.covariance_, self.means_.T)[0].T
         self.intercept_ = (-0.5 * np.diag(np.dot(self.means_, self.coef_.T)) +
                            np.log(self.priors_))
 
-    def _solve_eigen(self, X, y, shrinkage):
+    def _solve_eigen(
+            self, X, y, shrinkage,
+            covariance_estimator, standardize):
         """Eigenvalue solver.
 
         The eigenvalue solver computes the optimal solution of the Rayleigh
@@ -317,6 +440,18 @@ class scatter). This solver supports both classification and
               - 'auto': automatic shrinkage using the Ledoit-Wolf lemma.
               - float between 0 and 1: fixed shrinkage constant.
 
+        covariance_estimator : BaseEstimator
+            A covariance estimator, possible values
+                - None: shrinkage parameter drives the estimate
+                - BaseEstimator with a fit method and
+                a .covariance_ attribute like the estimators
+                in sklearn.covariance
+
+        standardize : bool
+            if True standardize the data before applying
+            the covariance estimator
+            Note: only used when covariance_estimator is not None
+
         Notes
         -----
         This solver is based on [1]_, section 3.8.3, pp. 121-124.
@@ -328,10 +463,20 @@ class scatter). This solver supports both classification and
            0-471-05669-3.
         """
         self.means_ = _class_means(X, y)
-        self.covariance_ = _class_cov(X, y, self.priors_, shrinkage)
+        self.covariance_ = _class_cov(
+                X,
+                y,
+                self.priors_,
+                shrinkage,
+                covariance_estimator,
+                standardize)
 
         Sw = self.covariance_  # within scatter
-        St = _cov(X, shrinkage)  # total scatter
+        St = _cov(
+                X,
+                shrinkage,
+                covariance_estimator,
+                standardize)  # total scatter
         Sb = St - Sw  # between scatter
 
         evals, evecs = linalg.eigh(Sb, Sw)
@@ -476,11 +621,26 @@ def fit(self, X, y):
         if self.solver == 'svd':
             if self.shrinkage is not None:
                 raise NotImplementedError('shrinkage not supported')
+            if self.covariance_estimator is not None:
+                raise NotImplementedError(
+                        'covariance estimator'
+                        'is not supported'
+                        'with svd solver. Try an other solver')
             self._solve_svd(X, y)
         elif self.solver == 'lsqr':
-            self._solve_lsqr(X, y, shrinkage=self.shrinkage)
+            self._solve_lsqr(
+                    X,
+                    y,
+                    shrinkage=self.shrinkage,
+                    covariance_estimator=self.covariance_estimator,
+                    standardize=self.standardize)
         elif self.solver == 'eigen':
-            self._solve_eigen(X, y, shrinkage=self.shrinkage)
+            self._solve_eigen(
+                    X,
+                    y,
+                    shrinkage=self.shrinkage,
+                    covariance_estimator=self.covariance_estimator,
+                    standardize=self.standardize)
         else:
             raise ValueError("unknown solver {} (valid solvers are 'svd', "
                              "'lsqr', and 'eigen').".format(self.solver))
@@ -571,6 +731,26 @@ class QuadraticDiscriminantAnalysis(BaseEstimator, ClassifierMixin):
 
     Parameters
     ----------
+    solver : string, optional
+        Solver to use, possible values:
+          - 'svd': Singular value decomposition (default).
+            Does not compute the covariance matrix, therefore this solver is
+            recommended for data with a large number of features.
+          - 'lsqr': Least squares solution, can be combined with
+          custom covariance estimator.
+
+    covariance_estimator : BaseEstimator
+        A covariance estimator, possible values
+            - None: shrinkage parameter drives the estimate
+            - BaseEstimator with a fit method and
+            a .covariance_ attribute like the estimators
+            in sklearn.covariance
+        Note: does not work when solver='svd'
+
+    standardize : bool
+        if True standardize the data before applying the covariance estimator
+        Note: only used when covariance_estimator is not None
+
     priors : array, optional, shape = [n_classes]
         Priors on classes
 
@@ -605,11 +785,13 @@ class QuadraticDiscriminantAnalysis(BaseEstimator, ClassifierMixin):
         ``n_k = min(n_features, number of elements in class k)``
         It is the rotation of the Gaussian distribution, i.e. its
         principal axis.
+        Note: only available when solver="svd"
 
     scalings_ : list of arrays
         For each class k an array of shape [n_k]. It contains the scaling
         of the Gaussian distributions along its principal axes, i.e. the
         variance in the rotated coordinate system.
+        Note: only available when solver="svd"
 
     Examples
     --------
@@ -629,22 +811,25 @@ class QuadraticDiscriminantAnalysis(BaseEstimator, ClassifierMixin):
         Discriminant Analysis
     """
 
-    def __init__(self, priors=None, reg_param=0., store_covariance=False,
-                 tol=1.0e-4):
+    def __init__(
+            self,
+            solver='svd',
+            covariance_estimator=None,
+            standardize=False,
+            priors=None,
+            reg_param=0.,
+            store_covariance=False,
+            tol=1.0e-4):
         self.priors = np.asarray(priors) if priors is not None else None
         self.reg_param = reg_param
         self.store_covariance = store_covariance
         self.tol = tol
+        self.solver = solver
+        self.covariance_estimator = covariance_estimator
+        self.standardize = standardize
 
-    def fit(self, X, y):
-        """Fit the model according to the given training data and parameters.
-
-            .. versionchanged:: 0.19
-               ``store_covariances`` has been moved to main constructor as
-               ``store_covariance``
-
-            .. versionchanged:: 0.19
-               ``tol`` has been moved to main constructor.
+    def _solve_svd(self, X, y):
+        """SVD solver
 
         Parameters
         ----------
@@ -655,18 +840,8 @@ def fit(self, X, y):
         y : array, shape = [n_samples]
             Target values (integers)
         """
-        X, y = check_X_y(X, y)
-        check_classification_targets(y)
-        self.classes_, y = np.unique(y, return_inverse=True)
         n_samples, n_features = X.shape
         n_classes = len(self.classes_)
-        if n_classes < 2:
-            raise ValueError('The number of classes has to be greater than'
-                             ' one; got %d class' % (n_classes))
-        if self.priors is None:
-            self.priors_ = np.bincount(y) / float(n_samples)
-        else:
-            self.priors_ = self.priors
 
         cov = None
         store_covariance = self.store_covariance
@@ -695,27 +870,173 @@ def fit(self, X, y):
                 cov.append(np.dot(S2 * Vt.T, Vt))
             scalings.append(S2)
             rotations.append(Vt.T)
+
         if self.store_covariance or store_covariance:
             self.covariance_ = cov
         self.means_ = np.asarray(means)
         self.scalings_ = scalings
         self.rotations_ = rotations
+
+        return self
+
+    def _solve_lsqr(
+            self,
+            X,
+            y,
+            reg_param,
+            covariance_estimator,
+            standardize):
+        """Least squares solver.
+
+        The least squares solver computes a straightforward solution of the
+        optimal decision rule based directly on the discriminant functions.
+        Parameters
+        ----------
+        X : array-like, shape (n_samples, n_features)
+            Training data.
+
+        y : array-like, shape (n_samples,) or (n_samples, n_classes)
+            Target values.
+
+        reg_param : float, optional
+            Regularizes the covariance estimate as
+            ``(1-reg_param)*Sigma + reg_param*np.eye(n_features)``
+
+        covariance_estimator : BaseEstimator
+            A covariance estimator, possible values
+                - None: shrinkage parameter drives the estimate
+                - BaseEstimator with a fit method and
+                a .covariance_ attribute like the estimators
+                in sklearn.covariance
+            Note: does not work when solver='svd'
+
+        standardize : bool
+            if True standardize the data before applying
+            the covariance estimator
+            Note: only used when covariance_estimator is not None
+
+       Notes
+        -----
+        This solver is based on [1]_, section 2.6.2, pp. 39-41.
+
+        References
+        ----------
+        .. [1] R. O. Duda, P. E. Hart, D. G. Stork. Pattern Classification
+           (Second Edition). John Wiley & Sons, Inc., New York, 2001. ISBN
+           0-471-05669-3.
+        """
+        self.means_ = _class_means(X, y)
+        self.covariance_ = _classes_cov(
+                X,
+                y,
+                reg_param,
+                covariance_estimator,
+                standardize)
+        return self
+
+    def fit(self, X, y):
+        """Fit the model according to the given training data and parameters.
+
+            .. versionchanged:: 0.19
+               ``store_covariances`` has been moved to main constructor as
+               ``store_covariance``
+
+            .. versionchanged:: 0.19
+               ``tol`` has been moved to main constructor.
+
+        Parameters
+        ----------
+        X : array-like, shape = [n_samples, n_features]
+            Training vector, where n_samples is the number of samples and
+            n_features is the number of features.
+
+        y : array, shape = [n_samples]
+            Target values (integers)
+        """
+        # FIXME: Future warning to be removed in 0.23
+        X, y = check_X_y(X, y, ensure_min_samples=2, estimator=self,
+                         dtype=[np.float64, np.float32])
+        check_classification_targets(y)
+        self.classes_, y = np.unique(y, return_inverse=True)
+        n_samples, _ = X.shape
+        n_classes = len(self.classes_)
+
+        if n_classes < 2:
+            raise ValueError('The number of classes has to be greater than'
+                             ' one; got %d class' % (n_classes))
+        if n_samples == n_classes:
+            raise ValueError("The number of samples must be more "
+                             "than the number of classes.")
+
+        if self.priors is None:  # estimate priors from sample
+            _, y_t = np.unique(y, return_inverse=True)  # non-negative ints
+            self.priors_ = np.bincount(y_t) / float(len(y))
+        else:
+            self.priors_ = np.asarray(self.priors)
+
+        if (self.priors_ < 0).any():
+            raise ValueError("priors must be non-negative")
+        if not np.isclose(self.priors_.sum(), 1.0):
+            warnings.warn("The priors do not sum to 1. Renormalizing",
+                          UserWarning)
+            self.priors_ = self.priors_ / self.priors_.sum()
+
+        if self.solver == 'svd':
+            if self.covariance_estimator is not None:
+                raise NotImplementedError(
+                    'covariance estimator is not'
+                    'supported with svd solver. Try an other solver')
+            self._solve_svd(X, y)
+        elif self.solver == 'lsqr':
+            self._solve_lsqr(
+                    X,
+                    y,
+                    reg_param=self.reg_param,
+                    covariance_estimator=self.covariance_estimator,
+                    standardize=self.standardize)
+        else:
+            raise ValueError("unknown solver {} (valid solvers are 'svd' "
+                             "and 'lsqr').".format(self.solver))
         return self
 
     def _decision_function(self, X):
         check_is_fitted(self, 'classes_')
-
         X = check_array(X)
-        norm2 = []
-        for i in range(len(self.classes_)):
-            R = self.rotations_[i]
-            S = self.scalings_[i]
-            Xm = X - self.means_[i]
-            X2 = np.dot(Xm, R * (S ** (-0.5)))
-            norm2.append(np.sum(X2 ** 2, 1))
-        norm2 = np.array(norm2).T  # shape = [len(X), n_classes]
-        u = np.asarray([np.sum(np.log(s)) for s in self.scalings_])
-        return (-0.5 * (norm2 + u) + np.log(self.priors_))
+
+        if self.solver == "svd":
+            norm2 = []
+            for i in range(len(self.classes_)):
+                R = self.rotations_[i]
+                S = self.scalings_[i]
+                Xm = X - self.means_[i]
+                X2 = np.dot(Xm, R * (S ** (-0.5)))
+                norm2.append(np.sum(X2 ** 2, 1))
+            norm2 = np.array(norm2).T  # shape = [len(X), n_classes]
+            u = np.asarray([np.sum(np.log(s)) for s in self.scalings_])
+            return -0.5 * (norm2 + u) + np.log(self.priors_)
+
+        if self.solver == "lsqr":
+            n_classes = len(self.classes_)
+            dec_func = []
+            for k in range(n_classes):
+                mean = self.means_[k]
+                covariance = self.covariance_[k]
+                centeredX = X - mean
+                sign, log_det_cov_k = np.linalg.slogdet(covariance)
+                if sign <= 0:
+                    raise ValueError(
+                        "Determinant of covariance matrix"
+                        " corresponding to class %s was"
+                        " <= 0" % self.classes_[k])
+                log_prior = np.log(self.priors_[k])
+                quadratic_part = np.sum(centeredX * linalg.lstsq(
+                    covariance,
+                    centeredX.T)[0].T, axis=1)
+                decision_function_k = - 0.5 * quadratic_part
+                decision_function_k += - 0.5 * log_det_cov_k
+                decision_function_k += log_prior
+                dec_func.append(decision_function_k)
+            return np.array(dec_func).T
 
     def decision_function(self, X):
         """Apply decision function to an array of samples.
@@ -751,6 +1072,7 @@ def predict(self, X):
         -------
         C : array, shape = [n_samples]
         """
+        check_is_fitted(self, '_decision_function')
         d = self._decision_function(X)
         y_pred = self.classes_.take(d.argmax(1))
         return y_pred
@@ -768,6 +1090,7 @@ def predict_proba(self, X):
         C : array, shape = [n_samples, n_classes]
             Posterior probabilities of classification per class.
         """
+        check_is_fitted(self, '_decision_function')
         values = self._decision_function(X)
         # compute the likelihood of the underlying gaussian models
         # up to a multiplicative constant.
diff --git a/sklearn/tests/test_discriminant_analysis.py b/sklearn/tests/test_discriminant_analysis.py
index 6a32c4ec15058..8f9c4148ff2d9 100644
--- a/sklearn/tests/test_discriminant_analysis.py
+++ b/sklearn/tests/test_discriminant_analysis.py
@@ -22,6 +22,9 @@
 from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
 from sklearn.discriminant_analysis import _cov
 
+from sklearn.covariance import ShrunkCovariance
+from sklearn.covariance import LedoitWolf
+from sklearn.covariance import EmpiricalCovariance
 
 # Data is just 6 separable points in the plane
 X = np.array([[-2, -1], [-1, -1], [-1, -2], [1, 1], [1, 2], [2, 1]], dtype='f')
@@ -48,6 +51,10 @@
 # One element class
 y4 = np.array([1, 1, 1, 1, 1, 1, 1, 1, 2])
 
+# n_samples == n_classes
+X8 = np.random.rand(5, 2)
+y8 = np.array([0, 1, 2, 3, 4])
+
 # Data with less samples in a class than n_features
 X5 = np.c_[np.arange(8), np.zeros((8, 3))]
 y5 = np.array([0, 0, 0, 0, 0, 1, 1, 1])
@@ -91,9 +98,15 @@ def test_lda_predict():
     assert_raises(ValueError, clf.fit, X, y)
     clf = LinearDiscriminantAnalysis(solver="svd", shrinkage="auto")
     assert_raises(NotImplementedError, clf.fit, X, y)
+    clf = LinearDiscriminantAnalysis(solver="lsqr", shrinkage=np.array([1, 2]))
+    assert_raises(TypeError, clf.fit, X, y)
     # Test unknown solver
     clf = LinearDiscriminantAnalysis(solver="dummy")
     assert_raises(ValueError, clf.fit, X, y)
+    clf = LinearDiscriminantAnalysis(
+            solver="svd",
+            covariance_estimator=LedoitWolf())
+    assert_raises(NotImplementedError, clf.fit, X, y)
 
 
 @pytest.mark.parametrize("n_classes", [2, 3])
@@ -330,6 +343,53 @@ def test_lda_store_covariance():
     )
 
 
+def test_lda_shrinkage():
+    """
+    Test that shrunk covariance estimator and
+    shrinkage parameter behave the same
+    """
+    X = np.random.rand(100, 10)
+    Y = np.random.randint(3, size=(100))
+    c1 = LinearDiscriminantAnalysis(
+            store_covariance=True,
+            shrinkage=0.5,
+            solver="lsqr")
+    c2 = LinearDiscriminantAnalysis(
+            store_covariance=True,
+            covariance_estimator=ShrunkCovariance(shrinkage=0.5),
+            solver="lsqr")
+    c1.fit(X, Y)
+    c2.fit(X, Y)
+    for i in range(2):
+        assert_array_almost_equal(c1.means_[i], c2.means_[i])
+    assert_array_almost_equal(c1.covariance_, c2.covariance_)
+
+
+def test_lda_ledoitwolf():
+    """
+    When shrinkage="auto" current implementation
+    uses ledoitwolf estimation of
+    covariance after standardizing the data.
+    This checks that it is indeed the case
+    """
+    X = np.random.rand(100, 10)
+    Y = np.random.randint(3, size=(100,))
+    c1 = LinearDiscriminantAnalysis(
+            store_covariance=True,
+            shrinkage="auto",
+            solver="lsqr")
+    c2 = LinearDiscriminantAnalysis(
+            store_covariance=True,
+            covariance_estimator=LedoitWolf(),
+            solver="lsqr",
+            standardize=True)
+    c1.fit(X, Y)
+    c2.fit(X, Y)
+    for i in range(2):
+        assert_array_almost_equal(c1.means_[i], c2.means_[i])
+    assert_array_almost_equal(c1.covariance_, c2.covariance_)
+
+
 @pytest.mark.parametrize('n_features', [3, 5])
 @pytest.mark.parametrize('n_classes', [5, 3])
 def test_lda_dimension_warning(n_classes, n_features):
@@ -392,11 +452,12 @@ def test_lda_numeric_consistency_float32_float64():
         assert_allclose(clf_32.coef_, clf_64.coef_, rtol=rtol)
 
 
-def test_qda():
+@pytest.mark.parametrize("solver", ["svd", "lsqr"])
+def test_qda(solver):
     # QDA classification.
     # This checks that QDA implements fit and predict and returns
     # correct values for a simple toy dataset.
-    clf = QuadraticDiscriminantAnalysis()
+    clf = QuadraticDiscriminantAnalysis(solver=solver)
     y_pred = clf.fit(X6, y6).predict(X6)
     assert_array_equal(y_pred, y6)
 
@@ -417,27 +478,56 @@ def test_qda():
     # Classes should have at least 2 elements
     assert_raises(ValueError, clf.fit, X6, y4)
 
+    # Number of sample should be bigger than number of classes
+    assert_raises(ValueError, clf.fit, X8, y8)
 
-def test_qda_priors():
-    clf = QuadraticDiscriminantAnalysis()
+
+@pytest.mark.parametrize("solver", ["svd", "lsqr"])
+def test_qda_priors(solver):
+    clf = QuadraticDiscriminantAnalysis(solver=solver)
     y_pred = clf.fit(X6, y6).predict(X6)
     n_pos = np.sum(y_pred == 2)
 
-    neg = 1e-10
-    clf = QuadraticDiscriminantAnalysis(priors=np.array([neg, 1 - neg]))
+    neg = 1e-15
+    clf = QuadraticDiscriminantAnalysis(
+            solver=solver,
+            priors=np.array([neg, 1 - neg])
+    )
     y_pred = clf.fit(X6, y6).predict(X6)
     n_pos2 = np.sum(y_pred == 2)
 
     assert n_pos2 > n_pos
 
+    # Test that negative priors raise an error
+    priors = np.array([0.5, -0.5])
+    clf = QuadraticDiscriminantAnalysis(solver=solver, priors=priors)
+    msg = "priors must be non-negative"
+    assert_raise_message(ValueError, msg, clf.fit, X, y)
+
+    # Test that priors passed as a list are correctly handled (run to see if
+    # failure)
+    clf = QuadraticDiscriminantAnalysis(solver=solver, priors=[0.5, 0.5])
+    clf.fit(X, y)
+
+    # Test that priors always sum to 1
+    priors = np.array([0.5, 0.6])
+    prior_norm = np.array([0.45, 0.55])
+    clf = QuadraticDiscriminantAnalysis(solver=solver, priors=priors)
+    assert_warns(UserWarning, clf.fit, X, y)
+    assert_array_almost_equal(clf.priors_, prior_norm, 2)
+
 
-def test_qda_store_covariance():
+@pytest.mark.parametrize("solver", ["svd"])
+def test_qda_store_covariance(solver):
     # The default is to not set the covariances_ attribute
-    clf = QuadraticDiscriminantAnalysis().fit(X6, y6)
+    clf = QuadraticDiscriminantAnalysis(solver=solver).fit(X6, y6)
     assert not hasattr(clf, 'covariance_')
 
     # Test the actual attribute:
-    clf = QuadraticDiscriminantAnalysis(store_covariance=True).fit(X6, y6)
+    clf = QuadraticDiscriminantAnalysis(
+            solver=solver,
+            store_covariance=True
+    ).fit(X6, y6)
     assert hasattr(clf, 'covariance_')
 
     assert_array_almost_equal(
@@ -451,7 +541,8 @@ def test_qda_store_covariance():
     )
 
 
-def test_qda_regularization():
+@pytest.mark.parametrize("solver", ["svd", "lsqr"])
+def test_qda_regularization(solver):
     # the default is reg_param=0. and will cause issues
     # when there is a constant variable
     clf = QuadraticDiscriminantAnalysis()
@@ -474,6 +565,68 @@ def test_qda_regularization():
     assert_array_equal(y_pred5, y5)
 
 
+def test_qda_ledoitwolf():
+    # test bad solver
+    clf = QuadraticDiscriminantAnalysis(
+            solver="svd",
+            covariance_estimator=LedoitWolf()
+    )
+    assert_raises(NotImplementedError, clf.fit, X, y)
+
+    # test with unkwown solver
+    clf = QuadraticDiscriminantAnalysis(
+            solver="dummy",
+            covariance_estimator=LedoitWolf()
+    )
+    assert_raises(ValueError, clf.fit, X, y)
+
+    clf = QuadraticDiscriminantAnalysis(
+            solver="lsqr",
+            covariance_estimator=LedoitWolf()
+    )
+    with ignore_warnings():
+        clf.fit(X2, y6)
+    y_pred = clf.predict(X2)
+    assert_array_equal(y_pred, y6)
+
+    # Case n_samples_in_a_class < n_features
+    clf = QuadraticDiscriminantAnalysis(
+            solver="lsqr",
+            covariance_estimator=LedoitWolf()
+    )
+    with ignore_warnings():
+        clf.fit(X5, y5)
+    y_pred5 = clf.predict(X5)
+    assert_array_equal(y_pred5, y5)
+
+
+def test_qda_covariance_sample_estimate():
+
+    # Bug with lsqr solver if covariance is degenerate
+    clf = QuadraticDiscriminantAnalysis(
+            solver="lsqr",
+            covariance_estimator=EmpiricalCovariance()
+    )
+    clf.fit(X2, y6)
+    assert_raises(ValueError, clf.predict, X2)
+
+    X = np.random.rand(100, 10)
+    Y = np.random.randint(3, size=(100))
+    c1 = QuadraticDiscriminantAnalysis(
+            solver="lsqr",
+            covariance_estimator=EmpiricalCovariance()
+    )
+    c2 = QuadraticDiscriminantAnalysis(
+            solver="lsqr",
+            store_covariance=True
+    )
+    c1.fit(X, Y)
+    c2.fit(X, Y)
+    for i in range(2):
+        assert_array_almost_equal(c1.means_[i], c2.means_[i])
+    assert_array_almost_equal(c1.covariance_, c2.covariance_)
+
+
 def test_covariance():
     x, y = make_blobs(n_samples=100, n_features=5,
                       centers=1, random_state=42)


Code review comments:
diff_hunk: @@ -25,47 +25,63 @@
 from .utils.extmath import softmax
 from .preprocessing import StandardScaler
 
-
 __all__ = ['LinearDiscriminantAnalysis', 'QuadraticDiscriminantAnalysis']
 
 
-def _cov(X, shrinkage=None):
-    """Estimate covariance matrix (using optional shrinkage).
-
+def _cov(X, shrinkage=None, covariance_estimator=None, standardize=False):
+    """Estimate covariance matrix (using covariance_estimator).
     Parameters
     ----------
     X : array-like, shape (n_samples, n_features)
         Input data.
-
+    covariance_estimator : BaseEstimator
path: sklearn/discriminant_analysis.py
text: @user1:
I'm not sure what to put here, a `BaseEstimator` is definitely not enough. I'd say `estimator, or None, default=None`

diff_hunk: @@ -25,47 +25,63 @@
 from .utils.extmath import softmax
 from .preprocessing import StandardScaler
 
-
 __all__ = ['LinearDiscriminantAnalysis', 'QuadraticDiscriminantAnalysis']
 
 
-def _cov(X, shrinkage=None):
-    """Estimate covariance matrix (using optional shrinkage).
-
+def _cov(X, shrinkage=None, covariance_estimator=None, standardize=False):
+    """Estimate covariance matrix (using covariance_estimator).
     Parameters
     ----------
     X : array-like, shape (n_samples, n_features)
         Input data.
-
+    covariance_estimator : BaseEstimator
+        A covariance estimator, possible values
+            - None: shrinkage parameter drives the estimate
+            - BaseEstimator with a fit method and
+            a .covariance_ attribute like the estimators
+            in sklearn.covariance
+    standardize : bool
+        if True standardize the data before applying the covariance estimator
+        Note: only used when covariance_estimator is not None
     shrinkage : string or float, optional
         Shrinkage parameter, possible values:
           - None or 'empirical': no shrinkage (default).
           - 'auto': automatic shrinkage using the Ledoit-Wolf lemma.
           - float between 0 and 1: fixed shrinkage parameter.
-
     Returns
     -------
     s : array, shape (n_features, n_features)
         Estimated covariance matrix.
     """
-    shrinkage = "empirical" if shrinkage is None else shrinkage
-    if isinstance(shrinkage, str):
-        if shrinkage == 'auto':
-            sc = StandardScaler()  # standardize features
+    if covariance_estimator is None:
+        shrinkage = "empirical" if shrinkage is None else shrinkage
+        if isinstance(shrinkage, str):
+            if shrinkage == 'auto':
+                sc = StandardScaler()  # standardize features
+                X = sc.fit_transform(X)
+                s = ledoit_wolf(X)[0]
+                # rescale
+                s = sc.scale_[:, np.newaxis] * s * sc.scale_[np.newaxis, :]
+            elif shrinkage == 'empirical':
+                s = empirical_covariance(X)
+            else:
+                raise ValueError('unknown shrinkage parameter')
+        elif isinstance(shrinkage, float) or isinstance(shrinkage, int):
+            if shrinkage < 0 or shrinkage > 1:
+                raise ValueError('shrinkage parameter must be between 0 and 1')
+            s = shrunk_covariance(empirical_covariance(X), shrinkage)
+        else:
+            raise TypeError('shrinkage must be of string or int type')
+    else:
+        if standardize:
+            sc = StandardScaler()
             X = sc.fit_transform(X)
-            s = ledoit_wolf(X)[0]
-            # rescale
+            covariance_estimator.fit(X)
+            s = covariance_estimator.covariance_
path: sklearn/discriminant_analysis.py
text: @user1:
It's probably a good idea to do some validation here and check if the fitted estimator actually has `covariance_` set. If not, raise the appropriate RuntimeError or something with a helpful message.