This PR will need to be rebased on top of the current master and the input validation need to be updated to use the new utilities, see:

• http://scikit-learn.org/stable/developers/utilities.html

Also as @larsmans said in the previous PR the narrative doc in ./doc/modules/lda_qda.rst needs to be updated to highlight this new option.

Finally it would be great to write and example in the examples folder to plot the impact of using shrinkage to improve performance accuracy and use that plot in the narrative documentation.

Once all this change are made please change the title of this PR from [WIP] to [MRG] for final review.

@cle1109 I've updated examples/plot_lda.py from the old PR to use the new LDA API: 2ef847f.
You can pull it from the slda branch in my fork.

Changes Unknown when pulling 5a18f80 on cle1109:slda into * on scikit-learn:master*.

OK, I think this should be it, everything should be done.

Changes Unknown when pulling e610f87 on cle1109:slda into * on scikit-learn:master*.

Changes Unknown when pulling 64dc88b on cle1109:slda into * on scikit-learn:master*.

Changes Unknown when pulling 807c83f on cle1109:slda into * on scikit-learn:master*.

OK, so here's my take on the discussions in #3500 and #3105 (might I suggest to move it here, because here's where our latest code is).

I really like the idea of introducing solver and alpha parameters. The default value for solver could be SVD (the current implementation with the bug fixed in the transform method). The other option could be covariance. The alpha parameter should work in both cases, but the value ledoit_wolf makes only sense in the latter.

I support the idea of introducing the possibility to optimize the parameter via cross-validation.

The other question is how close to the textbook we want to be in the covariance method. As @kazemakase pointed out, we don't need to solve the eigenvalue problem for classification. That's also how most algorithms are implemented, because inverting the between scatter matrix could be a problem (see also the discussion at Cross Validated).

So in principle we could go ahead and start implementing it according to all suggestions here. The only thing that needs to be clarified first is how we do both classification and transform in the covariance method.

You need some sort of eigen/singular value decomposition for the transform, but you can fit the classifier with a matrix inverse or an equation solver.

Instead of introducing a fit_transform option as I previously suggested, I think we could introduce two solvers covariance-lsqr and covariance-eigen. The transform method would only be available if covariance-eigen was used.

We could also add a covariance-cg solver which uses conjugate gradient instead of lsqr for solving the least squares problem / system of equations. This should make it possible to solve the problem without materializing the covariance matrix in memory (see the _solve_sparse_cg function in the Ridge module). This can of course be addressed in another PR.

(I took the liberty to change the title of this PR and add a TODO list.)

Actually, we can drop the covariance- prefix. This way, the names will be consistent with the Ridge module.

Just svd. The fact that lsqr and eigen use the covariance matrix can be documented in the docstring. In the Ridge module, we give a short description of each solver in the docstring:
https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/linear_model/ridge.py#L215.

OK! (I deleted my comment because I saw it in the todo list at the top)

i got a comment on one of the items in the todo list:

test that solvers return the same results (up to numerical errors)

Numerical errors can blow up and cause big differences in classification results if feature space is high dimensional. That is because these errors cause a small amount of uncertainty in the orientation of the hyperplane. The volume of this uncertainty increases exponentially with dimensionality. The bigger this volume the higher the chance of unseen samples lying on opposing sides of the hyperplane in different implementations.

In short, classification results can differ noticeably due to numerical errors. Something to keep in mind when designing the tests.

@agramfort, I am very confident now that the code does what it is supposed to do :-). The minor issue of differently scaled scalings_ should be addressed in a follow-up PR (it doesn't affect the functionality, but for consistency it would be nice if the same coefficients are obtained by both solvers. What do you think (in particular @ogrisel)?

I think we should merge this. Maybe squash the commits, there are quite a few.

Do you want me to squash the commits (to only one?) - or will you do that?

@amueller?

I don't mind pushing these commits as they are... let's merge this ASAP.

I don't mind pushing these commits as they are... let's merge this ASAP.

+1

Ok, it looks like we don't do squashing any more... I felt that was helpful for bugfix releases but ok....

I felt that was helpful for bugfix releases but ok....

I agree with you. I think that we should encourage it. But it don't feel
that it is mandatory.

🍻 @cle1109 !

😄 - awesome! Thank you guys, I really learned a lot while working on this PR. Looking forward to contributing more soon!

This fixes #1649, right? Or only for the new solvers?

I think so

@agramfort is that to "fixed #1649" or to "only for the new solvers"?