Repository: scikit-learn/scikit-learn
Pull request: ENH: Add regularization to the main NMF class
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at
the contribution guidelines: https://github.com/scikit-learn/scikit-learn/blob/master/CONTRIBUTING.md#pull-request-checklist
-->

#### Reference Issues/PRs
<!--
Example: Fixes #1234. See also #3456.
Please use keywords (e.g., Fixes) to create link to the issues or pull requests
you resolved, so that they will automatically be closed when your pull request
is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests
-->
Fixes #12520

#### What does this implement/fix? Explain your changes.
Introduces the `regularization` parameter to the main `NMF()` class instead of just the public `non_negative_factorization` public function

#### Any other comments?
Need guidance on how to add tests to the `test_nmf_regularization()`
Tagging @adrinjalali for PR review, based on the issue
<!--
Please be aware that we are a loose team of volunteers so patience is
necessary; assistance handling other issues is very welcome. We value
all user contributions, no matter how minor they are. If we are slow to
review, either the pull request needs some benchmarking, tinkering,
convincing, etc. or more likely the reviewers are simply busy. In either
case, we ask for your understanding during the review process.
For more information, see our FAQ on this topic:
http://scikit-learn.org/dev/faq.html#why-is-my-pull-request-not-getting-any-attention.

Thanks for contributing!
-->


Problem statement:
non_negative_factorization's 'regularize' parameter is lost in the OO interface
#### Description
The function 'non_negative_factorization' in decomposition/nmf.py has a nice parameter 'regularization' where the regularization can be specified towards one of the product factors. There is however no way to access this feature when using an NMF object. Could be solved easily by carrying over the parameter to the constructor of the NMF object and creating a field 'self.regularize'.

Hints from issue discussion:
Can I pick this up?

From what  I see in the code it seems like a reasonable fix; but I may be wrong. @konqr you can give it a shot, and I guess we'll see in the PR if it's reasonable to use those different options. Please bare in mind that you'll need to add some tests as well. It's interesting that currently only `regularization='both'` is covered in our tests of the function.
Thanks @adrinjalali. I'll give it a shot.

Commit message: ENH: Add regularization to the main NMF class

Patch to review:
diff --git a/doc/modules/decomposition.rst b/doc/modules/decomposition.rst
index def755f42b796..330dc59c759ed 100644
--- a/doc/modules/decomposition.rst
+++ b/doc/modules/decomposition.rst
@@ -752,9 +752,9 @@ and the regularized objective function is:
     + \frac{\alpha(1-\rho)}{2} ||W||_{\mathrm{Fro}} ^ 2
     + \frac{\alpha(1-\rho)}{2} ||H||_{\mathrm{Fro}} ^ 2
 
-:class:`NMF` regularizes both W and H. The public function
-:func:`non_negative_factorization` allows a finer control through the
-:attr:`regularization` attribute, and may regularize only W, only H, or both.
+:class:`NMF` regularizes both W and H by default. The :attr:`regularization`
+parameter allows for finer control, with which only W, only H,
+or both can be regularized.
 
 NMF with a beta-divergence
 --------------------------
diff --git a/sklearn/decomposition/_nmf.py b/sklearn/decomposition/_nmf.py
index f1385d21596e3..0c6846bb42b74 100644
--- a/sklearn/decomposition/_nmf.py
+++ b/sklearn/decomposition/_nmf.py
@@ -1163,6 +1163,12 @@ class NMF(TransformerMixin, BaseEstimator):
     max_iter : integer, default: 200
         Maximum number of iterations before timing out.
 
+    regularization : 'both' | 'components' | 'transformation' | None
+        Select whether the regularization affects the components (H), the
+        transformation (W), both or none of them. Defaults to 'both'.
+
+        .. versionadded:: 0.24
+
     random_state : int, RandomState instance, default=None
         Used for initialisation (when ``init`` == 'nndsvdar' or
         'random'), and in Coordinate Descent. Pass an int for reproducible
@@ -1236,14 +1242,15 @@ class NMF(TransformerMixin, BaseEstimator):
     @_deprecate_positional_args
     def __init__(self, n_components=None, *, init=None, solver='cd',
                  beta_loss='frobenius', tol=1e-4, max_iter=200,
-                 random_state=None, alpha=0., l1_ratio=0., verbose=0,
-                 shuffle=False):
+                 regularization='both', random_state=None, alpha=0.,
+                 l1_ratio=0., verbose=0, shuffle=False):
         self.n_components = n_components
         self.init = init
         self.solver = solver
         self.beta_loss = beta_loss
         self.tol = tol
         self.max_iter = max_iter
+        self.regularization = regularization
         self.random_state = random_state
         self.alpha = alpha
         self.l1_ratio = l1_ratio
@@ -1283,7 +1290,7 @@ def fit_transform(self, X, y=None, W=None, H=None):
             X=X, W=W, H=H, n_components=self.n_components, init=self.init,
             update_H=True, solver=self.solver, beta_loss=self.beta_loss,
             tol=self.tol, max_iter=self.max_iter, alpha=self.alpha,
