What is plotted, is it a norm of the matrices? If so which?

While I agree that with your change to add data dependent scaling the transfer effect seems to be smaller than without the scaling, I am worried that this kind of scaling could be non-standard in the literature. Could you please try to review popular implementation of NMF to check if they do this too?

Have you seen this kind of scaling in the literature?

The plots show the mean of all elements in W and H, i.e. the L1 norm normalized by the number of elements.

It all came from @vene's comment (#4852 (comment)), and to avoid blocking the PR, I opened a separate PR to further talk about it.
I have never seen any such scaling in the literature, yet as @vene's puts it (#4852 (comment)), "it's a usability detail, [...] But usability is important". I guess this is the difference between a paper's implementation and a general public package's implementation. At the same time, I understand the need to avoid non-standard behaviors.

As #4852 is merged now, what's the status of this PR ?

Status: Controversial
We are not set yet on what we prefer:

• standard (current) behavior: no scaling of alpha, but the regularization depends on the size of the data, and the regularization is not balanced between W and H.
• non-standard (proposed) behavior: scaling alpha, so the regularization does not depend on the size of the data, and is balanced between W and H.

New idea: another slightly more standard method could be to normalize H (aka the dictionary) to unit norm at each iteration, yet I am not sure of the consequences on regularization.

Coming back to this PR after all those years. It's true that the proposed normalization for the regularizer makes a lot of sense from the look of the experiment. Could you try to do the same experiments on some natural datases (e.g. TF-IDF vectors or scientific signal data) to see if the normal of the matrices shrink fast close to alpha=1 as is the case on your synthetic experiments?

However we would need a new option FutureWarning to select between the 2 regularizer defintiions (in the public API only) along with a FutureWarning that explains that the default behavior will change.

I agree this looks very convincing especially if as Olivier asks we see the pattern in real data too. This seems more than a usability issue since we don't provide the API to search for separate alphas for W and H, which would be needed to construct an equivalent problem, right?

We could have a scale_alpha=False|True arg and change its default value.

So I tried with two datasets:

• Olivetti faces, (400, 4096), dense. I get the same results than with randomly generated data (not shown here).
• 20 newsgroup, (11314, 39115), sparse. I get pretty different results (see below). The new scaling definitely improves the situation, but not quite as neatly as in the dense case. We probably need to add some heuristic that takes into account the sparsity level.

without scaling (master branch)

with scaling (this branch)

import itertools
import warnings

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import NMF
from sklearn.exceptions import ConvergenceWarning

warnings.simplefilter('ignore', DeprecationWarning)
warnings.simplefilter('ignore', ConvergenceWarning)


def load_20news():
    from sklearn.datasets import fetch_20newsgroups
    from sklearn.feature_extraction.text import TfidfVectorizer
    dataset = fetch_20newsgroups(shuffle=True, random_state=1,
                                 remove=('headers', 'footers', 'quotes'))
    vectorizer = TfidfVectorizer(max_df=0.95, min_df=2, stop_words='english')
    X = vectorizer.fit_transform(dataset.data)
    return X


def load_faces():
    from sklearn.datasets import fetch_olivetti_faces
    X = fetch_olivetti_faces(shuffle=True).data
    return X


def load_random(n_samples, n_features):
    rng = np.random.RandomState(0)
    X = np.abs(rng.randn(n_samples, n_features))
    return X


dataset = "random"
if dataset == "20news":
    X = load_20news()
    n_samples_list = [100, 1000, 10000]
    n_features_list = [100, 1000, 10000]
elif dataset == "faces":
    X = load_faces()
    n_samples_list = [4, 40, 400]
    n_features_list = [4, 40, 400]
elif dataset == "random":
    X = load_random(1000, 1000)
    n_samples_list = [10, 100, 1000]
    n_features_list = [10, 100, 1000]

alphas = np.logspace(-7, 3, num=100)

estimator = NMF(n_components=10, solver='cd', init='random',
                tol=1e-4, max_iter=10, random_state=0, alpha=0., l1_ratio=0.,
                verbose=0, shuffle=True)

fig, axes = plt.subplots(3, 3, figsize=(20, 15))
for (n_samples,
     n_features), ax in zip(itertools.product(n_samples_list, n_features_list),
                            axes.ravel()):
    print('%d - %d' % (n_samples, n_features))

    mean_w = np.zeros(len(alphas))
    mean_h = np.zeros(len(alphas))
    for ii, alpha in enumerate(alphas):
        # alpha *= np.sqrt(n_samples * n_features)
        estimator = estimator.set_params(alpha=alpha)
        W = estimator.fit_transform(X[:n_samples, :n_features])
        H = estimator.components_
        mean_w[ii] = W.mean()
        mean_h[ii] = H.mean()

    ax.set_xscale('log')
    ax.plot(alphas, mean_w, label='W')
    ax.plot(alphas, mean_h, label='H')
    ax.legend(loc=0)
    ax.set_title('X.shape = (%d, %d)' % (n_samples, n_features))

fig.suptitle('Mean coefficient, w.r.t alpha L2 regularisation', size=16)
fig.savefig(f'{dataset}_master.png')
plt.show()

Closed by #20512