""" Dictionary learning. """ # Author: Vlad Niculae, Gael Varoquaux, Alexandre Gramfort # License: BSD 3 clause import itertools import sys import time import warnings from math import ceil from numbers import Integral, Real import numpy as np from joblib import effective_n_jobs from scipy import linalg from ..base import ( BaseEstimator, ClassNamePrefixFeaturesOutMixin, TransformerMixin, _fit_context, ) from ..linear_model import Lars, Lasso, LassoLars, orthogonal_mp_gram from ..utils import check_array, check_random_state, gen_batches, gen_even_slices from ..utils._param_validation import Hidden, Interval, StrOptions, validate_params from ..utils.extmath import randomized_svd, row_norms, svd_flip from ..utils.parallel import Parallel, delayed from ..utils.validation import check_is_fitted def _check_positive_coding(method, positive): if positive and method in ["omp", "lars"]: raise ValueError( "Positive constraint not supported for '{}' coding method.".format(method) ) def _sparse_encode_precomputed( X, dictionary, *, gram=None, cov=None, algorithm="lasso_lars", regularization=None, copy_cov=True, init=None, max_iter=1000, verbose=0, positive=False, ): """Generic sparse coding with precomputed Gram and/or covariance matrices. Each row of the result is the solution to a Lasso problem. Parameters ---------- X : ndarray of shape (n_samples, n_features) Data matrix. dictionary : ndarray of shape (n_components, n_features) The dictionary matrix against which to solve the sparse coding of the data. Some of the algorithms assume normalized rows. gram : ndarray of shape (n_components, n_components), default=None Precomputed Gram matrix, `dictionary * dictionary'` gram can be `None` if method is 'threshold'. cov : ndarray of shape (n_components, n_samples), default=None Precomputed covariance, `dictionary * X'`. algorithm : {'lasso_lars', 'lasso_cd', 'lars', 'omp', 'threshold'}, \ default='lasso_lars' The algorithm used: * `'lars'`: uses the least angle regression method (`linear_model.lars_path`); * `'lasso_lars'`: uses Lars to compute the Lasso solution; * `'lasso_cd'`: uses the coordinate descent method to compute the Lasso solution (`linear_model.Lasso`). lasso_lars will be faster if the estimated components are sparse; * `'omp'`: uses orthogonal matching pursuit to estimate the sparse solution; * `'threshold'`: squashes to zero all coefficients less than regularization from the projection `dictionary * data'`. regularization : int or float, default=None The regularization parameter. It corresponds to alpha when algorithm is `'lasso_lars'`, `'lasso_cd'` or `'threshold'`. Otherwise it corresponds to `n_nonzero_coefs`. init : ndarray of shape (n_samples, n_components), default=None Initialization value of the sparse code. Only used if `algorithm='lasso_cd'`. max_iter : int, default=1000 Maximum number of iterations to perform if `algorithm='lasso_cd'` or `'lasso_lars'`. copy_cov : bool, default=True Whether to copy the precomputed covariance matrix; if `False`, it may be overwritten. verbose : int, default=0 Controls the verbosity; the higher, the more messages. positive: bool, default=False Whether to enforce a positivity constraint on the sparse code. .. versionadded:: 0.20 Returns ------- code : ndarray of shape (n_components, n_features) The sparse codes. """ n_samples, n_features = X.shape n_components = dictionary.shape[0] if algorithm == "lasso_lars": alpha = float(regularization) / n_features # account for scaling try: err_mgt = np.seterr(all="ignore") # Not passing in verbose=max(0, verbose-1) because Lars.fit already # corrects the verbosity level. lasso_lars = LassoLars( alpha=alpha, fit_intercept=False, verbose=verbose, precompute=gram, fit_path=False, positive=positive, max_iter=max_iter, ) lasso_lars.fit(dictionary.T, X.T, Xy=cov) new_code = lasso_lars.coef_ finally: np.seterr(**err_mgt) elif algorithm == "lasso_cd": alpha = float(regularization) / n_features # account for scaling # TODO: Make verbosity argument for Lasso? # sklearn.linear_model.coordinate_descent.enet_path has a verbosity # argument that we could pass in from Lasso. clf = Lasso( alpha=alpha, fit_intercept=False, precompute=gram, max_iter=max_iter, warm_start=True, positive=positive, ) if init is not None: # In some workflows using coordinate descent algorithms: # - users might provide NumPy arrays with read-only buffers # - `joblib` might memmap arrays making their buffer read-only # TODO: move this handling (which is currently too broad) # closer to the actual private function which need buffers to be writable. if not init.flags["WRITEABLE"]: init = np.array(init) clf.coef_ = init clf.fit(dictionary.T, X.T, check_input=False) new_code = clf.coef_ elif algorithm == "lars": try: err_mgt = np.seterr(all="ignore") # Not passing in verbose=max(0, verbose-1) because Lars.fit already # corrects the verbosity level. lars = Lars( fit_intercept=False, verbose=verbose, precompute=gram, n_nonzero_coefs=int(regularization), fit_path=False, ) lars.fit(dictionary.T, X.T, Xy=cov) new_code = lars.coef_ finally: np.seterr(**err_mgt) elif algorithm == "threshold": new_code = (np.sign(cov) * np.maximum(np.abs(cov) - regularization, 0)).T if positive: np.clip(new_code, 0, None, out=new_code) elif algorithm == "omp": new_code = orthogonal_mp_gram( Gram=gram, Xy=cov, n_nonzero_coefs=int(regularization), tol=None, norms_squared=row_norms(X, squared=True), copy_Xy=copy_cov, ).T return new_code.reshape(n_samples, n_components) @validate_params( { "X": ["array-like"], "dictionary": ["array-like"], "gram": ["array-like", None], "cov": ["array-like", None], "algorithm": [ StrOptions({"lasso_lars", "lasso_cd", "lars", "omp", "threshold"}) ], "n_nonzero_coefs": [Interval(Integral, 1, None, closed="left"), None], "alpha": [Interval(Real, 0, None, closed="left"), None], "copy_cov": ["boolean"], "init": ["array-like", None], "max_iter": [Interval(Integral, 0, None, closed="left")], "n_jobs": [Integral, None], "check_input": ["boolean"], "verbose": ["verbose"], "positive": ["boolean"], }, prefer_skip_nested_validation=True, ) # XXX : could be moved to the linear_model module def sparse_encode( X, dictionary, *, gram=None, cov=None, algorithm="lasso_lars", n_nonzero_coefs=None, alpha=None, copy_cov=True, init=None, max_iter=1000, n_jobs=None, check_input=True, verbose=0, positive=False, ): """Sparse coding. Each row of the result is the solution to a sparse coding problem. The goal is to find a sparse array `code` such that:: X ~= code * dictionary Read more in the :ref:`User Guide `. Parameters ---------- X : array-like of shape (n_samples, n_features) Data matrix. dictionary : array-like of shape (n_components, n_features) The dictionary matrix against which to solve the sparse coding of the data. Some of the algorithms assume normalized rows for meaningful output. gram : array-like of shape (n_components, n_components), default=None Precomputed Gram matrix, `dictionary * dictionary'`. cov : array-like of shape (n_components, n_samples), default=None Precomputed covariance, `dictionary' * X`. algorithm : {'lasso_lars', 'lasso_cd', 'lars', 'omp', 'threshold'}, \ default='lasso_lars' The algorithm used: * `'lars'`: uses the least angle regression method (`linear_model.lars_path`); * `'lasso_lars'`: uses Lars to compute the Lasso solution; * `'lasso_cd'`: uses the coordinate descent method to compute the Lasso solution (`linear_model.Lasso`). lasso_lars will be faster if the estimated components are sparse; * `'omp'`: uses orthogonal matching pursuit to estimate the sparse solution; * `'threshold'`: squashes to zero all coefficients less than regularization from the projection `dictionary * data'`. n_nonzero_coefs : int, default=None Number of nonzero coefficients to target in each column of the solution. This is only used by `algorithm='lars'` and `algorithm='omp'` and is overridden by `alpha` in the `omp` case. If `None`, then `n_nonzero_coefs=int(n_features / 10)`. alpha : float, default=None If `algorithm='lasso_lars'` or `algorithm='lasso_cd'`, `alpha` is the penalty applied to the L1 norm. If `algorithm='threshold'`, `alpha` is the absolute value of the threshold below which coefficients will be squashed to zero. If `algorithm='omp'`, `alpha` is the tolerance parameter: the value of the reconstruction error targeted. In this case, it overrides `n_nonzero_coefs`. If `None`, default to 1. copy_cov : bool, default=True Whether to copy the precomputed covariance matrix; if `False`, it may be overwritten. init : ndarray of shape (n_samples, n_components), default=None Initialization value of the sparse codes. Only used if `algorithm='lasso_cd'`. max_iter : int, default=1000 Maximum number of iterations to perform if `algorithm='lasso_cd'` or `'lasso_lars'`. n_jobs : int, default=None Number of parallel jobs to run. ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. ``-1`` means using all processors. See :term:`Glossary ` for more details. check_input : bool, default=True If `False`, the input arrays X and dictionary will not be checked. verbose : int, default=0 Controls the verbosity; the higher, the more messages. positive : bool, default=False Whether to enforce positivity when finding the encoding. .. versionadded:: 0.20 Returns ------- code : ndarray of shape (n_samples, n_components) The sparse codes. See Also -------- sklearn.linear_model.lars_path : Compute Least Angle Regression or Lasso path using LARS algorithm. sklearn.linear_model.orthogonal_mp : Solves Orthogonal Matching Pursuit problems. sklearn.linear_model.Lasso : Train Linear Model with L1 prior as regularizer. SparseCoder : Find a sparse representation of data from a fixed precomputed dictionary. """ if check_input: if algorithm == "lasso_cd": dictionary = check_array( dictionary, order="C", dtype=[np.float64, np.float32] ) X = check_array(X, order="C", dtype=[np.float64, np.float32]) else: dictionary = check_array(dictionary) X = check_array(X) if dictionary.shape[1] != X.shape[1]: raise ValueError( "Dictionary and X have different numbers of features:" "dictionary.shape: {} X.shape{}".format(dictionary.shape, X.shape) ) _check_positive_coding(algorithm, positive) return _sparse_encode( X, dictionary, gram=gram, cov=cov, algorithm=algorithm, n_nonzero_coefs=n_nonzero_coefs, alpha=alpha, copy_cov=copy_cov, init=init, max_iter=max_iter, n_jobs=n_jobs, verbose=verbose, positive=positive, ) def _sparse_encode( X, dictionary, *, gram=None, cov=None, algorithm="lasso_lars", n_nonzero_coefs=None, alpha=None, copy_cov=True, init=None, max_iter=1000, n_jobs=None, verbose=0, positive=False, ): """Sparse coding without input/parameter validation.""" n_samples, n_features = X.shape n_components = dictionary.shape[0] if algorithm in ("lars", "omp"): regularization = n_nonzero_coefs if regularization is None: regularization = min(max(n_features / 10, 1), n_components) else: regularization = alpha if regularization is None: regularization = 1.0 if gram is None and algorithm != "threshold": gram = np.dot(dictionary, dictionary.T) if cov is None and algorithm != "lasso_cd": copy_cov = False cov = np.dot(dictionary, X.T) if effective_n_jobs(n_jobs) == 1 or algorithm == "threshold": code = _sparse_encode_precomputed( X, dictionary, gram=gram, cov=cov, algorithm=algorithm, regularization=regularization, copy_cov=copy_cov, init=init, max_iter=max_iter, verbose=verbose, positive=positive, ) return code # Enter parallel code block n_samples = X.shape[0] n_components = dictionary.shape[0] code = np.empty((n_samples, n_components)) slices = list(gen_even_slices(n_samples, effective_n_jobs(n_jobs))) code_views = Parallel(n_jobs=n_jobs, verbose=verbose)( delayed(_sparse_encode_precomputed)( X[this_slice], dictionary, gram=gram, cov=cov[:, this_slice] if cov is not None else None, algorithm=algorithm, regularization=regularization, copy_cov=copy_cov, init=init[this_slice] if init is not None else None, max_iter=max_iter, verbose=verbose, positive=positive, ) for this_slice in slices ) for this_slice, this_view in zip(slices, code_views): code[this_slice] = this_view return code def _update_dict( dictionary, Y, code, A=None, B=None, verbose=False, random_state=None, positive=False, ): """Update the dense dictionary factor in place. Parameters ---------- dictionary : ndarray of shape (n_components, n_features) Value of the dictionary at the previous iteration. Y : ndarray of shape (n_samples, n_features) Data matrix. code : ndarray of shape (n_samples, n_components) Sparse coding of the data against which to optimize the dictionary. A : ndarray of shape (n_components, n_components), default=None Together with `B`, sufficient stats of the online model to update the dictionary. B : ndarray of shape (n_features, n_components), default=None Together with `A`, sufficient stats of the online model to update the dictionary. verbose: bool, default=False Degree of output the procedure will print. random_state : int, RandomState instance or None, default=None Used for randomly initializing the dictionary. Pass an int for reproducible results across multiple function calls. See :term:`Glossary `. positive : bool, default=False Whether to enforce positivity when finding the dictionary. .. versionadded:: 0.20 """ n_samples, n_components = code.shape random_state = check_random_state(random_state) if A is None: A = code.T @ code if B is None: B = Y.T @ code n_unused = 0 for k in range(n_components): if A[k, k] > 1e-6: # 1e-6 is arbitrary but consistent with the spams implementation dictionary[k] += (B[:, k] - A[k] @ dictionary) / A[k, k] else: # kth atom is almost never used -> sample a new one from the data newd = Y[random_state.choice(n_samples)] # add small noise to avoid making the sparse coding ill conditioned noise_level = 0.01 * (newd.std() or 1) # avoid 0 std noise = random_state.normal(0, noise_level, size=len(newd)) dictionary[k] = newd + noise code[:, k] = 0 n_unused += 1 if positive: np.clip(dictionary[k], 0, None, out=dictionary[k]) # Projection on the constraint set ||V_k|| <= 1 dictionary[k] /= max(linalg.norm(dictionary[k]), 1) if verbose and n_unused > 0: print(f"{n_unused} unused atoms resampled.") def _dict_learning( X, n_components, *, alpha, max_iter, tol, method, n_jobs, dict_init, code_init, callback, verbose, random_state, return_n_iter, positive_dict, positive_code, method_max_iter, ): """Main dictionary learning algorithm""" t0 = time.time() # Init the code and the dictionary with SVD of Y if code_init is not None and dict_init is not None: code = np.array(code_init, order="F") # Don't copy V, it will happen below dictionary = dict_init else: code, S, dictionary = linalg.svd(X, full_matrices=False) # flip the initial code's sign to enforce deterministic output code, dictionary = svd_flip(code, dictionary) dictionary = S[:, np.newaxis] * dictionary r = len(dictionary) if n_components <= r: # True even if n_components=None code = code[:, :n_components] dictionary = dictionary[:n_components, :] else: code = np.c_[code, np.zeros((len(code), n_components - r))] dictionary = np.r_[ dictionary, np.zeros((n_components - r, dictionary.shape[1])) ] # Fortran-order dict better suited for the sparse coding which is the # bottleneck of this algorithm. dictionary = np.asfortranarray(dictionary) errors = [] current_cost = np.nan if verbose == 1: print("[dict_learning]", end=" ") # If max_iter is 0, number of iterations returned should be zero ii = -1 for ii in range(max_iter): dt = time.time() - t0 if verbose == 1: sys.stdout.write(".") sys.stdout.flush() elif verbose: print( "Iteration % 3i (elapsed time: % 3is, % 4.1fmn, current cost % 7.3f)" % (ii, dt, dt / 60, current_cost) ) # Update code code = sparse_encode( X, dictionary, algorithm=method, alpha=alpha, init=code, n_jobs=n_jobs, positive=positive_code, max_iter=method_max_iter, verbose=verbose, ) # Update dictionary in place _update_dict( dictionary, X, code, verbose=verbose, random_state=random_state, positive=positive_dict, ) # Cost function current_cost = 0.5 * np.sum((X - code @ dictionary) ** 2) + alpha * np.sum( np.abs(code) ) errors.append(current_cost) if ii > 0: dE = errors[-2] - errors[-1] # assert(dE >= -tol * errors[-1]) if dE < tol * errors[-1]: if verbose == 1: # A line return print("") elif verbose: print("--- Convergence reached after %d iterations" % ii) break if ii % 5 == 0 and callback is not None: callback(locals()) if return_n_iter: return code, dictionary, errors, ii + 1 else: return code, dictionary, errors def _check_warn_deprecated(param, name, default, additional_message=None): if param != "deprecated": msg = ( f"'{name}' is deprecated in version 1.1 and will be removed in version 1.4." ) if additional_message: msg += f" {additional_message}" warnings.warn(msg, FutureWarning) return param else: return default def dict_learning_online( X, n_components=2, *, alpha=1, n_iter="deprecated", max_iter=None, return_code=True, dict_init=None, callback=None, batch_size=256, verbose=False, shuffle=True, n_jobs=None, method="lars", iter_offset="deprecated", random_state=None, return_inner_stats="deprecated", inner_stats="deprecated", return_n_iter="deprecated", positive_dict=False, positive_code=False, method_max_iter=1000, tol=1e-3, max_no_improvement=10, ): """Solve a dictionary learning matrix factorization problem online. Finds the best dictionary and the corresponding sparse code for approximating the data matrix X by solving:: (U^*, V^*) = argmin 0.5 || X - U V ||_Fro^2 + alpha * || U ||_1,1 (U,V) with || V_k ||_2 = 1 for all 0 <= k < n_components where V is the dictionary and U is the sparse code. ||.||_Fro stands for the Frobenius norm and ||.||_1,1 stands for the entry-wise matrix norm which is the sum of the absolute values of all the entries in the matrix. This is accomplished by repeatedly iterating over mini-batches by slicing the input data. Read more in the :ref:`User Guide `. Parameters ---------- X : ndarray of shape (n_samples, n_features) Data matrix. n_components : int or None, default=2 Number of dictionary atoms to extract. If None, then ``n_components`` is set to ``n_features``. alpha : float, default=1 Sparsity controlling parameter. n_iter : int, default=100 Number of mini-batch iterations to perform. .. deprecated:: 1.1 `n_iter` is deprecated in 1.1 and will be removed in 1.4. Use `max_iter` instead. max_iter : int, default=None Maximum number of iterations over the complete dataset before stopping independently of any early stopping criterion heuristics. If ``max_iter`` is not None, ``n_iter`` is ignored. .. versionadded:: 1.1 return_code : bool, default=True Whether to also return the code U or just the dictionary `V`. dict_init : ndarray of shape (n_components, n_features), default=None Initial values for the dictionary for warm restart scenarios. If `None`, the initial values for the dictionary are created with an SVD decomposition of the data via :func:`~sklearn.utils.extmath.randomized_svd`. callback : callable, default=None A callable that gets invoked at the end of each iteration. batch_size : int, default=256 The number of samples to take in each batch. .. versionchanged:: 1.3 The default value of `batch_size` changed from 3 to 256 in version 1.3. verbose : bool, default=False To control the verbosity of the procedure. shuffle : bool, default=True Whether to shuffle the data before splitting it in batches. n_jobs : int, default=None Number of parallel jobs to run. ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. ``-1`` means using all processors. See :term:`Glossary ` for more details. method : {'lars', 'cd'}, default='lars' * `'lars'`: uses the least angle regression method to solve the lasso problem (`linear_model.lars_path`); * `'cd'`: uses the coordinate descent method to compute the Lasso solution (`linear_model.Lasso`). Lars will be faster if the estimated components are sparse. iter_offset : int, default=0 Number of previous iterations completed on the dictionary used for initialization. .. deprecated:: 1.1 `iter_offset` serves internal purpose only and will be removed in 1.4. random_state : int, RandomState instance or None, default=None Used for initializing the dictionary when ``dict_init`` is not specified, randomly shuffling the data when ``shuffle`` is set to ``True``, and updating the dictionary. Pass an int for reproducible results across multiple function calls. See :term:`Glossary `. return_inner_stats : bool, default=False Return the inner statistics A (dictionary covariance) and B (data approximation). Useful to restart the algorithm in an online setting. If `return_inner_stats` is `True`, `return_code` is ignored. .. deprecated:: 1.1 `return_inner_stats` serves internal purpose only and will be removed in 1.4. inner_stats : tuple of (A, B) ndarrays, default=None Inner sufficient statistics that are kept by the algorithm. Passing them at initialization is useful in online settings, to avoid losing the history of the evolution. `A` `(n_components, n_components)` is the dictionary covariance matrix. `B` `(n_features, n_components)` is the data approximation matrix. .. deprecated:: 1.1 `inner_stats` serves internal purpose only and will be removed in 1.4. return_n_iter : bool, default=False Whether or not to return the number of iterations. .. deprecated:: 1.1 `return_n_iter` will be removed in 1.4 and n_iter will never be returned. positive_dict : bool, default=False Whether to enforce positivity when finding the dictionary. .. versionadded:: 0.20 positive_code : bool, default=False Whether to enforce positivity when finding the code. .. versionadded:: 0.20 method_max_iter : int, default=1000 Maximum number of iterations to perform when solving the lasso problem. .. versionadded:: 0.22 tol : float, default=1e-3 Control early stopping based on the norm of the differences in the dictionary between 2 steps. Used only if `max_iter` is not None. To disable early stopping based on changes in the dictionary, set `tol` to 0.0. .. versionadded:: 1.1 max_no_improvement : int, default=10 Control early stopping based on the consecutive number of mini batches that does not yield an improvement on the smoothed cost function. Used only if `max_iter` is not None. To disable convergence detection based on cost function, set `max_no_improvement` to None. .. versionadded:: 1.1 Returns ------- code : ndarray of shape (n_samples, n_components), The sparse code (only returned if `return_code=True`). dictionary : ndarray of shape (n_components, n_features), The solutions to the dictionary learning problem. n_iter : int Number of iterations run. Returned only if `return_n_iter` is set to `True`. See Also -------- dict_learning : Solve a dictionary learning matrix factorization problem. DictionaryLearning : Find a dictionary that sparsely encodes data. MiniBatchDictionaryLearning : A faster, less accurate, version of the dictionary learning algorithm. SparsePCA : Sparse Principal Components Analysis. MiniBatchSparsePCA : Mini-batch Sparse Principal Components Analysis. """ deps = (return_n_iter, return_inner_stats, iter_offset, inner_stats) if max_iter is not None and not all(arg == "deprecated" for arg in deps): raise ValueError( "The following arguments are incompatible with 'max_iter': " "return_n_iter, return_inner_stats, iter_offset, inner_stats" ) iter_offset = _check_warn_deprecated(iter_offset, "iter_offset", default=0) return_inner_stats = _check_warn_deprecated( return_inner_stats, "return_inner_stats", default=False, additional_message="From 1.4 inner_stats will never be returned.", ) inner_stats = _check_warn_deprecated(inner_stats, "inner_stats", default=None) return_n_iter = _check_warn_deprecated( return_n_iter, "return_n_iter", default=False, additional_message=( "From 1.4 'n_iter' will never be returned. Refer to the 'n_iter_' and " "'n_steps_' attributes of the MiniBatchDictionaryLearning object instead." ), ) if max_iter is not None: transform_algorithm = "lasso_" + method est = MiniBatchDictionaryLearning( n_components=n_components, alpha=alpha, n_iter=n_iter, n_jobs=n_jobs, fit_algorithm=method, batch_size=batch_size, shuffle=shuffle, dict_init=dict_init, random_state=random_state, transform_algorithm=transform_algorithm, transform_alpha=alpha, positive_code=positive_code, positive_dict=positive_dict, transform_max_iter=method_max_iter, verbose=verbose, callback=callback, tol=tol, max_no_improvement=max_no_improvement, ).fit(X) if not return_code: return est.components_ else: code = est.transform(X) return code, est.components_ # TODO(1.4) remove the whole old behavior # Fallback to old behavior n_iter = _check_warn_deprecated( n_iter, "n_iter", default=100, additional_message="Use 'max_iter' instead." ) if n_components is None: n_components = X.shape[1] if method not in ("lars", "cd"): raise ValueError("Coding method not supported as a fit algorithm.") _check_positive_coding(method, positive_code) method = "lasso_" + method t0 = time.time() n_samples, n_features = X.shape # Avoid integer division problems alpha = float(alpha) random_state = check_random_state(random_state) # Init V with SVD of X if dict_init is not None: dictionary = dict_init else: _, S, dictionary = randomized_svd(X, n_components, random_state=random_state) dictionary = S[:, np.newaxis] * dictionary r = len(dictionary) if n_components <= r: dictionary = dictionary[:n_components, :] else: dictionary = np.r_[ dictionary, np.zeros((n_components - r, dictionary.shape[1]), dtype=dictionary.dtype), ] if verbose == 1: print("[dict_learning]", end=" ") if shuffle: X_train = X.copy() random_state.shuffle(X_train) else: X_train = X X_train = check_array( X_train, order="C", dtype=[np.float64, np.float32], copy=False ) # Fortran-order dict better suited for the sparse coding which is the # bottleneck of this algorithm. dictionary = check_array(dictionary, order="F", dtype=X_train.dtype, copy=False) dictionary = np.require(dictionary, requirements="W") batches = gen_batches(n_samples, batch_size) batches = itertools.cycle(batches) # The covariance of the dictionary if inner_stats is None: A = np.zeros((n_components, n_components), dtype=X_train.dtype) # The data approximation B = np.zeros((n_features, n_components), dtype=X_train.dtype) else: A = inner_stats[0].copy() B = inner_stats[1].copy() # If n_iter is zero, we need to return zero. ii = iter_offset - 1 for ii, batch in zip(range(iter_offset, iter_offset + n_iter), batches): this_X = X_train[batch] dt = time.time() - t0 if verbose == 1: sys.stdout.write(".") sys.stdout.flush() elif verbose: if verbose > 10 or ii % ceil(100.0 / verbose) == 0: print( "Iteration % 3i (elapsed time: % 3is, % 4.1fmn)" % (ii, dt, dt / 60) ) this_code = sparse_encode( this_X, dictionary, algorithm=method, alpha=alpha, n_jobs=n_jobs, check_input=False, positive=positive_code, max_iter=method_max_iter, verbose=verbose, ) # Update the auxiliary variables if ii < batch_size - 1: theta = float((ii + 1) * batch_size) else: theta = float(batch_size**2 + ii + 1 - batch_size) beta = (theta + 1 - batch_size) / (theta + 1) A *= beta A += np.dot(this_code.T, this_code) B *= beta B += np.dot(this_X.T, this_code) # Update dictionary in place _update_dict( dictionary, this_X, this_code, A, B, verbose=verbose, random_state=random_state, positive=positive_dict, ) # Maybe we need a stopping criteria based on the amount of # modification in the dictionary if callback is not None: callback(locals()) if return_inner_stats: if return_n_iter: return dictionary, (A, B), ii - iter_offset + 1 else: return dictionary, (A, B) if return_code: if verbose > 1: print("Learning code...", end=" ") elif verbose == 1: print("|", end=" ") code = sparse_encode( X, dictionary, algorithm=method, alpha=alpha, n_jobs=n_jobs, check_input=False, positive=positive_code, max_iter=method_max_iter, verbose=verbose, ) if verbose > 1: dt = time.time() - t0 print("done (total time: % 3is, % 4.1fmn)" % (dt, dt / 60)) if return_n_iter: return code, dictionary, ii - iter_offset + 1 else: return code, dictionary if return_n_iter: return dictionary, ii - iter_offset + 1 else: return dictionary @validate_params( { "X": ["array-like"], "method": [StrOptions({"lars", "cd"})], "return_n_iter": ["boolean"], "method_max_iter": [Interval(Integral, 0, None, closed="left")], }, prefer_skip_nested_validation=False, ) def dict_learning( X, n_components, *, alpha, max_iter=100, tol=1e-8, method="lars", n_jobs=None, dict_init=None, code_init=None, callback=None, verbose=False, random_state=None, return_n_iter=False, positive_dict=False, positive_code=False, method_max_iter=1000, ): """Solve a dictionary learning matrix factorization problem. Finds the best dictionary and the corresponding sparse code for approximating the data matrix X by solving:: (U^*, V^*) = argmin 0.5 || X - U V ||_Fro^2 + alpha * || U ||_1,1 (U,V) with || V_k ||_2 = 1 for all 0 <= k < n_components where V is the dictionary and U is the sparse code. ||.||_Fro stands for the Frobenius norm and ||.||_1,1 stands for the entry-wise matrix norm which is the sum of the absolute values of all the entries in the matrix. Read more in the :ref:`User Guide `. Parameters ---------- X : array-like of shape (n_samples, n_features) Data matrix. n_components : int Number of dictionary atoms to extract. alpha : int or float Sparsity controlling parameter. max_iter : int, default=100 Maximum number of iterations to perform. tol : float, default=1e-8 Tolerance for the stopping condition. method : {'lars', 'cd'}, default='lars' The method used: * `'lars'`: uses the least angle regression method to solve the lasso problem (`linear_model.lars_path`); * `'cd'`: uses the coordinate descent method to compute the Lasso solution (`linear_model.Lasso`). Lars will be faster if the estimated components are sparse. n_jobs : int, default=None Number of parallel jobs to run. ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. ``-1`` means using all processors. See :term:`Glossary ` for more details. dict_init : ndarray of shape (n_components, n_features), default=None Initial value for the dictionary for warm restart scenarios. Only used if `code_init` and `dict_init` are not None. code_init : ndarray of shape (n_samples, n_components), default=None Initial value for the sparse code for warm restart scenarios. Only used if `code_init` and `dict_init` are not None. callback : callable, default=None Callable that gets invoked every five iterations. verbose : bool, default=False To control the verbosity of the procedure. random_state : int, RandomState instance or None, default=None Used for randomly initializing the dictionary. Pass an int for reproducible results across multiple function calls. See :term:`Glossary `. return_n_iter : bool, default=False Whether or not to return the number of iterations. positive_dict : bool, default=False Whether to enforce positivity when finding the dictionary. .. versionadded:: 0.20 positive_code : bool, default=False Whether to enforce positivity when finding the code. .. versionadded:: 0.20 method_max_iter : int, default=1000 Maximum number of iterations to perform. .. versionadded:: 0.22 Returns ------- code : ndarray of shape (n_samples, n_components) The sparse code factor in the matrix factorization. dictionary : ndarray of shape (n_components, n_features), The dictionary factor in the matrix factorization. errors : array Vector of errors at each iteration. n_iter : int Number of iterations run. Returned only if `return_n_iter` is set to True. See Also -------- dict_learning_online : Solve a dictionary learning matrix factorization problem online. DictionaryLearning : Find a dictionary that sparsely encodes data. MiniBatchDictionaryLearning : A faster, less accurate version of the dictionary learning algorithm. SparsePCA : Sparse Principal Components Analysis. MiniBatchSparsePCA : Mini-batch Sparse Principal Components Analysis. """ estimator = DictionaryLearning( n_components=n_components, alpha=alpha, max_iter=max_iter, tol=tol, fit_algorithm=method, n_jobs=n_jobs, dict_init=dict_init, callback=callback, code_init=code_init, verbose=verbose, random_state=random_state, positive_code=positive_code, positive_dict=positive_dict, transform_max_iter=method_max_iter, ) code = estimator.fit_transform(X) if return_n_iter: return ( code, estimator.components_, estimator.error_, estimator.n_iter_, ) return code, estimator.components_, estimator.error_ class _BaseSparseCoding(ClassNamePrefixFeaturesOutMixin, TransformerMixin): """Base class from SparseCoder and DictionaryLearning algorithms.""" def __init__( self, transform_algorithm, transform_n_nonzero_coefs, transform_alpha, split_sign, n_jobs, positive_code, transform_max_iter, ): self.transform_algorithm = transform_algorithm self.transform_n_nonzero_coefs = transform_n_nonzero_coefs self.transform_alpha = transform_alpha self.transform_max_iter = transform_max_iter self.split_sign = split_sign self.n_jobs = n_jobs self.positive_code = positive_code def _transform(self, X, dictionary): """Private method allowing to accommodate both DictionaryLearning and SparseCoder.""" X = self._validate_data(X, reset=False) if hasattr(self, "alpha") and self.transform_alpha is None: transform_alpha = self.alpha else: transform_alpha = self.transform_alpha code = sparse_encode( X, dictionary, algorithm=self.transform_algorithm, n_nonzero_coefs=self.transform_n_nonzero_coefs, alpha=transform_alpha, max_iter=self.transform_max_iter, n_jobs=self.n_jobs, positive=self.positive_code, ) if self.split_sign: # feature vector is split into a positive and negative side n_samples, n_features = code.shape split_code = np.empty((n_samples, 2 * n_features)) split_code[:, :n_features] = np.maximum(code, 0) split_code[:, n_features:] = -np.minimum(code, 0) code = split_code return code def transform(self, X): """Encode the data as a sparse combination of the dictionary atoms. Coding method is determined by the object parameter `transform_algorithm`. Parameters ---------- X : ndarray of shape (n_samples, n_features) Test data to be transformed, must have the same number of features as the data used to train the model. Returns ------- X_new : ndarray of shape (n_samples, n_components) Transformed data. """ check_is_fitted(self) return self._transform(X, self.components_) class SparseCoder(_BaseSparseCoding, BaseEstimator): """Sparse coding. Finds a sparse representation of data against a fixed, precomputed dictionary. Each row of the result is the solution to a sparse coding problem. The goal is to find a sparse array `code` such that:: X ~= code * dictionary Read more in the :ref:`User Guide `. Parameters ---------- dictionary : ndarray of shape (n_components, n_features) The dictionary atoms used for sparse coding. Lines are assumed to be normalized to unit norm. transform_algorithm : {'lasso_lars', 'lasso_cd', 'lars', 'omp', \ 'threshold'}, default='omp' Algorithm used to transform the data: - `'lars'`: uses the least angle regression method (`linear_model.lars_path`); - `'lasso_lars'`: uses Lars to compute the Lasso solution; - `'lasso_cd'`: uses the coordinate descent method to compute the Lasso solution (linear_model.Lasso). `'lasso_lars'` will be faster if the estimated components are sparse; - `'omp'`: uses orthogonal matching pursuit to estimate the sparse solution; - `'threshold'`: squashes to zero all coefficients less than alpha from the projection ``dictionary * X'``. transform_n_nonzero_coefs : int, default=None Number of nonzero coefficients to target in each column of the solution. This is only used by `algorithm='lars'` and `algorithm='omp'` and is overridden by `alpha` in the `omp` case. If `None`, then `transform_n_nonzero_coefs=int(n_features / 10)`. transform_alpha : float, default=None If `algorithm='lasso_lars'` or `algorithm='lasso_cd'`, `alpha` is the penalty applied to the L1 norm. If `algorithm='threshold'`, `alpha` is the absolute value of the threshold below which coefficients will be squashed to zero. If `algorithm='omp'`, `alpha` is the tolerance parameter: the value of the reconstruction error targeted. In this case, it overrides `n_nonzero_coefs`. If `None`, default to 1. split_sign : bool, default=False Whether to split the sparse feature vector into the concatenation of its negative part and its positive part. This can improve the performance of downstream classifiers. n_jobs : int, default=None Number of parallel jobs to run. ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. ``-1`` means using all processors. See :term:`Glossary ` for more details. positive_code : bool, default=False Whether to enforce positivity when finding the code. .. versionadded:: 0.20 transform_max_iter : int, default=1000 Maximum number of iterations to perform if `algorithm='lasso_cd'` or `lasso_lars`. .. versionadded:: 0.22 Attributes ---------- n_components_ : int Number of atoms. n_features_in_ : int Number of features seen during :term:`fit`. .. versionadded:: 0.24 feature_names_in_ : ndarray of shape (`n_features_in_`,) Names of features seen during :term:`fit`. Defined only when `X` has feature names that are all strings. .. versionadded:: 1.0 See Also -------- DictionaryLearning : Find a dictionary that sparsely encodes data. MiniBatchDictionaryLearning : A faster, less accurate, version of the dictionary learning algorithm. MiniBatchSparsePCA : Mini-batch Sparse Principal Components Analysis. SparsePCA : Sparse Principal Components Analysis. sparse_encode : Sparse coding where each row of the result is the solution to a sparse coding problem. Examples -------- >>> import numpy as np >>> from sklearn.decomposition import SparseCoder >>> X = np.array([[-1, -1, -1], [0, 0, 3]]) >>> dictionary = np.array( ... [[0, 1, 0], ... [-1, -1, 2], ... [1, 1, 1], ... [0, 1, 1], ... [0, 2, 1]], ... dtype=np.float64 ... ) >>> coder = SparseCoder( ... dictionary=dictionary, transform_algorithm='lasso_lars', ... transform_alpha=1e-10, ... ) >>> coder.transform(X) array([[ 0., 0., -1., 0., 0.], [ 0., 1., 1., 0., 0.]]) """ _required_parameters = ["dictionary"] def __init__( self, dictionary, *, transform_algorithm="omp", transform_n_nonzero_coefs=None, transform_alpha=None, split_sign=False, n_jobs=None, positive_code=False, transform_max_iter=1000, ): super().__init__( transform_algorithm, transform_n_nonzero_coefs, transform_alpha, split_sign, n_jobs, positive_code, transform_max_iter, ) self.dictionary = dictionary def fit(self, X, y=None): """Do nothing and return the estimator unchanged. This method is just there to implement the usual API and hence work in pipelines. Parameters ---------- X : Ignored Not used, present for API consistency by convention. y : Ignored Not used, present for API consistency by convention. Returns ------- self : object Returns the instance itself. """ return self def transform(self, X, y=None): """Encode the data as a sparse combination of the dictionary atoms. Coding method is determined by the object parameter `transform_algorithm`. Parameters ---------- X : ndarray of shape (n_samples, n_features) Training vector, where `n_samples` is the number of samples and `n_features` is the number of features. y : Ignored Not used, present for API consistency by convention. Returns ------- X_new : ndarray of shape (n_samples, n_components) Transformed data. """ return super()._transform(X, self.dictionary) def _more_tags(self): return { "requires_fit": False, "preserves_dtype": [np.float64, np.float32], } @property def n_components_(self): """Number of atoms.""" return self.dictionary.shape[0] @property def n_features_in_(self): """Number of features seen during `fit`.""" return self.dictionary.shape[1] @property def _n_features_out(self): """Number of transformed output features.""" return self.n_components_ class DictionaryLearning(_BaseSparseCoding, BaseEstimator): """Dictionary learning. Finds a dictionary (a set of atoms) that performs well at sparsely encoding the fitted data. Solves the optimization problem:: (U^*,V^*) = argmin 0.5 || X - U V ||_Fro^2 + alpha * || U ||_1,1 (U,V) with || V_k ||_2 <= 1 for all 0 <= k < n_components ||.||_Fro stands for the Frobenius norm and ||.||_1,1 stands for the entry-wise matrix norm which is the sum of the absolute values of all the entries in the matrix. Read more in the :ref:`User Guide `. Parameters ---------- n_components : int, default=None Number of dictionary elements to extract. If None, then ``n_components`` is set to ``n_features``. alpha : float, default=1.0 Sparsity controlling parameter. max_iter : int, default=1000 Maximum number of iterations to perform. tol : float, default=1e-8 Tolerance for numerical error. fit_algorithm : {'lars', 'cd'}, default='lars' * `'lars'`: uses the least angle regression method to solve the lasso problem (:func:`~sklearn.linear_model.lars_path`); * `'cd'`: uses the coordinate descent method to compute the Lasso solution (:class:`~sklearn.linear_model.Lasso`). Lars will be faster if the estimated components are sparse. .. versionadded:: 0.17 *cd* coordinate descent method to improve speed. transform_algorithm : {'lasso_lars', 'lasso_cd', 'lars', 'omp', \ 'threshold'}, default='omp' Algorithm used to transform the data: - `'lars'`: uses the least angle regression method (:func:`~sklearn.linear_model.lars_path`); - `'lasso_lars'`: uses Lars to compute the Lasso solution. - `'lasso_cd'`: uses the coordinate descent method to compute the Lasso solution (:class:`~sklearn.linear_model.Lasso`). `'lasso_lars'` will be faster if the estimated components are sparse. - `'omp'`: uses orthogonal matching pursuit to estimate the sparse solution. - `'threshold'`: squashes to zero all coefficients less than alpha from the projection ``dictionary * X'``. .. versionadded:: 0.17 *lasso_cd* coordinate descent method to improve speed. transform_n_nonzero_coefs : int, default=None Number of nonzero coefficients to target in each column of the solution. This is only used by `algorithm='lars'` and `algorithm='omp'`. If `None`, then `transform_n_nonzero_coefs=int(n_features / 10)`. transform_alpha : float, default=None If `algorithm='lasso_lars'` or `algorithm='lasso_cd'`, `alpha` is the penalty applied to the L1 norm. If `algorithm='threshold'`, `alpha` is the absolute value of the threshold below which coefficients will be squashed to zero. If `None`, defaults to `alpha`. .. versionchanged:: 1.2 When None, default value changed from 1.0 to `alpha`. n_jobs : int or None, default=None Number of parallel jobs to run. ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. ``-1`` means using all processors. See :term:`Glossary ` for more details. code_init : ndarray of shape (n_samples, n_components), default=None Initial value for the code, for warm restart. Only used if `code_init` and `dict_init` are not None. dict_init : ndarray of shape (n_components, n_features), default=None Initial values for the dictionary, for warm restart. Only used if `code_init` and `dict_init` are not None. callback : callable, default=None Callable that gets invoked every five iterations. .. versionadded:: 1.3 verbose : bool, default=False To control the verbosity of the procedure. split_sign : bool, default=False Whether to split the sparse feature vector into the concatenation of its negative part and its positive part. This can improve the performance of downstream classifiers. random_state : int, RandomState instance or None, default=None Used for initializing the dictionary when ``dict_init`` is not specified, randomly shuffling the data when ``shuffle`` is set to ``True``, and updating the dictionary. Pass an int for reproducible results across multiple function calls. See :term:`Glossary `. positive_code : bool, default=False Whether to enforce positivity when finding the code. .. versionadded:: 0.20 positive_dict : bool, default=False Whether to enforce positivity when finding the dictionary. .. versionadded:: 0.20 transform_max_iter : int, default=1000 Maximum number of iterations to perform if `algorithm='lasso_cd'` or `'lasso_lars'`. .. versionadded:: 0.22 Attributes ---------- components_ : ndarray of shape (n_components, n_features) dictionary atoms extracted from the data error_ : array vector of errors at each iteration n_features_in_ : int Number of features seen during :term:`fit`. .. versionadded:: 0.24 feature_names_in_ : ndarray of shape (`n_features_in_`,) Names of features seen during :term:`fit`. Defined only when `X` has feature names that are all strings. .. versionadded:: 1.0 n_iter_ : int Number of iterations run. See Also -------- MiniBatchDictionaryLearning: A faster, less accurate, version of the dictionary learning algorithm. MiniBatchSparsePCA : Mini-batch Sparse Principal Components Analysis. SparseCoder : Find a sparse representation of data from a fixed, precomputed dictionary. SparsePCA : Sparse Principal Components Analysis. References ---------- J. Mairal, F. Bach, J. Ponce, G. Sapiro, 2009: Online dictionary learning for sparse coding (https://www.di.ens.fr/sierra/pdfs/icml09.pdf) Examples -------- >>> import numpy as np >>> from sklearn.datasets import make_sparse_coded_signal >>> from sklearn.decomposition import DictionaryLearning >>> X, dictionary, code = make_sparse_coded_signal( ... n_samples=30, n_components=15, n_features=20, n_nonzero_coefs=10, ... random_state=42, ... ) >>> dict_learner = DictionaryLearning( ... n_components=15, transform_algorithm='lasso_lars', transform_alpha=0.1, ... random_state=42, ... ) >>> X_transformed = dict_learner.fit(X).transform(X) We can check the level of sparsity of `X_transformed`: >>> np.mean(X_transformed == 0) 0.52... We can compare the average squared euclidean norm of the reconstruction error of the sparse coded signal relative to the squared euclidean norm of the original signal: >>> X_hat = X_transformed @ dict_learner.components_ >>> np.mean(np.sum((X_hat - X) ** 2, axis=1) / np.sum(X ** 2, axis=1)) 0.05... """ _parameter_constraints: dict = { "n_components": [Interval(Integral, 1, None, closed="left"), None], "alpha": [Interval(Real, 0, None, closed="left")], "max_iter": [Interval(Integral, 0, None, closed="left")], "tol": [Interval(Real, 0, None, closed="left")], "fit_algorithm": [StrOptions({"lars", "cd"})], "transform_algorithm": [ StrOptions({"lasso_lars", "lasso_cd", "lars", "omp", "threshold"}) ], "transform_n_nonzero_coefs": [Interval(Integral, 1, None, closed="left"), None], "transform_alpha": [Interval(Real, 0, None, closed="left"), None], "n_jobs": [Integral, None], "code_init": [np.ndarray, None], "dict_init": [np.ndarray, None], "callback": [callable, None], "verbose": ["verbose"], "split_sign": ["boolean"], "random_state": ["random_state"], "positive_code": ["boolean"], "positive_dict": ["boolean"], "transform_max_iter": [Interval(Integral, 0, None, closed="left")], } def __init__( self, n_components=None, *, alpha=1, max_iter=1000, tol=1e-8, fit_algorithm="lars", transform_algorithm="omp", transform_n_nonzero_coefs=None, transform_alpha=None, n_jobs=None, code_init=None, dict_init=None, callback=None, verbose=False, split_sign=False, random_state=None, positive_code=False, positive_dict=False, transform_max_iter=1000, ): super().__init__( transform_algorithm, transform_n_nonzero_coefs, transform_alpha, split_sign, n_jobs, positive_code, transform_max_iter, ) self.n_components = n_components self.alpha = alpha self.max_iter = max_iter self.tol = tol self.fit_algorithm = fit_algorithm self.code_init = code_init self.dict_init = dict_init self.callback = callback self.verbose = verbose self.random_state = random_state self.positive_dict = positive_dict def fit(self, X, y=None): """Fit the model from data in X. Parameters ---------- X : array-like of shape (n_samples, n_features) Training vector, where `n_samples` is the number of samples and `n_features` is the number of features. y : Ignored Not used, present for API consistency by convention. Returns ------- self : object Returns the instance itself. """ self.fit_transform(X) return self @_fit_context(prefer_skip_nested_validation=True) def fit_transform(self, X, y=None): """Fit the model from data in X and return the transformed data. Parameters ---------- X : array-like of shape (n_samples, n_features) Training vector, where `n_samples` is the number of samples and `n_features` is the number of features. y : Ignored Not used, present for API consistency by convention. Returns ------- V : ndarray of shape (n_samples, n_components) Transformed data. """ _check_positive_coding(method=self.fit_algorithm, positive=self.positive_code) method = "lasso_" + self.fit_algorithm random_state = check_random_state(self.random_state) X = self._validate_data(X) if self.n_components is None: n_components = X.shape[1] else: n_components = self.n_components V, U, E, self.n_iter_ = _dict_learning( X, n_components, alpha=self.alpha, tol=self.tol, max_iter=self.max_iter, method=method, method_max_iter=self.transform_max_iter, n_jobs=self.n_jobs, code_init=self.code_init, dict_init=self.dict_init, callback=self.callback, verbose=self.verbose, random_state=random_state, return_n_iter=True, positive_dict=self.positive_dict, positive_code=self.positive_code, ) self.components_ = U self.error_ = E return V @property def _n_features_out(self): """Number of transformed output features.""" return self.components_.shape[0] def _more_tags(self): return { "preserves_dtype": [np.float64, np.float32], } class MiniBatchDictionaryLearning(_BaseSparseCoding, BaseEstimator): """Mini-batch dictionary learning. Finds a dictionary (a set of atoms) that performs well at sparsely encoding the fitted data. Solves the optimization problem:: (U^*,V^*) = argmin 0.5 || X - U V ||_Fro^2 + alpha * || U ||_1,1 (U,V) with || V_k ||_2 <= 1 for all 0 <= k < n_components ||.||_Fro stands for the Frobenius norm and ||.||_1,1 stands for the entry-wise matrix norm which is the sum of the absolute values of all the entries in the matrix. Read more in the :ref:`User Guide `. Parameters ---------- n_components : int, default=None Number of dictionary elements to extract. alpha : float, default=1 Sparsity controlling parameter. n_iter : int, default=1000 Total number of iterations over data batches to perform. .. deprecated:: 1.1 ``n_iter`` is deprecated in 1.1 and will be removed in 1.4. Use ``max_iter`` instead. max_iter : int, default=None Maximum number of iterations over the complete dataset before stopping independently of any early stopping criterion heuristics. If ``max_iter`` is not None, ``n_iter`` is ignored. .. versionadded:: 1.1 fit_algorithm : {'lars', 'cd'}, default='lars' The algorithm used: - `'lars'`: uses the least angle regression method to solve the lasso problem (`linear_model.lars_path`) - `'cd'`: uses the coordinate descent method to compute the Lasso solution (`linear_model.Lasso`). Lars will be faster if the estimated components are sparse. n_jobs : int, default=None Number of parallel jobs to run. ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. ``-1`` means using all processors. See :term:`Glossary ` for more details. batch_size : int, default=256 Number of samples in each mini-batch. .. versionchanged:: 1.3 The default value of `batch_size` changed from 3 to 256 in version 1.3. shuffle : bool, default=True Whether to shuffle the samples before forming batches. dict_init : ndarray of shape (n_components, n_features), default=None Initial value of the dictionary for warm restart scenarios. transform_algorithm : {'lasso_lars', 'lasso_cd', 'lars', 'omp', \ 'threshold'}, default='omp' Algorithm used to transform the data: - `'lars'`: uses the least angle regression method (`linear_model.lars_path`); - `'lasso_lars'`: uses Lars to compute the Lasso solution. - `'lasso_cd'`: uses the coordinate descent method to compute the Lasso solution (`linear_model.Lasso`). `'lasso_lars'` will be faster if the estimated components are sparse. - `'omp'`: uses orthogonal matching pursuit to estimate the sparse solution. - `'threshold'`: squashes to zero all coefficients less than alpha from the projection ``dictionary * X'``. transform_n_nonzero_coefs : int, default=None Number of nonzero coefficients to target in each column of the solution. This is only used by `algorithm='lars'` and `algorithm='omp'`. If `None`, then `transform_n_nonzero_coefs=int(n_features / 10)`. transform_alpha : float, default=None If `algorithm='lasso_lars'` or `algorithm='lasso_cd'`, `alpha` is the penalty applied to the L1 norm. If `algorithm='threshold'`, `alpha` is the absolute value of the threshold below which coefficients will be squashed to zero. If `None`, defaults to `alpha`. .. versionchanged:: 1.2 When None, default value changed from 1.0 to `alpha`. verbose : bool or int, default=False To control the verbosity of the procedure. split_sign : bool, default=False Whether to split the sparse feature vector into the concatenation of its negative part and its positive part. This can improve the performance of downstream classifiers. random_state : int, RandomState instance or None, default=None Used for initializing the dictionary when ``dict_init`` is not specified, randomly shuffling the data when ``shuffle`` is set to ``True``, and updating the dictionary. Pass an int for reproducible results across multiple function calls. See :term:`Glossary `. positive_code : bool, default=False Whether to enforce positivity when finding the code. .. versionadded:: 0.20 positive_dict : bool, default=False Whether to enforce positivity when finding the dictionary. .. versionadded:: 0.20 transform_max_iter : int, default=1000 Maximum number of iterations to perform if `algorithm='lasso_cd'` or `'lasso_lars'`. .. versionadded:: 0.22 callback : callable, default=None A callable that gets invoked at the end of each iteration. .. versionadded:: 1.1 tol : float, default=1e-3 Control early stopping based on the norm of the differences in the dictionary between 2 steps. Used only if `max_iter` is not None. To disable early stopping based on changes in the dictionary, set `tol` to 0.0. .. versionadded:: 1.1 max_no_improvement : int, default=10 Control early stopping based on the consecutive number of mini batches that does not yield an improvement on the smoothed cost function. Used only if `max_iter` is not None. To disable convergence detection based on cost function, set `max_no_improvement` to None. .. versionadded:: 1.1 Attributes ---------- components_ : ndarray of shape (n_components, n_features) Components extracted from the data. n_features_in_ : int Number of features seen during :term:`fit`. .. versionadded:: 0.24 feature_names_in_ : ndarray of shape (`n_features_in_`,) Names of features seen during :term:`fit`. Defined only when `X` has feature names that are all strings. .. versionadded:: 1.0 n_iter_ : int Number of iterations over the full dataset. n_steps_ : int Number of mini-batches processed. .. versionadded:: 1.1 See Also -------- DictionaryLearning : Find a dictionary that sparsely encodes data. MiniBatchSparsePCA : Mini-batch Sparse Principal Components Analysis. SparseCoder : Find a sparse representation of data from a fixed, precomputed dictionary. SparsePCA : Sparse Principal Components Analysis. References ---------- J. Mairal, F. Bach, J. Ponce, G. Sapiro, 2009: Online dictionary learning for sparse coding (https://www.di.ens.fr/sierra/pdfs/icml09.pdf) Examples -------- >>> import numpy as np >>> from sklearn.datasets import make_sparse_coded_signal >>> from sklearn.decomposition import MiniBatchDictionaryLearning >>> X, dictionary, code = make_sparse_coded_signal( ... n_samples=30, n_components=15, n_features=20, n_nonzero_coefs=10, ... random_state=42) >>> dict_learner = MiniBatchDictionaryLearning( ... n_components=15, batch_size=3, transform_algorithm='lasso_lars', ... transform_alpha=0.1, max_iter=20, random_state=42) >>> X_transformed = dict_learner.fit_transform(X) We can check the level of sparsity of `X_transformed`: >>> np.mean(X_transformed == 0) > 0.5 True We can compare the average squared euclidean norm of the reconstruction error of the sparse coded signal relative to the squared euclidean norm of the original signal: >>> X_hat = X_transformed @ dict_learner.components_ >>> np.mean(np.sum((X_hat - X) ** 2, axis=1) / np.sum(X ** 2, axis=1)) 0.052... """ _parameter_constraints: dict = { "n_components": [Interval(Integral, 1, None, closed="left"), None], "alpha": [Interval(Real, 0, None, closed="left")], "n_iter": [ Interval(Integral, 0, None, closed="left"), Hidden(StrOptions({"deprecated"})), ], "max_iter": [Interval(Integral, 0, None, closed="left"), None], "fit_algorithm": [StrOptions({"cd", "lars"})], "n_jobs": [None, Integral], "batch_size": [Interval(Integral, 1, None, closed="left")], "shuffle": ["boolean"], "dict_init": [None, np.ndarray], "transform_algorithm": [ StrOptions({"lasso_lars", "lasso_cd", "lars", "omp", "threshold"}) ], "transform_n_nonzero_coefs": [Interval(Integral, 1, None, closed="left"), None], "transform_alpha": [Interval(Real, 0, None, closed="left"), None], "verbose": ["verbose"], "split_sign": ["boolean"], "random_state": ["random_state"], "positive_code": ["boolean"], "positive_dict": ["boolean"], "transform_max_iter": [Interval(Integral, 0, None, closed="left")], "callback": [None, callable], "tol": [Interval(Real, 0, None, closed="left")], "max_no_improvement": [Interval(Integral, 0, None, closed="left"), None], } def __init__( self, n_components=None, *, alpha=1, n_iter="deprecated", max_iter=None, fit_algorithm="lars", n_jobs=None, batch_size=256, shuffle=True, dict_init=None, transform_algorithm="omp", transform_n_nonzero_coefs=None, transform_alpha=None, verbose=False, split_sign=False, random_state=None, positive_code=False, positive_dict=False, transform_max_iter=1000, callback=None, tol=1e-3, max_no_improvement=10, ): super().__init__( transform_algorithm, transform_n_nonzero_coefs, transform_alpha, split_sign, n_jobs, positive_code, transform_max_iter, ) self.n_components = n_components self.alpha = alpha self.n_iter = n_iter self.max_iter = max_iter self.fit_algorithm = fit_algorithm self.dict_init = dict_init self.verbose = verbose self.shuffle = shuffle self.batch_size = batch_size self.split_sign = split_sign self.random_state = random_state self.positive_dict = positive_dict self.callback = callback self.max_no_improvement = max_no_improvement self.tol = tol def _check_params(self, X): # n_components self._n_components = self.n_components if self._n_components is None: self._n_components = X.shape[1] # fit_algorithm _check_positive_coding(self.fit_algorithm, self.positive_code) self._fit_algorithm = "lasso_" + self.fit_algorithm # batch_size self._batch_size = min(self.batch_size, X.shape[0]) def _initialize_dict(self, X, random_state): """Initialization of the dictionary.""" if self.dict_init is not None: dictionary = self.dict_init else: # Init V with SVD of X _, S, dictionary = randomized_svd( X, self._n_components, random_state=random_state ) dictionary = S[:, np.newaxis] * dictionary if self._n_components <= len(dictionary): dictionary = dictionary[: self._n_components, :] else: dictionary = np.concatenate( ( dictionary, np.zeros( (self._n_components - len(dictionary), dictionary.shape[1]), dtype=dictionary.dtype, ), ) ) dictionary = check_array(dictionary, order="F", dtype=X.dtype, copy=False) dictionary = np.require(dictionary, requirements="W") return dictionary def _update_inner_stats(self, X, code, batch_size, step): """Update the inner stats inplace.""" if step < batch_size - 1: theta = (step + 1) * batch_size else: theta = batch_size**2 + step + 1 - batch_size beta = (theta + 1 - batch_size) / (theta + 1) self._A *= beta self._A += code.T @ code / batch_size self._B *= beta self._B += X.T @ code / batch_size def _minibatch_step(self, X, dictionary, random_state, step): """Perform the update on the dictionary for one minibatch.""" batch_size = X.shape[0] # Compute code for this batch code = _sparse_encode( X, dictionary, algorithm=self._fit_algorithm, alpha=self.alpha, n_jobs=self.n_jobs, positive=self.positive_code, max_iter=self.transform_max_iter, verbose=self.verbose, ) batch_cost = ( 0.5 * ((X - code @ dictionary) ** 2).sum() + self.alpha * np.sum(np.abs(code)) ) / batch_size # Update inner stats self._update_inner_stats(X, code, batch_size, step) # Update dictionary _update_dict( dictionary, X, code, self._A, self._B, verbose=self.verbose, random_state=random_state, positive=self.positive_dict, ) return batch_cost def _check_convergence( self, X, batch_cost, new_dict, old_dict, n_samples, step, n_steps ): """Helper function to encapsulate the early stopping logic. Early stopping is based on two factors: - A small change of the dictionary between two minibatch updates. This is controlled by the tol parameter. - No more improvement on a smoothed estimate of the objective function for a a certain number of consecutive minibatch updates. This is controlled by the max_no_improvement parameter. """ batch_size = X.shape[0] # counts steps starting from 1 for user friendly verbose mode. step = step + 1 # Ignore 100 first steps or 1 epoch to avoid initializing the ewa_cost with a # too bad value if step <= min(100, n_samples / batch_size): if self.verbose: print(f"Minibatch step {step}/{n_steps}: mean batch cost: {batch_cost}") return False # Compute an Exponentially Weighted Average of the cost function to # monitor the convergence while discarding minibatch-local stochastic # variability: https://en.wikipedia.org/wiki/Moving_average if self._ewa_cost is None: self._ewa_cost = batch_cost else: alpha = batch_size / (n_samples + 1) alpha = min(alpha, 1) self._ewa_cost = self._ewa_cost * (1 - alpha) + batch_cost * alpha if self.verbose: print( f"Minibatch step {step}/{n_steps}: mean batch cost: " f"{batch_cost}, ewa cost: {self._ewa_cost}" ) # Early stopping based on change of dictionary dict_diff = linalg.norm(new_dict - old_dict) / self._n_components if self.tol > 0 and dict_diff <= self.tol: if self.verbose: print(f"Converged (small dictionary change) at step {step}/{n_steps}") return True # Early stopping heuristic due to lack of improvement on smoothed # cost function if self._ewa_cost_min is None or self._ewa_cost < self._ewa_cost_min: self._no_improvement = 0 self._ewa_cost_min = self._ewa_cost else: self._no_improvement += 1 if ( self.max_no_improvement is not None and self._no_improvement >= self.max_no_improvement ): if self.verbose: print( "Converged (lack of improvement in objective function) " f"at step {step}/{n_steps}" ) return True return False @_fit_context(prefer_skip_nested_validation=True) def fit(self, X, y=None): """Fit the model from data in X. Parameters ---------- X : array-like of shape (n_samples, n_features) Training vector, where `n_samples` is the number of samples and `n_features` is the number of features. y : Ignored Not used, present for API consistency by convention. Returns ------- self : object Returns the instance itself. """ X = self._validate_data( X, dtype=[np.float64, np.float32], order="C", copy=False ) self._check_params(X) if self.n_iter != "deprecated": warnings.warn( ( "'n_iter' is deprecated in version 1.1 and will be removed " "in version 1.4. Use 'max_iter' and let 'n_iter' to its default " "value instead. 'n_iter' is also ignored if 'max_iter' is " "specified." ), FutureWarning, ) n_iter = self.n_iter self._random_state = check_random_state(self.random_state) dictionary = self._initialize_dict(X, self._random_state) old_dict = dictionary.copy() if self.shuffle: X_train = X.copy() self._random_state.shuffle(X_train) else: X_train = X n_samples, n_features = X_train.shape if self.verbose: print("[dict_learning]") # Inner stats self._A = np.zeros( (self._n_components, self._n_components), dtype=X_train.dtype ) self._B = np.zeros((n_features, self._n_components), dtype=X_train.dtype) if self.max_iter is not None: # Attributes to monitor the convergence self._ewa_cost = None self._ewa_cost_min = None self._no_improvement = 0 batches = gen_batches(n_samples, self._batch_size) batches = itertools.cycle(batches) n_steps_per_iter = int(np.ceil(n_samples / self._batch_size)) n_steps = self.max_iter * n_steps_per_iter i = -1 # to allow max_iter = 0 for i, batch in zip(range(n_steps), batches): X_batch = X_train[batch] batch_cost = self._minibatch_step( X_batch, dictionary, self._random_state, i ) if self._check_convergence( X_batch, batch_cost, dictionary, old_dict, n_samples, i, n_steps ): break # XXX callback param added for backward compat in #18975 but a common # unified callback API should be preferred if self.callback is not None: self.callback(locals()) old_dict[:] = dictionary self.n_steps_ = i + 1 self.n_iter_ = np.ceil(self.n_steps_ / n_steps_per_iter) else: # TODO remove this branch in 1.4 n_iter = 1000 if self.n_iter == "deprecated" else self.n_iter batches = gen_batches(n_samples, self._batch_size) batches = itertools.cycle(batches) for i, batch in zip(range(n_iter), batches): self._minibatch_step(X_train[batch], dictionary, self._random_state, i) trigger_verbose = self.verbose and i % ceil(100.0 / self.verbose) == 0 if self.verbose > 10 or trigger_verbose: print(f"{i} batches processed.") if self.callback is not None: self.callback(locals()) self.n_steps_ = n_iter self.n_iter_ = np.ceil(n_iter / int(np.ceil(n_samples / self._batch_size))) self.components_ = dictionary return self @_fit_context(prefer_skip_nested_validation=True) def partial_fit(self, X, y=None): """Update the model using the data in X as a mini-batch. Parameters ---------- X : array-like of shape (n_samples, n_features) Training vector, where `n_samples` is the number of samples and `n_features` is the number of features. y : Ignored Not used, present for API consistency by convention. Returns ------- self : object Return the instance itself. """ has_components = hasattr(self, "components_") X = self._validate_data( X, dtype=[np.float64, np.float32], order="C", reset=not has_components ) if not has_components: # This instance has not been fitted yet (fit or partial_fit) self._check_params(X) self._random_state = check_random_state(self.random_state) dictionary = self._initialize_dict(X, self._random_state) self.n_steps_ = 0 self._A = np.zeros((self._n_components, self._n_components), dtype=X.dtype) self._B = np.zeros((X.shape[1], self._n_components), dtype=X.dtype) else: dictionary = self.components_ self._minibatch_step(X, dictionary, self._random_state, self.n_steps_) self.components_ = dictionary self.n_steps_ += 1 return self @property def _n_features_out(self): """Number of transformed output features.""" return self.components_.shape[0] def _more_tags(self): return { "preserves_dtype": [np.float64, np.float32], }