-            l1_ratio=self.l1_ratio, regularization='both',
+            l1_ratio=self.l1_ratio, regularization=self.regularization,
             random_state=self.random_state, verbose=self.verbose,
             shuffle=self.shuffle)
 
@@ -1332,9 +1339,10 @@ def transform(self, X):
             X=X, W=None, H=self.components_, n_components=self.n_components_,
             init=self.init, update_H=False, solver=self.solver,
             beta_loss=self.beta_loss, tol=self.tol, max_iter=self.max_iter,
-            alpha=self.alpha, l1_ratio=self.l1_ratio, regularization='both',
-            random_state=self.random_state, verbose=self.verbose,
-            shuffle=self.shuffle)
+            alpha=self.alpha, l1_ratio=self.l1_ratio,
+            regularization=self.regularization,
+            random_state=self.random_state,
+            verbose=self.verbose, shuffle=self.shuffle)
 
         return W
 
diff --git a/sklearn/decomposition/tests/test_nmf.py b/sklearn/decomposition/tests/test_nmf.py
index a7ef1243d8e25..50b0e180e36af 100644
--- a/sklearn/decomposition/tests/test_nmf.py
+++ b/sklearn/decomposition/tests/test_nmf.py
@@ -20,12 +20,14 @@
 
 
 @pytest.mark.parametrize('solver', ['cd', 'mu'])
-def test_convergence_warning(solver):
+@pytest.mark.parametrize('regularization',
+                         [None, 'both', 'components', 'transformation'])
+def test_convergence_warning(solver, regularization):
     convergence_warning = ("Maximum number of iterations 1 reached. "
                            "Increase it to improve convergence.")
     A = np.ones((2, 2))
     with pytest.warns(ConvergenceWarning, match=convergence_warning):
-        NMF(solver=solver, max_iter=1).fit(A)
+        NMF(solver=solver, regularization=regularization, max_iter=1).fit(A)
 
 
 def test_initialize_nn_output():
@@ -44,6 +46,8 @@ def test_parameter_checking():
     assert_raise_message(ValueError, msg, NMF(solver=name).fit, A)
     msg = "Invalid init parameter: got 'spam' instead of one of"
     assert_raise_message(ValueError, msg, NMF(init=name).fit, A)
+    msg = "Invalid regularization parameter: got 'spam' instead of one of"
+    assert_raise_message(ValueError, msg, NMF(regularization=name).fit, A)
     msg = "Invalid beta_loss parameter: got 'spam' instead of one"
     assert_raise_message(ValueError, msg, NMF(solver='mu',
                                               beta_loss=name).fit, A)
@@ -103,30 +107,36 @@ def test_nmf_fit_nn_output():
               5. + np.arange(1, 6)]
     for solver in ('cd', 'mu'):
         for init in (None, 'nndsvd', 'nndsvda', 'nndsvdar', 'random'):
-            model = NMF(n_components=2, solver=solver, init=init,
-                        random_state=0)
-            transf = model.fit_transform(A)
-            assert not((model.components_ < 0).any() or
-                       (transf < 0).any())
+            for regularization in (None, 'both',
+                                   'components', 'transformation'):
+                model = NMF(n_components=2, solver=solver, init=init,
+                            regularization=regularization, random_state=0)
+                transf = model.fit_transform(A)
+                assert not((model.components_ < 0).any() or
+                           (transf < 0).any())
 
 
 @pytest.mark.parametrize('solver', ('cd', 'mu'))
-def test_nmf_fit_close(solver):
+@pytest.mark.parametrize('regularization',
+                         (None, 'both', 'components', 'transformation'))
+def test_nmf_fit_close(solver, regularization):
     rng = np.random.mtrand.RandomState(42)
     # Test that the fit is not too far away
     pnmf = NMF(5, solver=solver, init='nndsvdar', random_state=0,
-               max_iter=600)
+               regularization=regularization, max_iter=600)
     X = np.abs(rng.randn(6, 5))
     assert pnmf.fit(X).reconstruction_err_ < 0.1
 
 
 @pytest.mark.parametrize('solver', ('cd', 'mu'))
-def test_nmf_transform(solver):
+@pytest.mark.parametrize('regularization',
+                         (None, 'both', 'components', 'transformation'))
+def test_nmf_transform(solver, regularization):
     # Test that NMF.transform returns close values
     rng = np.random.mtrand.RandomState(42)
     A = np.abs(rng.randn(6, 5))
     m = NMF(solver=solver, n_components=3, init='random',
-            random_state=0, tol=1e-5)
+            regularization=regularization, random_state=0, tol=1e-5)
     ft = m.fit_transform(A)
     t = m.transform(A)
     assert_array_almost_equal(ft, t, decimal=2)
@@ -148,12 +158,14 @@ def test_nmf_transform_custom_init():
 
 
 @pytest.mark.parametrize('solver', ('cd', 'mu'))
