Aluode/PerceptionLabPortable
0
1"""Non-negative matrix factorization."""
2
3# Authors: The scikit-learn developers
4# SPDX-License-Identifier: BSD-3-Clause
5
6import itertools
7import time
8import warnings
9from abc import ABC
10from math import sqrt
11from numbers import Integral, Real
12
13import numpy as np
14import scipy.sparse as sp
15from scipy import linalg
16
17from .._config import config_context
18from ..base import (
19 BaseEstimator,
20 ClassNamePrefixFeaturesOutMixin,
21 TransformerMixin,
22 _fit_context,
23)
24from ..exceptions import ConvergenceWarning
25from ..utils import check_array, check_random_state, gen_batches
26from ..utils._param_validation import (
27 Interval,
28 StrOptions,
29 validate_params,
30)
31from ..utils.extmath import _randomized_svd, safe_sparse_dot, squared_norm
32from ..utils.validation import (
33 check_is_fitted,
34 check_non_negative,
35 validate_data,
36)
37from ._cdnmf_fast import _update_cdnmf_fast
38
39EPSILON = np.finfo(np.float32).eps
40
41
42def norm(x):
43 """Dot product-based Euclidean norm implementation.
44
45 See: http://fa.bianp.net/blog/2011/computing-the-vector-norm/
46
47 Parameters
48 ----------
49 x : array-like
50 Vector for which to compute the norm.
51 """
52 return sqrt(squared_norm(x))
53
54
55def trace_dot(X, Y):
56 """Trace of np.dot(X, Y.T).
57
58 Parameters
59 ----------
60 X : array-like
61 First matrix.
62 Y : array-like
63 Second matrix.
64 """
65 return np.dot(X.ravel(), Y.ravel())
66
67
68def _check_init(A, shape, whom):
69 A = check_array(A)
70 if shape[0] != "auto" and A.shape[0] != shape[0]:
71 raise ValueError(
72 f"Array with wrong first dimension passed to {whom}. Expected {shape[0]}, "
73 f"but got {A.shape[0]}."
74 )
75 if shape[1] != "auto" and A.shape[1] != shape[1]:
76 raise ValueError(
77 f"Array with wrong second dimension passed to {whom}. Expected {shape[1]}, "
78 f"but got {A.shape[1]}."
79 )
80 check_non_negative(A, whom)
81 if np.max(A) == 0:
82 raise ValueError(f"Array passed to {whom} is full of zeros.")
83
84
85def _beta_divergence(X, W, H, beta, square_root=False):
86 """Compute the beta-divergence of X and dot(W, H).
87
88 Parameters
89 ----------
90 X : float or array-like of shape (n_samples, n_features)
91
92 W : float or array-like of shape (n_samples, n_components)
93
94 H : float or array-like of shape (n_components, n_features)
95
96 beta : float or {'frobenius', 'kullback-leibler', 'itakura-saito'}
97 Parameter of the beta-divergence.
98 If beta == 2, this is half the Frobenius *squared* norm.
99 If beta == 1, this is the generalized Kullback-Leibler divergence.
100 If beta == 0, this is the Itakura-Saito divergence.
101 Else, this is the general beta-divergence.
102
103 square_root : bool, default=False
104 If True, return np.sqrt(2 * res)
105 For beta == 2, it corresponds to the Frobenius norm.
106
107 Returns
108 -------
109 res : float
110 Beta divergence of X and np.dot(X, H).
111 """
112 beta = _beta_loss_to_float(beta)
113
114 # The method can be called with scalars
115 if not sp.issparse(X):
116 X = np.atleast_2d(X)
117 W = np.atleast_2d(W)
118 H = np.atleast_2d(H)
119
120 # Frobenius norm
121 if beta == 2:
122 # Avoid the creation of the dense np.dot(W, H) if X is sparse.
123 if sp.issparse(X):
124 norm_X = np.dot(X.data, X.data)
125 norm_WH = trace_dot(np.linalg.multi_dot([W.T, W, H]), H)
126 cross_prod = trace_dot((X @ H.T), W)
127 res = (norm_X + norm_WH - 2.0 * cross_prod) / 2.0
128 else:
129 res = squared_norm(X - np.dot(W, H)) / 2.0
130
131 if square_root:
132 return np.sqrt(res * 2)
133 else:
134 return res
135
136 if sp.issparse(X):
137 # compute np.dot(W, H) only where X is nonzero
138 WH_data = _special_sparse_dot(W, H, X).data
139 X_data = X.data
140 else:
141 WH = np.dot(W, H)
142 WH_data = WH.ravel()
143 X_data = X.ravel()
144
145 # do not affect the zeros: here 0 ** (-1) = 0 and not infinity
146 indices = X_data > EPSILON
147 WH_data = WH_data[indices]
148 X_data = X_data[indices]
149
150 # used to avoid division by zero
151 WH_data[WH_data < EPSILON] = EPSILON
152
153 # generalized Kullback-Leibler divergence
154 if beta == 1:
155 # fast and memory efficient computation of np.sum(np.dot(W, H))
156 sum_WH = np.dot(np.sum(W, axis=0), np.sum(H, axis=1))
157 # computes np.sum(X * log(X / WH)) only where X is nonzero
158 div = X_data / WH_data
159 res = np.dot(X_data, np.log(div))
160 # add full np.sum(np.dot(W, H)) - np.sum(X)
161 res += sum_WH - X_data.sum()
162
163 # Itakura-Saito divergence
164 elif beta == 0:
165 div = X_data / WH_data
166 res = np.sum(div) - np.prod(X.shape) - np.sum(np.log(div))
167
168 # beta-divergence, beta not in (0, 1, 2)
169 else:
170 if sp.issparse(X):
171 # slow loop, but memory efficient computation of :
172 # np.sum(np.dot(W, H) ** beta)
173 sum_WH_beta = 0
174 for i in range(X.shape[1]):
175 sum_WH_beta += np.sum(np.dot(W, H[:, i]) ** beta)
176
177 else:
178 sum_WH_beta = np.sum(WH**beta)
179
180 sum_X_WH = np.dot(X_data, WH_data ** (beta - 1))
181 res = (X_data**beta).sum() - beta * sum_X_WH
182 res += sum_WH_beta * (beta - 1)
183 res /= beta * (beta - 1)
184
185 if square_root:
186 res = max(res, 0) # avoid negative number due to rounding errors
187 return np.sqrt(2 * res)
188 else:
189 return res
190
191
192def _special_sparse_dot(W, H, X):
193 """Computes np.dot(W, H), only where X is non zero."""
194 if sp.issparse(X):
195 ii, jj = X.nonzero()
196 n_vals = ii.shape[0]
197 dot_vals = np.empty(n_vals)
198 n_components = W.shape[1]
199
200 batch_size = max(n_components, n_vals // n_components)
201 for start in range(0, n_vals, batch_size):
202 batch = slice(start, start + batch_size)
203 dot_vals[batch] = np.multiply(W[ii[batch], :], H.T[jj[batch], :]).sum(
204 axis=1
205 )
206
207 WH = sp.coo_matrix((dot_vals, (ii, jj)), shape=X.shape)
208 return WH.tocsr()
209 else:
210 return np.dot(W, H)
211
212
213def _beta_loss_to_float(beta_loss):
214 """Convert string beta_loss to float."""
215 beta_loss_map = {"frobenius": 2, "kullback-leibler": 1, "itakura-saito": 0}
216 if isinstance(beta_loss, str):
217 beta_loss = beta_loss_map[beta_loss]
218 return beta_loss
219
220
221def _initialize_nmf(X, n_components, init=None, eps=1e-6, random_state=None):
222 """Algorithms for NMF initialization.
223
224 Computes an initial guess for the non-negative
225 rank k matrix approximation for X: X = WH.
226
227 Parameters
228 ----------
229 X : array-like of shape (n_samples, n_features)
230 The data matrix to be decomposed.
231
232 n_components : int
233 The number of components desired in the approximation.
234
235 init : {'random', 'nndsvd', 'nndsvda', 'nndsvdar'}, default=None
236 Method used to initialize the procedure.
237 Valid options:
238
239 - None: 'nndsvda' if n_components <= min(n_samples, n_features),
240 otherwise 'random'.
241
242 - 'random': non-negative random matrices, scaled with:
243 sqrt(X.mean() / n_components)
244
245 - 'nndsvd': Nonnegative Double Singular Value Decomposition (NNDSVD)
246 initialization (better for sparseness)
247
248 - 'nndsvda': NNDSVD with zeros filled with the average of X
249 (better when sparsity is not desired)
250
251 - 'nndsvdar': NNDSVD with zeros filled with small random values
252 (generally faster, less accurate alternative to NNDSVDa
253 for when sparsity is not desired)
254
255 - 'custom': use custom matrices W and H
256
257 .. versionchanged:: 1.1
258 When `init=None` and n_components is less than n_samples and n_features
259 defaults to `nndsvda` instead of `nndsvd`.
260
261 eps : float, default=1e-6
262 Truncate all values less then this in output to zero.
263
264 random_state : int, RandomState instance or None, default=None
265 Used when ``init`` == 'nndsvdar' or 'random'. Pass an int for
266 reproducible results across multiple function calls.
267 See :term:`Glossary <random_state>`.
268
269 Returns
270 -------
271 W : array-like of shape (n_samples, n_components)
272 Initial guesses for solving X ~= WH.
273
274 H : array-like of shape (n_components, n_features)
275 Initial guesses for solving X ~= WH.
276
277 References
278 ----------
279 C. Boutsidis, E. Gallopoulos: SVD based initialization: A head start for
280 nonnegative matrix factorization - Pattern Recognition, 2008
281 http://tinyurl.com/nndsvd
282 """
283 check_non_negative(X, "NMF initialization")
284 n_samples, n_features = X.shape
285
286 if (
287 init is not None
288 and init != "random"
289 and n_components > min(n_samples, n_features)
290 ):
291 raise ValueError(
292 "init = '{}' can only be used when "
293 "n_components <= min(n_samples, n_features)".format(init)
294 )
295
296 if init is None:
297 if n_components <= min(n_samples, n_features):
298 init = "nndsvda"
299 else:
300 init = "random"
301
302 # Random initialization
303 if init == "random":
304 avg = np.sqrt(X.mean() / n_components)
305 rng = check_random_state(random_state)
306 H = avg * rng.standard_normal(size=(n_components, n_features)).astype(
307 X.dtype, copy=False
308 )
309 W = avg * rng.standard_normal(size=(n_samples, n_components)).astype(
310 X.dtype, copy=False
311 )
312 np.abs(H, out=H)
313 np.abs(W, out=W)
314 return W, H
315
316 # NNDSVD initialization
317 U, S, V = _randomized_svd(X, n_components, random_state=random_state)
318 W = np.zeros_like(U)
319 H = np.zeros_like(V)
320
321 # The leading singular triplet is non-negative
322 # so it can be used as is for initialization.
323 W[:, 0] = np.sqrt(S[0]) * np.abs(U[:, 0])
324 H[0, :] = np.sqrt(S[0]) * np.abs(V[0, :])
325
326 for j in range(1, n_components):
327 x, y = U[:, j], V[j, :]
328
329 # extract positive and negative parts of column vectors
330 x_p, y_p = np.maximum(x, 0), np.maximum(y, 0)
331 x_n, y_n = np.abs(np.minimum(x, 0)), np.abs(np.minimum(y, 0))
332
333 # and their norms
334 x_p_nrm, y_p_nrm = norm(x_p), norm(y_p)
335 x_n_nrm, y_n_nrm = norm(x_n), norm(y_n)
336
337 m_p, m_n = x_p_nrm * y_p_nrm, x_n_nrm * y_n_nrm
338
339 # choose update
340 if m_p > m_n:
341 u = x_p / x_p_nrm
342 v = y_p / y_p_nrm
343 sigma = m_p
344 else:
345 u = x_n / x_n_nrm
346 v = y_n / y_n_nrm
347 sigma = m_n
348
349 lbd = np.sqrt(S[j] * sigma)
350 W[:, j] = lbd * u
351 H[j, :] = lbd * v
352
353 W[W < eps] = 0
354 H[H < eps] = 0
355
356 if init == "nndsvd":
357 pass
358 elif init == "nndsvda":
359 avg = X.mean()
360 W[W == 0] = avg
361 H[H == 0] = avg
362 elif init == "nndsvdar":
363 rng = check_random_state(random_state)
364 avg = X.mean()
365 W[W == 0] = abs(avg * rng.standard_normal(size=len(W[W == 0])) / 100)
366 H[H == 0] = abs(avg * rng.standard_normal(size=len(H[H == 0])) / 100)
367 else:
368 raise ValueError(
369 "Invalid init parameter: got %r instead of one of %r"
370 % (init, (None, "random", "nndsvd", "nndsvda", "nndsvdar"))
371 )
372
373 return W, H
374
375
376def _update_coordinate_descent(X, W, Ht, l1_reg, l2_reg, shuffle, random_state):
377 """Helper function for _fit_coordinate_descent.
378
379 Update W to minimize the objective function, iterating once over all
380 coordinates. By symmetry, to update H, one can call
381 _update_coordinate_descent(X.T, Ht, W, ...).
382
383 """
384 n_components = Ht.shape[1]
385
386 HHt = np.dot(Ht.T, Ht)
387 XHt = safe_sparse_dot(X, Ht)
388
389 # L2 regularization corresponds to increase of the diagonal of HHt
390 if l2_reg != 0.0:
391 # adds l2_reg only on the diagonal
392 HHt.flat[:: n_components + 1] += l2_reg
393 # L1 regularization corresponds to decrease of each element of XHt
394 if l1_reg != 0.0:
395 XHt -= l1_reg
396
397 if shuffle:
398 permutation = random_state.permutation(n_components)
399 else:
400 permutation = np.arange(n_components)
401 # The following seems to be required on 64-bit Windows w/ Python 3.5.
402 permutation = np.asarray(permutation, dtype=np.intp)
403 return _update_cdnmf_fast(W, HHt, XHt, permutation)
404
405
406def _fit_coordinate_descent(
407 X,
408 W,
409 H,
410 tol=1e-4,
411 max_iter=200,
412 l1_reg_W=0,
413 l1_reg_H=0,
414 l2_reg_W=0,
415 l2_reg_H=0,
416 update_H=True,
417 verbose=0,
418 shuffle=False,
419 random_state=None,
420):
421 """Compute Non-negative Matrix Factorization (NMF) with Coordinate Descent
422
423 The objective function is minimized with an alternating minimization of W
424 and H. Each minimization is done with a cyclic (up to a permutation of the
425 features) Coordinate Descent.
426
427 Parameters
428 ----------
429 X : array-like of shape (n_samples, n_features)
430 Constant matrix.
431
432 W : array-like of shape (n_samples, n_components)
433 Initial guess for the solution.
434
435 H : array-like of shape (n_components, n_features)
436 Initial guess for the solution.
437
438 tol : float, default=1e-4
439 Tolerance of the stopping condition.
440
441 max_iter : int, default=200
442 Maximum number of iterations before timing out.
443
444 l1_reg_W : float, default=0.
445 L1 regularization parameter for W.
446
447 l1_reg_H : float, default=0.
448 L1 regularization parameter for H.
449
450 l2_reg_W : float, default=0.
451 L2 regularization parameter for W.
452
453 l2_reg_H : float, default=0.
454 L2 regularization parameter for H.
455
456 update_H : bool, default=True
457 Set to True, both W and H will be estimated from initial guesses.
458 Set to False, only W will be estimated.
459
460 verbose : int, default=0
461 The verbosity level.
462
463 shuffle : bool, default=False
464 If true, randomize the order of coordinates in the CD solver.
465
466 random_state : int, RandomState instance or None, default=None
467 Used to randomize the coordinates in the CD solver, when
468 ``shuffle`` is set to ``True``. Pass an int for reproducible
469 results across multiple function calls.
470 See :term:`Glossary <random_state>`.
471
472 Returns
473 -------
474 W : ndarray of shape (n_samples, n_components)
475 Solution to the non-negative least squares problem.
476
477 H : ndarray of shape (n_components, n_features)
478 Solution to the non-negative least squares problem.
479
480 n_iter : int
481 The number of iterations done by the algorithm.
482
483 References
484 ----------
485 .. [1] :doi:`"Fast local algorithms for large scale nonnegative matrix and tensor
486 factorizations" <10.1587/transfun.E92.A.708>`
487 Cichocki, Andrzej, and P. H. A. N. Anh-Huy. IEICE transactions on fundamentals
488 of electronics, communications and computer sciences 92.3: 708-721, 2009.
489 """
490 # so W and Ht are both in C order in memory
491 Ht = check_array(H.T, order="C")
492 X = check_array(X, accept_sparse="csr")
493
494 rng = check_random_state(random_state)
495
496 for n_iter in range(1, max_iter + 1):
497 violation = 0.0
498
499 # Update W
500 violation += _update_coordinate_descent(
501 X, W, Ht, l1_reg_W, l2_reg_W, shuffle, rng
502 )
503 # Update H
504 if update_H:
505 violation += _update_coordinate_descent(
506 X.T, Ht, W, l1_reg_H, l2_reg_H, shuffle, rng
507 )
508
509 if n_iter == 1:
510 violation_init = violation
511
512 if violation_init == 0:
513 break
514
515 if verbose:
516 print("violation:", violation / violation_init)
517
518 if violation / violation_init <= tol:
519 if verbose:
520 print("Converged at iteration", n_iter + 1)
521 break
522
523 return W, Ht.T, n_iter
524
525
526def _multiplicative_update_w(
527 X,
528 W,
529 H,
530 beta_loss,
531 l1_reg_W,
532 l2_reg_W,
533 gamma,
534 H_sum=None,
535 HHt=None,
536 XHt=None,
537 update_H=True,
538):
539 """Update W in Multiplicative Update NMF."""
540 if beta_loss == 2:
541 # Numerator
542 if XHt is None:
543 XHt = safe_sparse_dot(X, H.T)
544 if update_H:
545 # avoid a copy of XHt, which will be re-computed (update_H=True)
546 numerator = XHt
547 else:
548 # preserve the XHt, which is not re-computed (update_H=False)
549 numerator = XHt.copy()
550
551 # Denominator
552 if HHt is None:
553 HHt = np.dot(H, H.T)
554 denominator = np.dot(W, HHt)
555
556 else:
557 # Numerator
558 # if X is sparse, compute WH only where X is non zero
559 WH_safe_X = _special_sparse_dot(W, H, X)
560 if sp.issparse(X):
561 WH_safe_X_data = WH_safe_X.data
562 X_data = X.data
563 else:
564 WH_safe_X_data = WH_safe_X
565 X_data = X
566 # copy used in the Denominator
567 WH = WH_safe_X.copy()
568 if beta_loss - 1.0 < 0:
569 WH[WH < EPSILON] = EPSILON
570
571 # to avoid taking a negative power of zero
572 if beta_loss - 2.0 < 0:
573 WH_safe_X_data[WH_safe_X_data < EPSILON] = EPSILON
574
575 if beta_loss == 1:
576 np.divide(X_data, WH_safe_X_data, out=WH_safe_X_data)
577 elif beta_loss == 0:
578 # speeds up computation time
579 # refer to /numpy/numpy/issues/9363
580 WH_safe_X_data **= -1
581 WH_safe_X_data **= 2
582 # element-wise multiplication
583 WH_safe_X_data *= X_data
584 else:
585 WH_safe_X_data **= beta_loss - 2
586 # element-wise multiplication
587 WH_safe_X_data *= X_data
588
589 # here numerator = dot(X * (dot(W, H) ** (beta_loss - 2)), H.T)
590 numerator = safe_sparse_dot(WH_safe_X, H.T)
591
592 # Denominator
593 if beta_loss == 1:
594 if H_sum is None:
595 H_sum = np.sum(H, axis=1) # shape(n_components, )
596 denominator = H_sum[np.newaxis, :]
597
598 else:
599 # computation of WHHt = dot(dot(W, H) ** beta_loss - 1, H.T)
600 if sp.issparse(X):
601 # memory efficient computation
602 # (compute row by row, avoiding the dense matrix WH)
603 WHHt = np.empty(W.shape)
604 for i in range(X.shape[0]):
605 WHi = np.dot(W[i, :], H)
606 if beta_loss - 1 < 0:
607 WHi[WHi < EPSILON] = EPSILON
608 WHi **= beta_loss - 1
609 WHHt[i, :] = np.dot(WHi, H.T)
610 else:
611 WH **= beta_loss - 1
612 WHHt = np.dot(WH, H.T)
613 denominator = WHHt
614
615 # Add L1 and L2 regularization
616 if l1_reg_W > 0:
617 denominator += l1_reg_W
618 if l2_reg_W > 0:
619 denominator = denominator + l2_reg_W * W
620 denominator[denominator == 0] = EPSILON
621
622 numerator /= denominator
623 delta_W = numerator
624
625 # gamma is in ]0, 1]
626 if gamma != 1:
627 delta_W **= gamma
628
629 W *= delta_W
630
631 return W, H_sum, HHt, XHt
632
633
634def _multiplicative_update_h(
635 X, W, H, beta_loss, l1_reg_H, l2_reg_H, gamma, A=None, B=None, rho=None
636):
637 """update H in Multiplicative Update NMF."""
638 if beta_loss == 2:
639 numerator = safe_sparse_dot(W.T, X)
640 denominator = np.linalg.multi_dot([W.T, W, H])
641
642 else:
643 # Numerator
644 WH_safe_X = _special_sparse_dot(W, H, X)
645 if sp.issparse(X):
646 WH_safe_X_data = WH_safe_X.data
647 X_data = X.data
648 else:
649 WH_safe_X_data = WH_safe_X
650 X_data = X
651 # copy used in the Denominator
652 WH = WH_safe_X.copy()
653 if beta_loss - 1.0 < 0:
654 WH[WH < EPSILON] = EPSILON
655
656 # to avoid division by zero
657 if beta_loss - 2.0 < 0:
658 WH_safe_X_data[WH_safe_X_data < EPSILON] = EPSILON
659
660 if beta_loss == 1:
661 np.divide(X_data, WH_safe_X_data, out=WH_safe_X_data)
662 elif beta_loss == 0:
663 # speeds up computation time
664 # refer to /numpy/numpy/issues/9363
665 WH_safe_X_data **= -1
666 WH_safe_X_data **= 2
667 # element-wise multiplication
668 WH_safe_X_data *= X_data
669 else:
670 WH_safe_X_data **= beta_loss - 2
671 # element-wise multiplication
672 WH_safe_X_data *= X_data
673
674 # here numerator = dot(W.T, (dot(W, H) ** (beta_loss - 2)) * X)
675 numerator = safe_sparse_dot(W.T, WH_safe_X)
676
677 # Denominator
678 if beta_loss == 1:
679 W_sum = np.sum(W, axis=0) # shape(n_components, )
680 W_sum[W_sum == 0] = 1.0
681 denominator = W_sum[:, np.newaxis]
682
683 # beta_loss not in (1, 2)
684 else:
685 # computation of WtWH = dot(W.T, dot(W, H) ** beta_loss - 1)
686 if sp.issparse(X):
687 # memory efficient computation
688 # (compute column by column, avoiding the dense matrix WH)
689 WtWH = np.empty(H.shape)
690 for i in range(X.shape[1]):
691 WHi = np.dot(W, H[:, i])
692 if beta_loss - 1 < 0:
693 WHi[WHi < EPSILON] = EPSILON
694 WHi **= beta_loss - 1
695 WtWH[:, i] = np.dot(W.T, WHi)
696 else:
697 WH **= beta_loss - 1
698 WtWH = np.dot(W.T, WH)
699 denominator = WtWH
700
701 # Add L1 and L2 regularization
702 if l1_reg_H > 0:
703 denominator += l1_reg_H
704 if l2_reg_H > 0:
705 denominator = denominator + l2_reg_H * H
706 denominator[denominator == 0] = EPSILON
707
708 if A is not None and B is not None:
709 # Updates for the online nmf
710 if gamma != 1:
711 H **= 1 / gamma
712 numerator *= H
713 A *= rho
714 B *= rho
715 A += numerator
716 B += denominator
717 H = A / B
718
719 if gamma != 1:
720 H **= gamma
721 else:
722 delta_H = numerator
723 delta_H /= denominator
724 if gamma != 1:
725 delta_H **= gamma
726 H *= delta_H
727
728 return H
729
730
731def _fit_multiplicative_update(
732 X,
733 W,
734 H,
735 beta_loss="frobenius",
736 max_iter=200,
737 tol=1e-4,
738 l1_reg_W=0,
739 l1_reg_H=0,
740 l2_reg_W=0,
741 l2_reg_H=0,
742 update_H=True,
743 verbose=0,
744):
745 """Compute Non-negative Matrix Factorization with Multiplicative Update.
746
747 The objective function is _beta_divergence(X, WH) and is minimized with an
748 alternating minimization of W and H. Each minimization is done with a
749 Multiplicative Update.
750
751 Parameters
752 ----------
753 X : array-like of shape (n_samples, n_features)
754 Constant input matrix.
755
756 W : array-like of shape (n_samples, n_components)
757 Initial guess for the solution.
758
759 H : array-like of shape (n_components, n_features)
760 Initial guess for the solution.
761
762 beta_loss : float or {'frobenius', 'kullback-leibler', \
763 'itakura-saito'}, default='frobenius'
764 String must be in {'frobenius', 'kullback-leibler', 'itakura-saito'}.
765 Beta divergence to be minimized, measuring the distance between X
766 and the dot product WH. Note that values different from 'frobenius'
767 (or 2) and 'kullback-leibler' (or 1) lead to significantly slower
768 fits. Note that for beta_loss <= 0 (or 'itakura-saito'), the input
769 matrix X cannot contain zeros.
770
771 max_iter : int, default=200
772 Number of iterations.
773
774 tol : float, default=1e-4
775 Tolerance of the stopping condition.
776
777 l1_reg_W : float, default=0.
778 L1 regularization parameter for W.
779
780 l1_reg_H : float, default=0.
781 L1 regularization parameter for H.
782
783 l2_reg_W : float, default=0.
784 L2 regularization parameter for W.
785
786 l2_reg_H : float, default=0.
787 L2 regularization parameter for H.
788
789 update_H : bool, default=True
790 Set to True, both W and H will be estimated from initial guesses.
791 Set to False, only W will be estimated.
792
793 verbose : int, default=0
794 The verbosity level.
795
796 Returns
797 -------
798 W : ndarray of shape (n_samples, n_components)
799 Solution to the non-negative least squares problem.
800
801 H : ndarray of shape (n_components, n_features)
802 Solution to the non-negative least squares problem.
803
804 n_iter : int
805 The number of iterations done by the algorithm.
806
807 References
808 ----------
809 Lee, D. D., & Seung, H., S. (2001). Algorithms for Non-negative Matrix
810 Factorization. Adv. Neural Inform. Process. Syst.. 13.
811 Fevotte, C., & Idier, J. (2011). Algorithms for nonnegative matrix
812 factorization with the beta-divergence. Neural Computation, 23(9).
813 """
814 start_time = time.time()
815
816 beta_loss = _beta_loss_to_float(beta_loss)
817
818 # gamma for Maximization-Minimization (MM) algorithm [Fevotte 2011]
819 if beta_loss < 1:
820 gamma = 1.0 / (2.0 - beta_loss)
821 elif beta_loss > 2:
822 gamma = 1.0 / (beta_loss - 1.0)
823 else:
824 gamma = 1.0
825
826 # used for the convergence criterion
827 error_at_init = _beta_divergence(X, W, H, beta_loss, square_root=True)
828 previous_error = error_at_init
829
830 H_sum, HHt, XHt = None, None, None
831 for n_iter in range(1, max_iter + 1):
832 # update W
833 # H_sum, HHt and XHt are saved and reused if not update_H
834 W, H_sum, HHt, XHt = _multiplicative_update_w(
835 X,
836 W,
837 H,
838 beta_loss=beta_loss,
839 l1_reg_W=l1_reg_W,
840 l2_reg_W=l2_reg_W,
841 gamma=gamma,
842 H_sum=H_sum,
843 HHt=HHt,
844 XHt=XHt,
845 update_H=update_H,
846 )
847
848 # necessary for stability with beta_loss < 1
849 if beta_loss < 1:
850 W[W < np.finfo(np.float64).eps] = 0.0
851
852 # update H (only at fit or fit_transform)
853 if update_H:
854 H = _multiplicative_update_h(
855 X,
856 W,
857 H,
858 beta_loss=beta_loss,
859 l1_reg_H=l1_reg_H,
860 l2_reg_H=l2_reg_H,
861 gamma=gamma,
862 )
863
864 # These values will be recomputed since H changed
865 H_sum, HHt, XHt = None, None, None
866
867 # necessary for stability with beta_loss < 1
868 if beta_loss <= 1:
869 H[H < np.finfo(np.float64).eps] = 0.0
870
871 # test convergence criterion every 10 iterations
872 if tol > 0 and n_iter % 10 == 0:
873 error = _beta_divergence(X, W, H, beta_loss, square_root=True)
874
875 if verbose:
876 iter_time = time.time()
877 print(
878 "Epoch %02d reached after %.3f seconds, error: %f"
879 % (n_iter, iter_time - start_time, error)
880 )
881
882 if (previous_error - error) / error_at_init < tol:
883 break
884 previous_error = error
885
886 # do not print if we have already printed in the convergence test
887 if verbose and (tol == 0 or n_iter % 10 != 0):
888 end_time = time.time()
889 print(
890 "Epoch %02d reached after %.3f seconds." % (n_iter, end_time - start_time)
891 )
892
893 return W, H, n_iter
894
895
896@validate_params(
897 {
898 "X": ["array-like", "sparse matrix"],
899 "W": ["array-like", None],
900 "H": ["array-like", None],
901 "update_H": ["boolean"],
902 },
903 prefer_skip_nested_validation=False,
904)
905def non_negative_factorization(
906 X,
907 W=None,
908 H=None,
909 n_components="auto",
910 *,
911 init=None,
912 update_H=True,
913 solver="cd",
914 beta_loss="frobenius",
915 tol=1e-4,
916 max_iter=200,
917 alpha_W=0.0,
918 alpha_H="same",
919 l1_ratio=0.0,
920 random_state=None,
921 verbose=0,
922 shuffle=False,
923):
924 """Compute Non-negative Matrix Factorization (NMF).
925
926 Find two non-negative matrices (W, H) whose product approximates the non-
927 negative matrix X. This factorization can be used for example for
928 dimensionality reduction, source separation or topic extraction.
929
930 The objective function is:
931
932 .. math::
933
934 L(W, H) &= 0.5 * ||X - WH||_{loss}^2
935
936 &+ alpha\\_W * l1\\_ratio * n\\_features * ||vec(W)||_1
937
938 &+ alpha\\_H * l1\\_ratio * n\\_samples * ||vec(H)||_1
939
940 &+ 0.5 * alpha\\_W * (1 - l1\\_ratio) * n\\_features * ||W||_{Fro}^2
941
942 &+ 0.5 * alpha\\_H * (1 - l1\\_ratio) * n\\_samples * ||H||_{Fro}^2,
943
944 where :math:`||A||_{Fro}^2 = \\sum_{i,j} A_{ij}^2` (Frobenius norm) and
945 :math:`||vec(A)||_1 = \\sum_{i,j} abs(A_{ij})` (Elementwise L1 norm)
946
947 The generic norm :math:`||X - WH||_{loss}^2` may represent
948 the Frobenius norm or another supported beta-divergence loss.
949 The choice between options is controlled by the `beta_loss` parameter.
950
951 The regularization terms are scaled by `n_features` for `W` and by `n_samples` for
952 `H` to keep their impact balanced with respect to one another and to the data fit
953 term as independent as possible of the size `n_samples` of the training set.
954
955 The objective function is minimized with an alternating minimization of W
956 and H. If H is given and update_H=False, it solves for W only.
957
958 Note that the transformed data is named W and the components matrix is named H. In
959 the NMF literature, the naming convention is usually the opposite since the data
960 matrix X is transposed.
961
962 Parameters
963 ----------
964 X : {array-like, sparse matrix} of shape (n_samples, n_features)
965 Constant matrix.
966
967 W : array-like of shape (n_samples, n_components), default=None
968 If `init='custom'`, it is used as initial guess for the solution.
969 If `update_H=False`, it is initialised as an array of zeros, unless
970 `solver='mu'`, then it is filled with values calculated by
971 `np.sqrt(X.mean() / self._n_components)`.
972 If `None`, uses the initialisation method specified in `init`.
973
974 H : array-like of shape (n_components, n_features), default=None
975 If `init='custom'`, it is used as initial guess for the solution.
976 If `update_H=False`, it is used as a constant, to solve for W only.
977 If `None`, uses the initialisation method specified in `init`.
978
979 n_components : int or {'auto'} or None, default='auto'
980 Number of components. If `None`, all features are kept.
981 If `n_components='auto'`, the number of components is automatically inferred
982 from `W` or `H` shapes.
983
984 .. versionchanged:: 1.4
985 Added `'auto'` value.
986
987 .. versionchanged:: 1.6
988 Default value changed from `None` to `'auto'`.
989
990 init : {'random', 'nndsvd', 'nndsvda', 'nndsvdar', 'custom'}, default=None
991 Method used to initialize the procedure.
992
993 Valid options:
994
995 - None: 'nndsvda' if n_components < n_features, otherwise 'random'.
996 - 'random': non-negative random matrices, scaled with:
997 `sqrt(X.mean() / n_components)`
998 - 'nndsvd': Nonnegative Double Singular Value Decomposition (NNDSVD)
999 initialization (better for sparseness)
1000 - 'nndsvda': NNDSVD with zeros filled with the average of X
1001 (better when sparsity is not desired)
1002 - 'nndsvdar': NNDSVD with zeros filled with small random values
1003 (generally faster, less accurate alternative to NNDSVDa
1004 for when sparsity is not desired)
1005 - 'custom': If `update_H=True`, use custom matrices W and H which must both
1006 be provided. If `update_H=False`, then only custom matrix H is used.
1007
1008 .. versionchanged:: 0.23
1009 The default value of `init` changed from 'random' to None in 0.23.
1010
1011 .. versionchanged:: 1.1
1012 When `init=None` and n_components is less than n_samples and n_features
1013 defaults to `nndsvda` instead of `nndsvd`.
1014
1015 update_H : bool, default=True
1016 Set to True, both W and H will be estimated from initial guesses.
1017 Set to False, only W will be estimated.
1018
1019 solver : {'cd', 'mu'}, default='cd'
1020 Numerical solver to use:
1021
1022 - 'cd' is a Coordinate Descent solver that uses Fast Hierarchical
1023 Alternating Least Squares (Fast HALS).
1024 - 'mu' is a Multiplicative Update solver.
1025
1026 .. versionadded:: 0.17
1027 Coordinate Descent solver.
1028
1029 .. versionadded:: 0.19
1030 Multiplicative Update solver.
1031
1032 beta_loss : float or {'frobenius', 'kullback-leibler', \
1033 'itakura-saito'}, default='frobenius'
1034 Beta divergence to be minimized, measuring the distance between X
1035 and the dot product WH. Note that values different from 'frobenius'
1036 (or 2) and 'kullback-leibler' (or 1) lead to significantly slower
1037 fits. Note that for beta_loss <= 0 (or 'itakura-saito'), the input
1038 matrix X cannot contain zeros. Used only in 'mu' solver.
1039
1040 .. versionadded:: 0.19
1041
1042 tol : float, default=1e-4
1043 Tolerance of the stopping condition.
1044
1045 max_iter : int, default=200
1046 Maximum number of iterations before timing out.
1047
1048 alpha_W : float, default=0.0
1049 Constant that multiplies the regularization terms of `W`. Set it to zero
1050 (default) to have no regularization on `W`.
1051
1052 .. versionadded:: 1.0
1053
1054 alpha_H : float or "same", default="same"
1055 Constant that multiplies the regularization terms of `H`. Set it to zero to
1056 have no regularization on `H`. If "same" (default), it takes the same value as
1057 `alpha_W`.
1058
1059 .. versionadded:: 1.0
1060
1061 l1_ratio : float, default=0.0
1062 The regularization mixing parameter, with 0 <= l1_ratio <= 1.
1063 For l1_ratio = 0 the penalty is an elementwise L2 penalty
1064 (aka Frobenius Norm).
1065 For l1_ratio = 1 it is an elementwise L1 penalty.
1066 For 0 < l1_ratio < 1, the penalty is a combination of L1 and L2.
1067
1068 random_state : int, RandomState instance or None, default=None
1069 Used for NMF initialisation (when ``init`` == 'nndsvdar' or
1070 'random'), and in Coordinate Descent. Pass an int for reproducible
1071 results across multiple function calls.
1072 See :term:`Glossary <random_state>`.
1073
1074 verbose : int, default=0
1075 The verbosity level.
1076
1077 shuffle : bool, default=False
1078 If true, randomize the order of coordinates in the CD solver.
1079
1080 Returns
1081 -------
1082 W : ndarray of shape (n_samples, n_components)
1083 Solution to the non-negative least squares problem.
1084
1085 H : ndarray of shape (n_components, n_features)
1086 Solution to the non-negative least squares problem.
1087
1088 n_iter : int
1089 Actual number of iterations.
1090
1091 References
1092 ----------
1093 .. [1] :doi:`"Fast local algorithms for large scale nonnegative matrix and tensor
1094 factorizations" <10.1587/transfun.E92.A.708>`
1095 Cichocki, Andrzej, and P. H. A. N. Anh-Huy. IEICE transactions on fundamentals
1096 of electronics, communications and computer sciences 92.3: 708-721, 2009.
1097
1098 .. [2] :doi:`"Algorithms for nonnegative matrix factorization with the
1099 beta-divergence" <10.1162/NECO_a_00168>`
1100 Fevotte, C., & Idier, J. (2011). Neural Computation, 23(9).
1101
1102 Examples
1103 --------
1104 >>> import numpy as np
1105 >>> X = np.array([[1,1], [2, 1], [3, 1.2], [4, 1], [5, 0.8], [6, 1]])
1106 >>> from sklearn.decomposition import non_negative_factorization
1107 >>> W, H, n_iter = non_negative_factorization(
1108 ... X, n_components=2, init='random', random_state=0)
1109 """
1110 est = NMF(
1111 n_components=n_components,
1112 init=init,
1113 solver=solver,
1114 beta_loss=beta_loss,
1115 tol=tol,
1116 max_iter=max_iter,
1117 random_state=random_state,
1118 alpha_W=alpha_W,
1119 alpha_H=alpha_H,
1120 l1_ratio=l1_ratio,
1121 verbose=verbose,
1122 shuffle=shuffle,
1123 )
1124 est._validate_params()
1125
1126 X = check_array(X, accept_sparse=("csr", "csc"), dtype=[np.float64, np.float32])
1127
1128 with config_context(assume_finite=True):
1129 W, H, n_iter = est._fit_transform(X, W=W, H=H, update_H=update_H)
1130
1131 return W, H, n_iter
1132
1133
1134class _BaseNMF(ClassNamePrefixFeaturesOutMixin, TransformerMixin, BaseEstimator, ABC):
1135 """Base class for NMF and MiniBatchNMF."""
1136
1137 _parameter_constraints: dict = {
1138 "n_components": [
1139 Interval(Integral, 1, None, closed="left"),
1140 None,
1141 StrOptions({"auto"}),
1142 ],
1143 "init": [
1144 StrOptions({"random", "nndsvd", "nndsvda", "nndsvdar", "custom"}),
1145 None,
1146 ],
1147 "beta_loss": [
1148 StrOptions({"frobenius", "kullback-leibler", "itakura-saito"}),
1149 Real,
1150 ],
1151 "tol": [Interval(Real, 0, None, closed="left")],
1152 "max_iter": [Interval(Integral, 1, None, closed="left")],
1153 "random_state": ["random_state"],
1154 "alpha_W": [Interval(Real, 0, None, closed="left")],
1155 "alpha_H": [Interval(Real, 0, None, closed="left"), StrOptions({"same"})],
1156 "l1_ratio": [Interval(Real, 0, 1, closed="both")],
1157 "verbose": ["verbose"],
1158 }
1159
1160 def __init__(
1161 self,
1162 n_components="auto",
1163 *,
1164 init=None,
1165 beta_loss="frobenius",
1166 tol=1e-4,
1167 max_iter=200,
1168 random_state=None,
1169 alpha_W=0.0,
1170 alpha_H="same",
1171 l1_ratio=0.0,
1172 verbose=0,
1173 ):
1174 self.n_components = n_components
1175 self.init = init
1176 self.beta_loss = beta_loss
1177 self.tol = tol
1178 self.max_iter = max_iter
1179 self.random_state = random_state
1180 self.alpha_W = alpha_W
1181 self.alpha_H = alpha_H
1182 self.l1_ratio = l1_ratio
1183 self.verbose = verbose
1184
1185 def _check_params(self, X):
1186 # n_components
1187 self._n_components = self.n_components
1188 if self._n_components is None:
1189 self._n_components = X.shape[1]
1190
1191 # beta_loss
1192 self._beta_loss = _beta_loss_to_float(self.beta_loss)
1193
1194 def _check_w_h(self, X, W, H, update_H):
1195 """Check W and H, or initialize them."""
1196 n_samples, n_features = X.shape
1197
1198 if self.init == "custom" and update_H:
1199 _check_init(H, (self._n_components, n_features), "NMF (input H)")
1200 _check_init(W, (n_samples, self._n_components), "NMF (input W)")