[MRG+1] Add multiplicative-update solver in NMF, with all beta-divergence

[TomDLT]

This PR is a second part of what have been discussed in #4811. (first part, #4852, merged)
It includes :

• a Multiplicative-Update solver in NMF. This solver is generally slower than the Coordinate-Descent solver, but it handles all beta-divergences (including Frobenius, generalized Kullback_Leibler and Itakura-Saito).
• a plot to visualize the beta-divergence for several values of beta.
• benchmarks below

Link to the rendered doc

---

plot_beta_divergence.py

---

This change is

[mblondel]

I would favor using CD for beta divergences as well. We decided to remove the slow pg solver. I would prefer not add a new slow solver if possible. The KDD paper shows that for generalized KL, we can just use the element-wise Newton update without line search. For other divergences, this will need some verification (unless there are more recent works that cover this?) but in worst case we can implement a simple line search. As a start, it would be nice to implement the CD update for generalized KL to compare.

[tttthomasssss]

Could anybody comment on when this is planned to be merged into master? I'm quite keen on using NMF with generalized KL divergence.

[agramfort]

@mblondel maybe we can do better than multiplicative updates but I think we should be pragmatic here. We can deprecate after when we have something better. wdyt?

[TomDLT]

Sorry for the long silence, I have not as much time for this as before.
@mblondel I spent some time trying to implement the coordinate descent for generalized Kullback Leibler divergence, as described in Coordinate descent with greedy selection, at least in the dense case.
However, I did not manage to be faster than multiplicative update, which is not what is indicated in the paper, and probably underlines my own limits (or proves that my implementation of multiplicative update is very good :)). More seriously, even if we manage to be faster, the extension to all beta-divergence is not trivial and would require much more work.

I understand that multiplicative update is an old method and probably not the fastest, yet it is still a good way to improve the NMF in scikit-learn, extending to a lot of different loss functions.

[TomDLT]

I compared my implementation of multiplicative update with #2540 (Frobenius norm) and #1348 (KL divergence).

---

Comparing to #2540, the results are identical, provided some slight modifications:
To avoid division by zero, what do we prefer?:

1. A / (B + eps) (like in [WIP] NMF estimator based on the Lee and Seung algorithm #2540)
2. (A + eps) / (B + eps) (like in NMF for Kullback-Leibler divergence #1348)
3. B[B == 0] = eps then A / B
   I have a preference for the last one, for it does not affect update with nonzero elements. WDYT ?

---

Comparing to #1348, the results are different since it forgets a division in the update
(see Eq 5 here), and adds a normalization of H.

[TomDLT]

Some benchmarking:
20_news - sparse (11314, 39116)

Faces - dense (400, 4096)

MNIST - dense (70000, 784)

RCV1 - sparse (804414, 47236)

The MU solver is not very fast, but it was implemented not for its speed, but for it handles all beta-divergences. Note that the poor performances with 'nndsvd' initialization (top-left) are expected, since this initialization have a lot of zeros that cannot be modified by a Multiplicative Update.

[tguillemot]

@TomDLT @mblondel @agramfort This PR is mergeable ?
I will do a first round of review but I don't know those methods. Someone can have a look with me ?

[tguillemot]

We can use pandas for benchs ?

[TomDLT]

Using it in a benchmark does not make it a dependency of the package, as for matplotlib.

[agramfort]

@mblondel what's your take on this PR?

[mblondel]

I spent some time trying to implement the coordinate descent for generalized Kullback Leibler divergence, as described in Coordinate descent with greedy selection, at least in the dense case.

Thanks for the investigation.

More seriously, even if we manage to be faster, the extension to all beta-divergence is not trivial and would require much more work.

Agreed.

@mblondel what's your take on this PR?

+1 on my side and sorry for the late reply. As a big fan of CD, I would be potentially interested in the CD code for GKL divergence as it could be faster in the sparse case.

[agramfort]

cool. Looks like we can proceed then :) let's make it a WIP->MRG then

[TomDLT]

let's make it a WIP->MRG then

it is, all reviews will be greatly appreciated.

[tguillemot]

I would change a bit the formulation :
"Other distance functions can be used in NMF as, for example, the (generalized) Kullback-Leibler (KL) divergence, also referred as I-divergence:

..math::
d_{KL}(X, Y) = \sum_{i,j} (X_{ij} * log(\frac{X_{ij}}{Y_{ij}}) - X_{ij} + Y_{ij})

Or, the the Itakura-Saito (IS) divergence:"
See if you prefer.

[TomDLT]

Done

[amueller]

Thanks. Hm so @mblondel I saw a +1 from you up there, but that wasn't a full review, was it? @agramfort did you review? I see +1 from @ogrisel

[agramfort]

no I did not. No time :(

[VictorBst]

Hi, @TomDLT and @agramfort invited to have a look at the code.

I mostly looked at the math aspect of the code and it seems to be correct. After quick experiments on my data, everything seems to work as intended and gives similar results to other simple NMF python codes I tested. As a side note, this one was always either equal or the fastest in my tests, especially for higher dimensions.

[amueller]

hm I'm not sure if this should be merged before #7927, there are some conflicts, I think.... opinions @TomDLT ?

[TomDLT]

ok let's wait for #7927

[ogrisel]

@TomDLT now that #7927 has been merged, this needs a rebase.

[TomDLT]

Done

[ogrisel]

Thanks @TomDLT (and @VictorBst for the review). I squash-merged. 🍻

[TomDLT]

🍾 Cool ! Thanks a lot for the reviews ! 🍾

[tguillemot]

Yeah !!!

[amueller]

Awesome! Congrats!