From 07980bc1410c90ae826b492df31248238f7c4e13 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Fri, 6 Feb 2026 12:36:54 -0500 Subject: [PATCH 1/2] feat: skip update_stretch() when rho = 0, reducing program to regular NMF --- src/diffpy/snmf/snmf_class.py | 60 +++++++++++++++++++---------------- 1 file changed, 33 insertions(+), 27 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 13076d3..4fd4310 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -69,9 +69,9 @@ def __init__( init_weights=None, init_components=None, init_stretch=None, - max_iter=500, + max_iter=100, min_iter=20, - tol=5e-7, + tol=1e-6, n_components=None, random_state=None, show_plots=False, @@ -125,8 +125,8 @@ def __init__( n_components is not None and init_weights is not None ): raise ValueError( - "Conflicting source for n_components. Must provide either init_weights or n_components " - "directly, but not both." + "Conflicting or missing source for n_components. Must provide either init_weights " + "or n_components directly, but not both." ) # Initialize weights and determine number of components @@ -181,6 +181,7 @@ def fit(self, rho=0, eta=0, reset=True): rho : float Optional Default = 0 The stretching factor that influences the decomposition. Zero corresponds to no stretching present. Relatively insensitive and typically adjusted in powers of 10. + If equal to zero, program acts like standard NMF. eta : int Optional Default = 0 The sparsity factor that influences the decomposition. Should be set to zero for non-sparse data such as PDF. Can be used to improve results for sparse data such @@ -200,6 +201,10 @@ def fit(self, rho=0, eta=0, reset=True): self.rho = rho self.eta = eta + # If rho = 0, set stretching matrix to identity + if self.rho == 0: + self.stretch_ = np.ones_like(self.stretch_) + # Set up residual matrix, objective function, and history self.residuals = self.get_residual_matrix() self.objective_function = self.get_objective_function() @@ -386,30 +391,31 @@ def outer_loop(self): ): break - self.update_stretch() - self.residuals = self.get_residual_matrix() - self.objective_function = self.get_objective_function() - print( - f"Objective function after update_stretch: {self.objective_function:.5e}" - ) - self._objective_history.append(self.objective_function) - self.objective_difference = ( - self._objective_history[-2] - self._objective_history[-1] - ) - if self.objective_function < self.best_objective: - self.best_objective = self.objective_function - self.best_matrices = [ - self.components_.copy(), - self.weights_.copy(), - self.stretch_.copy(), - ] - if self.plotter is not None: - self.plotter.update( - components=self.components_, - weights=self.weights_, - stretch=self.stretch_, - update_tag="stretch", + if not self.rho == 0: # Don't update stretch if rho = 0 + self.update_stretch() + self.residuals = self.get_residual_matrix() + self.objective_function = self.get_objective_function() + print( + f"Objective function after update_stretch: {self.objective_function:.5e}" + ) + self._objective_history.append(self.objective_function) + self.objective_difference = ( + self._objective_history[-2] - self._objective_history[-1] ) + if self.objective_function < self.best_objective: + self.best_objective = self.objective_function + self.best_matrices = [ + self.components_.copy(), + self.weights_.copy(), + self.stretch_.copy(), + ] + if self.plotter is not None: + self.plotter.update( + components=self.components_, + weights=self.weights_, + stretch=self.stretch_, + update_tag="stretch", + ) def get_residual_matrix(self, components=None, weights=None, stretch=None): """ From 6facf5f0891427a298de271e0b17962875923970 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Fri, 6 Feb 2026 21:00:10 -0500 Subject: [PATCH 2/2] chore: fix line lengths not updated by black --- src/diffpy/snmf/snmf_class.py | 312 ++++++++++++++++++++-------------- 1 file changed, 184 insertions(+), 128 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 4fd4310..a4d3a5a 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -7,25 +7,27 @@ class SNMFOptimizer: - """An implementation of stretched NMF (sNMF), including sparse stretched NMF. + """An implementation of stretched NMF (sNMF), + including sparse stretched NMF. - Instantiating the SNMFOptimizer class prepares initial guesses and sets up the - optimization. It can then be run using fit(). + Instantiating the SNMFOptimizer class prepares initial guesses and + sets up the optimization. It can then be run using fit(). The results matrices can be accessed as instance attributes of the class (components_, weights_, and stretch_). For more information on sNMF, please reference: - Gu, R., Rakita, Y., Lan, L. et al. Stretched non-negative matrix factorization. - npj Comput Mater 10, 193 (2024). https://doi.org/10.1038/s41524-024-01377-5 + Gu, R., Rakita, Y., Lan, L. et al. + Stretched non-negative matrix factorization. + npj Compu Mater 10, 193 (2024). https://doi.org/10.1038/s41524-024-01377-5 Attributes ---------- source_matrix : ndarray - The original, unmodified data to be decomposed and later, compared against. - Shape is (length_of_signal, number_of_signals). + The original, unmodified data to be decomposed and later, + compared against. Shape is (length_of_signal, number_of_signals). stretch_ : ndarray - The best guess (or while running, the current guess) for the stretching - factor matrix. + The best guess (or while running, the current guess) for the + stretching factor matrix. components_ : ndarray The best guess (or while running, the current guess) for the matrix of component intensities. @@ -33,34 +35,38 @@ class SNMFOptimizer: The best guess (or while running, the current guess) for the matrix of component weights. rho : float - The stretching factor that influences the decomposition. Zero corresponds to no - stretching present. Relatively insensitive and typically adjusted in powers of 10. + The stretching factor that influences the decomposition. Zero + corresponds to no stretching present. Relatively insensitive and + typically adjusted in powers of 10. eta : float - The sparsity factor that influences the decomposition. Should be set to zero for - non-sparse data such as PDF. Can be used to improve results for sparse data such - as XRD, but due to instability, should be used only after first selecting the - best value for rho. Suggested adjustment is by powers of 2. + The sparsity factor that influences the decomposition. Should be set + to zero for non-sparse data such as PDF. Can be used to improve + results for sparse data such as XRD, but due to instability, should + be used only after first selecting the best value for rho. Suggested + adjustment is by powers of 2. max_iter : int - The maximum number of times to update each of stretch, components, and weights before stopping - the optimization. + The maximum number of times to update each of stretch, components, + and weights before stopping the optimization. min_iter : int - The minimum number of times to update each of stretch, components, and weights before terminating - the optimization due to low/no improvement. + The minimum number of times to update each of stretch, components, and + weights before terminating the optimization due to low/no improvement. tol : float - The convergence threshold. This is the minimum fractional improvement in the - objective function to allow without terminating the optimization. + The convergence threshold. This is the minimum fractional improvement + in the objective function to allow without terminating the + optimization. n_components : int - The number of components to extract from source_matrix. Must be provided when and only when - init_weights is not provided. + The number of components to extract from source_matrix. Must be + provided when and only when init_weights is not provided. random_state : int - The seed for the initial guesses at the matrices (stretch, components, and weights) created by - the decomposition. + The seed for the initial guesses at the matrices (stretch, components, + and weights) created by the decomposition. num_updates : int - The total number of times that any of (stretch, components, and weights) have had their values changed. - If not terminated by other means, this value is used to stop when reaching max_iter. + The total number of times that any of (stretch, components, + and weights) have had their values changed. If not terminated by + other means, this value is used to stop when reaching max_iter. objective_difference : float - The change in the objective function value since the last update. A negative value - means that the result improved. + The change in the objective function value since the last update. + A negative value means that the result improved. """ def __init__( @@ -69,9 +75,9 @@ def __init__( init_weights=None, init_components=None, init_stretch=None, - max_iter=100, + max_iter=500, min_iter=20, - tol=1e-6, + tol=5e-7, n_components=None, random_state=None, show_plots=False, @@ -81,31 +87,37 @@ def __init__( Parameters ---------- source_matrix : ndarray - The data to be decomposed. Shape is (length_of_signal, number_of_conditions). - init_weights : ndarray Optional Default = rng.beta(a=2.0, b=2.0, size=(n_components, n_signals)) - The initial guesses for the component weights at each stretching condition. - Shape is (number_of_components, number_of_signals) Must provide exactly one - of this or n_components. - init_components : ndarray Optional Default = rng.random((self.signal_length, self.n_components)) + The data to be decomposed. Shape is (length_of_signal, + number_of_conditions). + init_weights : ndarray Optional Default = rng.beta(a=2.0, b=2.0, + size=(n_components, n_signals)) + The initial guesses for the component weights at each stretching + condition. Shape is (number_of_components, number_of_signals) + Must provide exactly one of this or n_components. + init_components : ndarray Optional Default = rng.random( + (self.signal_length, self.n_components)) The initial guesses for the intensities of each component per - row/sample/angle. Shape is (length_of_signal, number_of_components). - init_stretch : ndarray Optional Default = np.ones((self.n_components, self.n_signals)) + self._rng.normal( + row/sample/angle. Shape is (length_of_signal, number_of_components) + init_stretch : ndarray Optional Default = np.ones((self.n_components, + self.n_signals)) + self._rng.normal( 0, 1e-3, size=(self.n_components, self.n_signals) - The initial guesses for the stretching factor for each component, at each - condition (for each signal). Shape is (number_of_components, number_of_signals). + The initial guesses for the stretching factor for each component, + at each condition (for each signal). + Shape is (number_of_components, number_of_signals). max_iter : int Optional Default = 500 - The maximum number of times to update each of A, X, and Y before stopping - the optimization. + The maximum number of times to update each of A, X, and Y before + stopping the optimization. tol : float Optional Default = 5e-7 - The convergence threshold. This is the minimum fractional improvement in the - objective function to allow without terminating the optimization. Note that - a minimum of 20 updates are run before this parameter is checked. + The convergence threshold. This is the minimum fractional + improvement in the objective function to allow without terminating + the optimization. Note that a minimum of 20 updates are run before + this parameter is checked. n_components : int Optional Default = None - The number of components to extract from source_matrix. Must be provided when and only when - Y0 is not provided. + The number of components to extract from source_matrix. Must be + provided when and only when Y0 is not provided. random_state : int Optional Default = None - The seed for the initial guesses at the matrices (A, X, and Y) created by - the decomposition. + The seed for the initial guesses at the matrices (A, X, and Y) + created by the decomposition. show_plots : boolean Optional Default = False Enables plotting at each step of the decomposition. """ @@ -125,8 +137,8 @@ def __init__( n_components is not None and init_weights is not None ): raise ValueError( - "Conflicting or missing source for n_components. Must provide either init_weights " - "or n_components directly, but not both." + "Conflicting source for n_components. Must provide either " + "init_weights or n_components directly, but not both." ) # Initialize weights and determine number of components @@ -166,7 +178,7 @@ def __init__( self.init_weights = self.weights_.copy() self.init_stretch = self.stretch_.copy() - # Second-order spline: Tridiagonal (-2 on diagonal, 1 on sub/superdiagonals) + # Second-order spline: Tridiagonal (-2 on diag, 1 on sub/superdiags) self._spline_smooth_operator = 0.25 * diags( [1, -2, 1], offsets=[0, 1, 2], @@ -174,23 +186,26 @@ def __init__( ) def fit(self, rho=0, eta=0, reset=True): - """Run the sNMF optimization with the given parameters, using the setup from __init__. + """Run the sNMF optimization with the given parameters, + using the setup from __init__. Parameters ---------- rho : float Optional Default = 0 - The stretching factor that influences the decomposition. Zero corresponds to no - stretching present. Relatively insensitive and typically adjusted in powers of 10. - If equal to zero, program acts like standard NMF. + The stretching factor that influences the decomposition. Zero + corresponds to no stretching present. Relatively insensitive and + typically adjusted in powers of 10. eta : int Optional Default = 0 - The sparsity factor that influences the decomposition. Should be set to zero for - non-sparse data such as PDF. Can be used to improve results for sparse data such - as XRD, but due to instability, should be used only after first selecting the + The sparsity factor that influences the decomposition. Should be + set to zero for non-sparse data such as PDF. Can be used to + improve results for sparse data such as XRD, but due to + instability, should be used only after first selecting the best value for rho. Suggested adjustment is by powers of 2. reset : boolean Optional Default = True - Whether to return to the initial set of components_, weights_, and stretch_ before - running the optimization. When set to False, sequential calls to fit() will use the - output of the previous fit() as their input. + Whether to return to the initial set of components_, weights_, + and stretch_ before running the optimization. When set to False, + sequential calls to fit() will use the output of the previous + fit() as their input. """ if reset: @@ -201,10 +216,6 @@ def fit(self, rho=0, eta=0, reset=True): self.rho = rho self.eta = eta - # If rho = 0, set stretching matrix to identity - if self.rho == 0: - self.stretch_ = np.ones_like(self.stretch_) - # Set up residual matrix, objective function, and history self.residuals = self.get_residual_matrix() self.objective_function = self.get_objective_function() @@ -235,7 +246,8 @@ def fit(self, rho=0, eta=0, reset=True): ) # Square root penalty print( f"Start, Objective function: {self.objective_function:.5e}" - f", Obj - reg/sparse: {self.objective_function - regularization_term - sparsity_term:.5e}" + f", Obj - reg/sparse: {self.objective_function - + regularization_term - sparsity_term:.5e}" ) # Main optimization loop @@ -256,11 +268,13 @@ def fit(self, rho=0, eta=0, reset=True): ) # Square root penalty print( f"Obj fun: {self.objective_function:.5e}, " - f"Obj - reg/sparse: {self.objective_function - regularization_term - sparsity_term:.5e}, " + f"Obj - reg/sparse: {self.objective_function + - regularization_term + - sparsity_term:.5e}, " f"Iter: {self.outiter}" ) - # Convergence check: Stop if diffun is small and at least min_iter iterations have passed + # Conv. check: Stop if diffun small and min_iter iterations passed print( "Checking if ", self.objective_difference, @@ -289,7 +303,7 @@ def normalize_results(self): stretch_row_max = np.max(self.stretch_, axis=1, keepdims=True) self.stretch_ = self.stretch_ / stretch_row_max - # effectively just re-running with component updates only vs normalized weights/stretch + # just re-running with component updates only vs norm. weights/stretch self._grad_components = np.zeros_like( self.components_ ) # Gradient of X (zeros for now) @@ -309,7 +323,8 @@ def normalize_results(self): self.residuals = self.get_residual_matrix() self.objective_function = self.get_objective_function() print( - f"Objective function after normalize_components: {self.objective_function:.5e}" + f"Objective function after normalize_components: " + f"{self.objective_function:.5e}" ) self._objective_history.append(self.objective_function) self.objective_difference = ( @@ -336,7 +351,8 @@ def outer_loop(self): self.residuals = self.get_residual_matrix() self.objective_function = self.get_objective_function() print( - f"Objective function after update_components: {self.objective_function:.5e}" + f"Objective function after update_components: " + f"{self.objective_function:.5e}" ) self._objective_history.append(self.objective_function) self.objective_difference = ( @@ -361,7 +377,8 @@ def outer_loop(self): self.residuals = self.get_residual_matrix() self.objective_function = self.get_objective_function() print( - f"Objective function after update_weights: {self.objective_function:.5e}" + f"Objective function after update_weights: " + f"{self.objective_function:.5e}" ) self._objective_history.append(self.objective_function) self.objective_difference = ( @@ -391,35 +408,36 @@ def outer_loop(self): ): break - if not self.rho == 0: # Don't update stretch if rho = 0 - self.update_stretch() - self.residuals = self.get_residual_matrix() - self.objective_function = self.get_objective_function() - print( - f"Objective function after update_stretch: {self.objective_function:.5e}" - ) - self._objective_history.append(self.objective_function) - self.objective_difference = ( - self._objective_history[-2] - self._objective_history[-1] + self.update_stretch() + self.residuals = self.get_residual_matrix() + self.objective_function = self.get_objective_function() + print( + f"Objective function after update_stretch: " + f"{self.objective_function:.5e}" + ) + self._objective_history.append(self.objective_function) + self.objective_difference = ( + self._objective_history[-2] - self._objective_history[-1] + ) + if self.objective_function < self.best_objective: + self.best_objective = self.objective_function + self.best_matrices = [ + self.components_.copy(), + self.weights_.copy(), + self.stretch_.copy(), + ] + if self.plotter is not None: + self.plotter.update( + components=self.components_, + weights=self.weights_, + stretch=self.stretch_, + update_tag="stretch", ) - if self.objective_function < self.best_objective: - self.best_objective = self.objective_function - self.best_matrices = [ - self.components_.copy(), - self.weights_.copy(), - self.stretch_.copy(), - ] - if self.plotter is not None: - self.plotter.update( - components=self.components_, - weights=self.weights_, - stretch=self.stretch_, - update_tag="stretch", - ) def get_residual_matrix(self, components=None, weights=None, stretch=None): """ - Return the residuals (difference) between the source matrix and its reconstruction. + Return the residuals (difference) between the source matrix and + its reconstruction. Parameters ---------- @@ -474,10 +492,12 @@ def compute_stretched_components( self, components=None, weights=None, stretch=None ): """ - Interpolates each component along its sample axis according to per-(component, signal) - stretch factors, then applies per-(component, signal) weights. Also computes the - first and second derivatives with respect to stretch. Left and right, respectively, - refer to the sample prior to and subsequent to the interpolated sample's position. + Interpolates each component along its sample axis according to + per-(component, signal) stretch factors, then applies + per-(component, signal) weights. Also computes the first and second + derivatives with respect to stretch. Left and right, respectively, + refer to the sample prior to and subsequent to the interpolated + sample's position. Inputs ------ @@ -490,11 +510,14 @@ def compute_stretched_components( Outputs ------- - stretched_components : array, shape (signal_len, n_components * n_signals) + stretched_components : array, + shape (signal_len, n_components * n_signals) Interpolated and weighted components. - d_stretched_components : array, shape (signal_len, n_components * n_signals) + d_stretched_components : array, shape + (signal_len, n_components * n_signals) First derivatives with respect to stretch. - dd_stretched_components : array, shape (signal_len, n_components * n_signals) + dd_stretched_components : array, shape + (signal_len, n_components * n_signals) Second derivatives with respect to stretch. """ @@ -516,18 +539,18 @@ def compute_stretched_components( stretch = np.clip(stretch, eps, None) stretch_inv = 1.0 / stretch - # Apply stretching to the original sample indices, represented as a "time-stretch" + # Apply stretching to original sample indices as a "time-stretch" t = ( np.arange(signal_len, dtype=float)[:, None, None] * stretch_inv[None, :, :] ) # has shape (signal_len, n_components, n_signals) - # For each stretched coordinate, find its prior integer (original) index and their difference + # For each str coordinate, find prior integer index and their diff i0 = np.floor(t).astype(np.int64) # prior original index alpha = t - i0.astype(float) # fractional distance between left/right - # Clip indices to valid range (0, signal_len - 1) to maintain original size + # Clip to valid range (0, signal_len - 1) to maintain original size max_idx = signal_len - 1 i0 = np.clip(i0, 0, max_idx) i1 = np.clip(i0 + 1, 0, max_idx) @@ -535,7 +558,7 @@ def compute_stretched_components( # Gather sample values comps_3d = components[ :, :, None - ] # expand components by a dimension for broadcasting across n_signals + ] # expand components by a dimension for broadcast across n_signals c0 = np.take_along_axis(comps_3d, i0, axis=0) # left sample values c1 = np.take_along_axis(comps_3d, i1, axis=0) # right sample values @@ -555,7 +578,7 @@ def compute_stretched_components( d_weighted = d_unweighted * weights[None, :, :] dd_weighted = dd_unweighted * weights[None, :, :] - # Flatten back to expected shape (signal_len, n_components * n_signals) + # Flatten to expected shape (signal_len, n_components * n_signals) return ( interp_weighted.reshape(signal_len, n_components * n_signals), d_weighted.reshape(signal_len, n_components * n_signals), @@ -646,7 +669,8 @@ def solve_quadratic_program(self, t, m): Parameters: - t: (N, k) ndarray - - source_matrix_col: (N,) column of source_matrix for the corresponding m + - source_matrix_col: (N,) column of source_matrix for the + corresponding m Returns: - y: (k,) optimal solution @@ -686,7 +710,8 @@ def solve_quadratic_program(self, t, m): def update_components(self): """ - Updates `components` using gradient-based optimization with adaptive step size. + Updates `components` using gradient-based optimization with + adaptive step size. """ # Compute stretched components using the interpolation function stretched_components, _, _ = ( @@ -767,15 +792,17 @@ def update_components(self): def update_weights(self): """ - Updates weights by building the stretched component matrix `stretched_comps` with np.interp - and solving a quadratic program for each signal. + Updates weights by building the stretched component matrix + `stretched_comps` with np.interp and solving a quadratic program + for each signal. """ sample_indices = np.arange(self.signal_length) for signal in range(self.n_signals): # Stretch factors for this signal across components: this_stretch = self.stretch_[:, signal] - # Build stretched_comps[:, k] by interpolating component at frac. pos. index / this_stretch[comp] + # Build stretched_comps[:, k] by interpolating component + # at frac. pos. index / this_stretch[comp] stretched_comps = np.empty( (self.signal_length, self.n_components), dtype=self.components_.dtype, @@ -832,13 +859,17 @@ def regularize_function(self, stretch=None): def update_stretch(self): """ - Updates stretching matrix using constrained optimization (equivalent to fmincon in MATLAB). + Updates stretching matrix using constrained optimization + (equivalent to fmincon in MATLAB). """ - # Flatten stretch for compatibility with the optimizer (since SciPy expects 1D input) + # Flatten stretch for compatibility with the optimizer + # (since SciPy expects 1D input) stretch_flat_initial = self.stretch_.flatten() # Define the optimization function + cache = {"x": None, "fun": None, "grad": None} + def objective(stretch_vec): stretch_matrix = stretch_vec.reshape( self.stretch_.shape @@ -847,6 +878,26 @@ def objective(stretch_vec): gra = gra.flatten() return fun, gra + def fun_only(stretch_vec): + if cache["x"] is None or not np.array_equal( + stretch_vec, cache["x"] + ): + fun, grad = objective(stretch_vec) + cache["x"] = stretch_vec.copy() + cache["fun"] = fun + cache["grad"] = grad + return cache["fun"] + + def jac_only(stretch_vec): + if cache["x"] is None or not np.array_equal( + stretch_vec, cache["x"] + ): + fun, grad = objective(stretch_vec) + cache["x"] = stretch_vec.copy() + cache["fun"] = fun + cache["grad"] = grad + return cache["grad"] + # Optimization constraints: lower bound 0.1, no upper bound bounds = [ (0.1, None) @@ -854,10 +905,10 @@ def objective(stretch_vec): # Solve optimization problem (equivalent to fmincon) result = minimize( - fun=lambda stretch_vec: objective(stretch_vec)[0], + fun=fun_only, x0=stretch_flat_initial, method="trust-constr", # Substitute for 'trust-region-reflective' - jac=lambda stretch_vec: objective(stretch_vec)[1], # Gradient + jac=jac_only, # Gradient bounds=bounds, ) @@ -869,7 +920,8 @@ def _compute_objective_function( components, residuals, stretch, rho, eta, spline_smooth_operator ): r""" - Computes the objective function used in stretched non-negative matrix factorization. + Computes the objective function used in stretched + non-negative matrix factorization. Parameters ---------- @@ -878,13 +930,15 @@ def _compute_objective_function( residuals : ndarray Difference between reconstructed and observed data. stretch : ndarray - Stretching factors :math:`A` applied to each component across samples. + Stretching factors :math:`A` applied to each component + across samples. rho : float Regularization parameter enforcing smooth variation in :math:`A`. eta : float Sparsity-promoting regularization parameter applied to :math:`X`. spline_smooth_operator : ndarray - Linear operator :math:`L` penalizing non-smooth changes in :math:`A`. + Linear operator :math:`L` penalizing non-smooth changes + in :math:`A`. Returns ------- @@ -902,13 +956,14 @@ def _compute_objective_function( + \tfrac{\rho}{2} \lVert L A \rVert_F^2 + \eta \sum_{i,j} \sqrt{X_{ij}} \,, - where :math:`Z` is the data matrix, :math:`Y` contains the non-negative - weights, :math:`S(A)` denotes the spline-interpolated stretching operator, - and :math:`\lVert \cdot \rVert_F` is the Frobenius norm. + where :math:`Z` is the data matrix, :math:`Y` contains + the non-negative weights, :math:`S(A)` denotes the spline-interpolated + stretching operator, and :math:`\lVert \cdot \rVert_F` is the + Frobenius norm. Special cases ------------- - - :math:`\rho = 0` — no smoothness regularization on stretching factors. + - :math:`\rho = 0` — no smoothness regularization on stretch factors - :math:`\eta = 0` — no sparsity promotion on components. - :math:`\rho = \eta = 0` — reduces to the classical NMF least-squares objective :math:`\tfrac{1}{2} \lVert Z - YX \rVert_F^2`. @@ -926,7 +981,8 @@ def _compute_objective_function( def cubic_largest_real_root(p, q): """ - Solves x^3 + p*x + q = 0 element-wise for matrices, returning the largest real root. + Solves x^3 + p*x + q = 0 element-wise for matrices, returning the + largest real root. """ # Handle special case where q == 0 y = np.where(