This issue discovered in #16948 has 2 parts.

• The default init for NMF is None, which corresponds to "nndsvd". However there's this warning

if solver == 'mu' and init == 'nndsvd':
        warnings.warn("The multiplicative update ('mu') solver cannot update "
                      "zeros present in the initialization, and so leads to "
                      "poorer results when used jointly with init='nndsvd'. "
                      "You may try init='nndsvda' or init='nndsvdar' instead.",
                      UserWarning)

which is not issued when init is left to None. If it's discouraged to use nndsvd with the multiplicative update, don't you think we should change the default to be nndsvda when solver == 'mu' ? We could keep None or name it "auto".

• There's a discrepancy between fit.transform and fit_transform, especially with 'mu' solver.
  The current default for the init being nndsvd, W likely has zeros in its initialization which are propagated all the way using fit_transform. Using fit.transform is different because for transform, W is initialized as a constant matrix.
  Maybe changing the default init to nndsvda can be considered enough. Otherwise I wonder if transform should initialize W using the same function as fit. wdyt ?
  edit: changing the default init still fails the common test fit.transform vs fit_transform for the multiplicative update solver. Using the same init at transform also does not work.

ping @cmarmo @TomDLT