-def test_nmf_inverse_transform(solver):
+@pytest.mark.parametrize('regularization',
+                         (None, 'both', 'components', 'transformation'))
+def test_nmf_inverse_transform(solver, regularization):
     # Test that NMF.inverse_transform returns close values
     random_state = np.random.RandomState(0)
     A = np.abs(random_state.randn(6, 4))
     m = NMF(solver=solver, n_components=4, init='random', random_state=0,
-            max_iter=1000)
+            regularization=regularization, max_iter=1000)
     ft = m.fit_transform(A)
     A_new = m.inverse_transform(ft)
     assert_array_almost_equal(A, A_new, decimal=2)
@@ -176,9 +188,11 @@ def test_nmf_sparse_input():
     A_sparse = csc_matrix(A)
 
     for solver in ('cd', 'mu'):
-        est1 = NMF(solver=solver, n_components=5, init='random',
-                   random_state=0, tol=1e-2)
-        est2 = clone(est1)
+        for regularization in (None, 'both', 'components', 'transformation'):
+            est1 = NMF(solver=solver, n_components=5, init='random',
+                       regularization=regularization, random_state=0,
+                       tol=1e-2)
+            est2 = clone(est1)
 
     W1 = est1.fit_transform(A)
     W2 = est2.fit_transform(A_sparse)
@@ -213,16 +227,20 @@ def test_non_negative_factorization_consistency():
 
     for init in ['random', 'nndsvd']:
         for solver in ('cd', 'mu'):
-            W_nmf, H, _ = non_negative_factorization(
-                A, init=init, solver=solver, random_state=1, tol=1e-2)
-            W_nmf_2, _, _ = non_negative_factorization(
-                A, H=H, update_H=False, init=init, solver=solver,
-                random_state=1, tol=1e-2)
-
-            model_class = NMF(init=init, solver=solver, random_state=1,
-                              tol=1e-2)
-            W_cls = model_class.fit_transform(A)
-            W_cls_2 = model_class.transform(A)
+            for regularization in (None, 'both',
+                                   'components', 'transformation'):
+                W_nmf, H, _ = non_negative_factorization(
+                    A, init=init, solver=solver,
+                    regularization=regularization, random_state=1, tol=1e-2)
+                W_nmf_2, _, _ = non_negative_factorization(
+                    A, H=H, update_H=False, init=init, solver=solver,
+                    regularization=regularization, random_state=1, tol=1e-2)
+
+                model_class = NMF(init=init, solver=solver,
+                                  regularization=regularization,
+                                  random_state=1, tol=1e-2)
+                W_cls = model_class.fit_transform(A)
+                W_cls_2 = model_class.transform(A)
 
             assert_array_almost_equal(W_nmf, W_cls, decimal=10)
             assert_array_almost_equal(W_nmf_2, W_cls_2, decimal=10)
@@ -515,11 +533,13 @@ def test_nmf_underflow():
     (np.int32, np.float64),
     (np.int64, np.float64)])
 @pytest.mark.parametrize("solver", ["cd", "mu"])
-def test_nmf_dtype_match(dtype_in, dtype_out, solver):
+@pytest.mark.parametrize("regularization",
+                         (None, "both", "components", "transformation"))
+def test_nmf_dtype_match(dtype_in, dtype_out, solver, regularization):
     # Check that NMF preserves dtype (float32 and float64)
     X = np.random.RandomState(0).randn(20, 15).astype(dtype_in, copy=False)
     np.abs(X, out=X)
-    nmf = NMF(solver=solver)
+    nmf = NMF(solver=solver, regularization=regularization)
 
     assert nmf.fit(X).transform(X).dtype == dtype_out
     assert nmf.fit_transform(X).dtype == dtype_out
@@ -527,13 +547,15 @@ def test_nmf_dtype_match(dtype_in, dtype_out, solver):
 
 
 @pytest.mark.parametrize("solver", ["cd", "mu"])
-def test_nmf_float32_float64_consistency(solver):
+@pytest.mark.parametrize("regularization",
+                         (None, "both", "components", "transformation"))
+def test_nmf_float32_float64_consistency(solver, regularization):
     # Check that the result of NMF is the same between float32 and float64
     X = np.random.RandomState(0).randn(50, 7)
     np.abs(X, out=X)
-    nmf32 = NMF(solver=solver, random_state=0)
+    nmf32 = NMF(solver=solver, regularization=regularization, random_state=0)
     W32 = nmf32.fit_transform(X.astype(np.float32))
-    nmf64 = NMF(solver=solver, random_state=0)
+    nmf64 = NMF(solver=solver, regularization=regularization, random_state=0)
     W64 = nmf64.fit_transform(X)
 
     assert_allclose(W32, W64, rtol=1e-6, atol=1e-5)


Code review comments:
diff_hunk: @@ -1163,6 +1163,12 @@ class NMF(TransformerMixin, BaseEstimator):
     max_iter : integer, default: 200
         Maximum number of iterations before timing out.
 
+    regularization : 'both' | 'components' | 'transformation' | None
path: sklearn/decomposition/_nmf.py
text: @user1:
```suggestion
    regularization : {'both', 'components', 'transformation', None}, default=None
```

Should we use `'none'` instead of `None`?

@author:
This should be your call @user1 , but by convention in the rest of the codebase I still see `default=None` being used when the input is a string type, so I retained this 😅