Positive matrix factorization of large real-time atmospheric mass spectrometry datasets using error-weighted randomized hierarchical alternating least squares

Weighted positive matrix factorization (PMF) has been used by scientists to find small sets of underlying factors in environmental data. However, as the size of the data has grown, increasing computational costs have made it impractical to use traditional methods for this factorization. In this paper, we present a new external weighting method to dramatically decrease computational costs for these traditional algorithms. The external weighting scheme, along with the randomized hierarchical alternating least squares (RHALS) algorithm, was applied to the Southern Oxidant and Aerosol Study (SOAS 2013) dataset of gaseous highly oxidized multifunctional molecules (HOMs). The modified RHALS algorithm successfully reproduced six previously identified interpretable factors, with the total computation time of the nonoptimized code showing potential improvements of the order of 1 to 2 orders of magnitude compared to competing algorithms. We also investigate rotational ambiguity in the solution and present a simple “pulling” method to rotate a set of factors. This method is shown to find alternative solutions and, in some cases, lower the weighted residual error of the algorithm.

1 Introduction

1.1 Problem statement

Low-rank matrix factorization has been widely used in data science to explain underlying factors in large datasets (Xie et al., 1998; Kim and Hopke, 2007; Wei et al., 2016). The process considers a data matrix, A, of size m×n, which is decomposed into two smaller matrices, W of size m×k and H of size k×n, where $k\ll \min (m,n)$ and A≈WH. Traditionally, principal component analysis (PCA) and singular value decomposition (SVD) have been used to find these factors (Kumar, 2017; Wei et al., 2016). PCA finds the eigenvectors of the covariance matrix A^TA, which are called the principal components, representing the directions of maximum variance in the data. The vectors are ordered by how much variance they explain, and only the most important vectors are kept, which are identified as the underlying factors in the data. Closely related, the SVD finds the factorization A=USV^T, where U and V contain the left and right singular vectors, respectively, and S is a diagonal matrix containing the singular values of A in decreasing order. Thus, to find a low-rank approximation of A, one could keep only the k most significant singular values and vectors to form the truncated SVD, ${\mathbf{\text{U}}}_k{\mathbf{\text{S}}}_k{\mathbf{\text{V}}}_k^T$. Mathematically, the truncated SVD is the most optimal rank k factorization of A for minimizing squared error (Eckart and Young, 1936). However, the SVD is not appropriate for all factorization problems for several reasons (Paatero and Tapper, 1994):

1. The SVD produces factors with negative values. For some factor analysis problems, such as finding chemical sources for air pollution data, SVD results can be difficult to interpret as chemical concentrations can only be nonnegative.

2. The SVD produces orthogonal factors. Many factor analysis problems are not constrained by requirements of orthogonality between factors.

3. The SVD is not fit to solve the following weighted least squares problem. Suppose that accompanying the dataset A is an equally sized (m×n) matrix Σ with Σ_ij=σ_ij representing the uncertainty of the measurement for A_ij. If rank(Σ)>1, the SVD cannot be scaled to find a solution minimizing the weighted residual error, which is defined as ${\sum }_{i=\mathrm{1}}^m{\sum }_{j=\mathrm{1}}^n\frac{({\mathbf{\text{A}}}_{ij}-{\sum }_{l=\mathrm{1}}^k{\mathbf{\text{W}}}_{il}{\mathbf{\text{H}}}_{lj}{)}^{\mathrm{2}}}{{\mathit{\sigma }}_{ij}^{\mathrm{2}}}$ (Paatero and Tapper, 1994).

Positive matrix factorization (PMF) was introduced by Paatero and Tapper (1993) to address these concerns. Weighted PMF attempts to find two factor matrices, W and H, by minimizing the equation

$$\begin{array}{}\text{(1)}& \left|\right|(\mathbf{\text{A}}-\mathbf{\text{W}}\mathbf{\text{H}})\oslash \mathbf{\Sigma }|{|}_F^{\mathrm{2}}\text{ with }\mathbf{\text{W}}\ge \mathrm{0}\text{,}\mathbf{\text{H}}\ge \mathrm{0}.\end{array}$$

In Eq. (1), ⊘ represents elementwise division, the norm $\left|\right|\cdot |{|}_F$ is the Frobenius norm, and all elements of W and H are constrained to be nonnegative. Further, we note that for consistency with nomenclature in the literature related to the use of this algorithm for factorization of real-time atmospheric mass spectrometry datasets, we refer to this approach as “positive” matrix factorization (i.e., PMF) while recognizing that a more precise name would be nonnegative matrix factorization (NMF).

Traditional factor analysis methods are known to be computationally expensive. Steps to speed up factor analysis have been explored, such as randomization and the use of graphical processing units (GPUs) (Halko et al., 2011; Tan et al., 2018). Developing efficient algorithms is especially critical in atmospheric mass spectrometry, as improvements in instrumentation and increases in the duration of their use in field campaigns have led to intractably large datasets. Currently, analysis of these datasets requires sacrificing data resolution or extensive manual preprocessing to operate within existing PMF software tools, and full analysis can routinely take days or weeks of computation time (Hopke et al., 2023). As a result, a variety of approaches have emerged for efficient source apportionment of atmospheric mass spectrometry data. Algorithms to solve the nonconvex optimization posed by PMF range from gradient descent and block coordinate descent to projected gradient methods (Guo et al., 2024). Attempts at using supervised ensemble machine learning approaches have been shown to be capable of replicating results from traditional (unsupervised) factorization methods while reducing computation time (Zhang et al., 2025). Recently, Erichson et al. (2018) applied randomization to PMF and introduced a new method, randomized hierarchical alternating least squares (RHALS), to solve the unweighted PMF problem. In this paper, we test the application of RHALS to atmospheric concentration data that contain uncertainties. Accounting for these uncertainties as regression weights, we introduce a method of externally weighting and unweighting the data, which to our knowledge is novel in its application to RHALS. We consider the accuracy and the reduced computational costs compared to other PMF algorithms commonly used in the field of atmospheric science.

1.2 Background

1.2.1 PMF2 and Paatero

The first widely accepted algorithm for PMF was derived in Paatero (1997) using the Gauss–Newton method. This algorithm, called PMF2, is commonly used with environmental data (Kim and Hopke, 2007; Ulbrich et al., 2009; Massoli et al., 2018). Paatero (1997) defines an enhanced objective function and attempts to find factor matrices W and H that minimize the cost function Q in

$$\begin{array}{}\text{(2)}& \begin{array}{rl}Q& =\sum _{i=\mathrm{1}}^m\sum _{j=\mathrm{1}}^n{\displaystyle \frac{({\mathbf{\text{A}}}_{ij}-{\sum }_{l=\mathrm{1}}^k{\mathbf{\text{W}}}_{il}{\mathbf{\text{H}}}_{lj}{)}^{\mathrm{2}}}{{\mathit{\sigma }}_{ij}^{\mathrm{2}}}}\\ & -\mathit{\alpha }\sum _{i=\mathrm{1}}^m\sum _{l=\mathrm{1}}^k\log \left({\mathbf{\text{W}}}_{il}\right)-\mathit{\beta }\sum _{l=\mathrm{1}}^k\sum _{j=\mathrm{1}}^n\log \left({\mathbf{\text{H}}}_{lj}\right)\\ & +\mathit{\gamma }\sum _{i=\mathrm{1}}^m\sum _{l=\mathrm{1}}^k{\mathbf{\text{W}}}_{il}^{\mathrm{2}}+\mathit{\delta }\sum _{l=\mathrm{1}}^k\sum _{j=\mathrm{1}}^n{\mathbf{\text{H}}}_{lj}^{\mathrm{2}}.\end{array}\end{array}$$

In Eq. (2), logarithmic penalty terms are added to penalize factor values that become too close to zero (and therefore potentially negative), and L2 regularization is added to smooth out the factors and avoid overfitting (Paatero, 1997). Here, α and β control the strength of the penalty terms, while γ and δ control the strength of L2 regularization. To date, Paatero's exact algorithmic approach to solving Eq. (2) remains unpublished. However, pseudo-code for using the Gauss–Newton method to solve Eq. (2) is detailed in Lu and Wu (2004).

1.2.2 Multiplicative update

An alternative method for PMF was developed in Lee and Seung (1999). The multiplicative update (MU) method utilizes a special case of gradient descent where the learning rates are chosen to avoid subtraction in the gradient (Gillis, 2020). A multiplicative update estimate of a parameter θ (either W or H) is found by updates of the form (Gillis, 2020)

$$\begin{array}{}\text{(3)}& \mathit{\theta }=\mathit{\theta }\odot \left({\mathrm{\nabla }}_{\mathit{\theta }}^{-}Q\right(\mathit{\theta })\oslash {\mathrm{\nabla }}_{\mathit{\theta }}^{+}Q(\mathit{\theta }\left)\right),\end{array}$$

where Q(θ) is a cost function to be minimized, ${\mathrm{\nabla }}_{\mathit{\theta }}^{-}$ consists of the negative terms of the gradient of the cost function, ${\mathrm{\nabla }}_{\mathit{\theta }}^{+}$ consists of the positive terms of the gradient, and ⊙ denotes elementwise multiplication. Here, θ is initialized with all positive entries as the MU cannot update an entry θ_ij if it is equal to zero (Gillis, 2020). As the data matrix A, the uncertainties (σ_ij), and the factor matrices are all nonnegative at each step of this algorithm, the factor matrices in the subsequent step are guaranteed to be positive because Eq. (3) only deals with the multiplication and division of positive numbers.

It is possible to perform PMF using other forms of gradient descent – for example, the projected gradient method (PGM) sets the step size to the inverse of the maximum eigenvalue of the Hessian of the cost function and may lead to faster convergence than MU (Gillis, 2020). However, we choose to only test MU due to its widespread use and flexibility (Gillis, 2020).

1.2.3 Alternating least squares and hierarchical ALS

Alternating least squares (ALS) methods solve for the factor matrices W and H by iteratively updating each matrix until convergence is reached (Cichocki et al., 2009). The cost function Q, containing $\left|\right|(\mathbf{\text{A}}-\mathbf{\text{W}}\mathbf{\text{H}})\oslash \mathbf{\Sigma }|{|}_F^{\mathrm{2}}$, is minimized by setting the partial derivatives $\frac{\partial Q}{\partial \mathbf{\text{W}}}$ and $\frac{\partial Q}{\partial \mathbf{\text{H}}}$ to zero and solving for W and H. To satisfy the positivity constraint, negative elements in the factors are set to zero.

Nonnegative ALS has no theoretical convergence guarantee and, in some problems, may fail to converge to a feasible solution (Gillis, 2020). For this reason, the alternating nonnegative least squares (ANLS) method and the alternating direction method of multipliers (ADMM) are interesting alternatives. In ANLS, indices of an “active set” are set to zero, and the rest are updated via an unconstrained optimization (Kim and Park, 2011). The active set is then updated to contain the indices with the new negative factor elements. In ADMM, an auxiliary factor matrix Y is formed, and an additional term is added to the cost function, which penalizes the distance between the target factor matrix (W or H) and Y (Gillis, 2020). Both of these methods may lead to faster and better convergence than nonnegative ALS (Gillis, 2020). However, we find that the simple nonnegative ALS almost always converges to feasible solutions for our dataset, and we do not explore these alternative methods.

In recent years, hierarchical alternating least squares (HALS) has become increasingly popular as an efficient method for PMF (Cichocki and Phan, 2009). Instead of minimizing with respect to the entire factor matrices W and H, HALS minimizes the cost function with respect to one block, or an outer product, of individual factors at a time. The main component of the cost function is redefined as $Q_j=\left|\right|({\mathbf{\text{R}}}_j-{\mathbf{\text{W}}}_{(:,j)}{\mathbf{\text{H}}}_{(j,:)})\oslash \mathbf{\Sigma }|{|}_F^{\mathrm{2}}$, where ${\mathbf{\text{W}}}_{(:,j)}$ and ${\mathbf{\text{H}}}_{(j,:)}$ are the jth factors. Q_j can be minimized for each factor j by setting the partial derivatives $\frac{\partial Q_j}{\partial {\mathbf{\text{W}}}_{(:,j)}}$ and $\frac{\partial Q_j}{\partial {\mathbf{\text{H}}}_{(j,:)}}$ to zero and solving for ${\mathbf{\text{W}}}_{(:,j)}$ and ${\mathbf{\text{H}}}_{(j,:)}$.

The derivation of the ALS update rules is detailed in Appendix A, while the derivation of HALS is detailed in Sect. 2.

1.3 Random projections

To reduce the computational costs of a matrix factorization algorithm for large datasets, randomization methods have been used as a dimension reduction technique (Erichson et al., 2018; Halko et al., 2011; Kaloorazi and Chen, 2019). Below, we present a brief overview of the theory and results laid out in Halko et al. (2011).

When performing randomization techniques, we hope that much of the relevant information about the column space of the data matrix A can be stored in a much smaller subset of vectors that we can sample. This is only true if the “effective rank” of A is low (A only has a few nonnegligible singular values), but that is generally assumed to be the case in any PMF problem (Erichson et al., 2018). Mathematically, we seek the approximation

$$\begin{array}{}\text{(4)}& \mathbf{\text{A}}\approx \mathbf{\text{P}}{\mathbf{\text{P}}}^T\mathbf{\text{A}},\end{array}$$

where the relatively small number of columns in the matrix P are orthonormal and form an approximate basis of A. Choosing the columns of P to be the left singular vectors of A would minimize the L2 error: if the first k singular vectors were chosen, the error term $\left|\right|\mathbf{\text{A}}-\mathbf{\text{P}}{\mathbf{\text{P}}}^T\mathbf{\text{A}}|{|}_{\mathrm{2}}={\mathit{\sigma }}_{k+\mathrm{1}}$, with σ_k+1 being the singular value of A with index (k+1) (Halko et al., 2011). However, random sampling from the column space of A can also produce a suitable basis.

The assumption that there are k underlying factors within A implies that the effective rank of A is k. It then appears reasonable to use k random samples from the column space of A to form a basis. However, with underlying uncertainties, we can write $\mathbf{\text{A}}=\mathbf{\text{B}}+\mathbf{\text{E}}$, where B is the rank k matrix spanned by the factors for which we wish to find a basis, and E is a perturbation matrix filled with the noise in A (Halko et al., 2011). Suppose we were to sample from the column space of A – that is, form the vector $\mathit{y}=\mathbf{\text{A}}\mathit{\omega }=\mathbf{\text{B}}\mathit{\omega }+\mathbf{\text{E}}\mathit{\omega }$. Each vector y is slightly pushed out of the column space of B by the term Eω. Thus, to increase the likelihood of spanning the column space of B, an additional p vectors are sampled from A. In practice, choosing p to be 10 or 20 is sufficient (Erichson et al., 2018).

To construct this low-rank approximation, k+p random normal samples, stored as columns of the matrix Ω (dimensions $n\times (k+p)$) of the column space of A, are taken and stored in Y:

$$\begin{array}{}\text{(5)}& \underset{m\times (k+p)}{\mathbf{\text{Y}}}=\underset{m\times n}{\mathbf{\text{A}}}\cdot \underset{n\times (k+p)}{\mathbf{\Omega }}.\end{array}$$

Next, the columns of Y are orthonormalized using a QR decomposition to form our projection matrix P. The algorithm can now be run on the lower dimensional matrix B=P^TA.

1.4 Nonuniqueness of solutions

Unlike the SVD, there is no guarantee of uniqueness for the factor matrices W and H in PMF. That is, the factorization A=WH can also be expressed as $\mathbf{\text{A}}=\mathbf{\text{W}}\mathbf{\text{T}}{\mathbf{\text{T}}}^{-\mathrm{1}}\mathbf{\text{H}}$, where T is a “rotational” matrix and $\widehat{\mathbf{\text{W}}}=\mathbf{\text{W}}\mathbf{\text{T}}$ and $\widehat{\mathbf{\text{H}}}={\mathbf{\text{T}}}^{-\mathrm{1}}\mathbf{\text{H}}$ are the new rotated factors. We note that T does not necessarily represent a true rotation in a mathematical form, which would require T to be orthogonal. To span the space of feasible solutions, previous approaches such as PMF2 have aimed to find new solutions $\widehat{\mathbf{\text{W}}}$ and $\widehat{\mathbf{\text{H}}}$ by varying T (Paatero, 1997) or varying initializations (Ulbrich et al., 2009).

The rotational matrix T is a k×k matrix where t_ii, a diagonal element of T, represents a scaling of the ith factor, and t_ij represents a rotation of the jth factor towards the ith factor in W and a rotation of the ith factor away from the jth factor in H. For example, consider the elementary rotation matrices (Paatero et al., 2002)

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\mathbf{\text{T}}}_E=\left(\begin{array}{ccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& r& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{1}\end{array}\right)\text{, }\\ & {\mathbf{\text{T}}}_E^{-\mathrm{1}}=\left(\begin{array}{ccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& -r& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{1}\end{array}\right).\end{array}\end{array}$$

All factors remain the same, except ${\widehat{\mathbf{\text{W}}}}_{(:,\mathrm{3})}={\mathbf{\text{W}}}_{(:,\mathrm{3})}+r{\mathbf{\text{W}}}_{(:,\mathrm{2})}$ and ${\widehat{\mathbf{\text{H}}}}_{(\mathrm{2},:)}={\mathbf{\text{H}}}_{(\mathrm{2},:)}-r{\mathbf{\text{H}}}_{(\mathrm{3},:)}$.

Regardless of whether r is positive or negative, values in either W or H will be pulled towards negative values. Thus if a large proportion of the factors are filled with zeros, there may be little to no pure rotations or rotations that do not change the residual $\left|\right|\mathbf{\text{A}}-\mathbf{\text{W}}\mathbf{\text{H}}|{|}_F^{\mathrm{2}}$. Algorithms such as those developed by Paatero, where a logarithmic penalty term is added to push factor values to be more positive, will have few zeros in the factor matrices and will thus have more rotational ambiguity (Paatero et al., 2005). However, RHALS enforces nonnegativity merely by setting negative elements to zero, and thus many values in the factor matrices may end up being zero. Thus for RHALS, a perfect rotational matrix T will almost certainly not exist, and only “approximate” rotations can be studied, in which the rotation will alter the value of the weighted error.

It is not feasible to span all possible variants that T can take. Thus, the problem is often simplified to considering only positive rotations (values of T greater than zero) and negative rotations (values of T less than zero). A rotational program in PMF2 called FPEAK uses the parameter ϕ to denote the rotation strength, with positive values leading to positive rotations in W (Paatero, 1997). Paatero further improved this method in the multilinear engine (ME) algorithm, where the strength of rotation is allowed to vary between factors (Paatero and Hopke, 2009). The pulling algorithm presented in Paatero and Hopke (2009) is a sophisticated rotational method; more rudimentary pulling methods that mimic varying the regularization of the factor matrices are presented in Paatero (1997) and Paatero et al. (2002). Recent attempts at controlling for rotational ambiguity have involved additional factorization of the time-series matrix W into a matrix incorporating shape regularization to reflect known diurnal patterns of factors and a diagonal scaling matrix (Nanra et al., 2024).

1.5 Scaling with uncertainties

Recall that to account for inaccuracies in real data, we measure the squared error Q of the algorithm by dividing the residual by the standard deviation of the uncertainty of each measurement:

$$\begin{array}{}\text{(7)}& Q=\sum _{i=\mathrm{1}}^m\sum _{j=\mathrm{1}}^n{\left({\displaystyle \frac{{\mathbf{\text{A}}}_{ij}-{\sum }_{l=\mathrm{1}}^k{\mathbf{\text{W}}}_{il}{\mathbf{\text{H}}}_{lj}}{{\mathit{\sigma }}_{ij}}}\right)}^{\mathrm{2}}.\end{array}$$

To account for these uncertainties, one could incorporate them into each update rule of the factor matrices, as is done in PMF2 (Paatero, 1997). In this paper, we refer to this as “internal weighting”. However, this is computationally expensive due to repeated elementwise operations with the uncertainty matrix Σ. Elementwise operations of large arrays are inefficient processes compared to other operations of the same computational complexity, such as matrix–vector multiplication, due to the large allocation of memory towards intermediary results (Jia et al., 2020). We thus introduce an alternate approach to weighted PMF where the data are prescaled by the uncertainty matrix, the unweighted algorithm is applied to the scaled data, and the converged factors are scaled by the uncertainties. This approach, which we refer to as “external weighting”, dramatically reduces computational costs and allows for dimensionality reduction as weights are not included in the update rules.

It is noted in Paatero and Tapper (1993) that if the weighting matrix is rank 1, that is if $\frac{\mathrm{1}}{{\mathit{\sigma }}_{ij}}={\mathbf{\text{B}}}_{(i,:)}{\mathbf{\text{C}}}_{(:,j)}$ for vectors of B and C, then an optimal scaling can be found. By forming the diagonal matrices D_L (with the elements of ${\mathbf{\text{B}}}_{(i,:)}$) and D_R (with the elements of ${\mathbf{\text{C}}}_{(:,j)}$), if $\mathbf{\text{A}}=\mathbf{\text{W}}\mathbf{\text{H}}+\mathbf{\text{E}}$, then ${\mathbf{\text{D}}}_L\mathbf{\text{A}}{\mathbf{\text{D}}}_R={\mathbf{\text{D}}}_L\mathbf{\text{W}}\mathbf{\text{H}}{\mathbf{\text{D}}}_R+{\mathbf{\text{D}}}_L\mathbf{\text{E}}{\mathbf{\text{D}}}_R=\widehat{\mathbf{\text{W}}}\widehat{\mathbf{\text{H}}}+\widehat{\mathbf{\text{E}}}$. Thus, by first finding $\widehat{\mathbf{\text{A}}}={\mathbf{\text{D}}}_L\mathbf{\text{A}}{\mathbf{\text{D}}}_R$ and then running a PMF algorithm without weights on $\widehat{\mathbf{\text{A}}}$, one can produce estimated factor matrices $\widehat{\mathbf{\text{W}}}$ and $\widehat{\mathbf{\text{H}}}$, where the unscaled estimates are $\mathbf{\text{W}}={\mathbf{\text{D}}}_L^{-\mathrm{1}}\widehat{\mathbf{\text{W}}}$ and $\mathbf{\text{H}}=\widehat{\mathbf{\text{H}}}{\mathbf{\text{D}}}_R^{-\mathrm{1}}$.

If the weighting matrix Σ is not rank 1 – which is likely for environmental data – and $\frac{\mathrm{1}}{{\mathit{\sigma }}_{i,j}}$ cannot be estimated as the outer product of two vectors, there is no scaling of the previous form that can be applied to the data matrix A (Paatero and Tapper, 1993). To address this, Paatero and Tapper (1993) presented a simple algorithm to find an approximate rank 1 factorization of Σ. We present a different method where the data are scaled by the full-rank matrix Σ and then unscaled after the algorithm is complete. This algorithm is described in detail in Sect. 2.3.

Expectation maximization

The expectation maximization (EM) approach was first designed for matrix factorization problems associated with missing entries (Zhang et al., 2006). Specifically, if A^o is the observed data and A^u is the unknown data within A, then the EM approach seeks to find factors W and H that satisfy the following (Zhang et al., 2006):

$$\begin{array}{}\text{(8)}& {\mathrm{argmax}}_{\mathbf{\text{WH}}}\mathbb{E}\left(\log \right(\mathbb{P}(A^o,A^u|\mathbf{\text{W}}\mathbf{\text{H}}\left)\right)|A^o,\mathbf{\text{W}}{\mathbf{\text{H}}}^{(t-\mathrm{1})}),\end{array}$$

where WH^(t−1) is the product of the previous estimates of W and H and ℙ is a probability measure. This problem is equivalent to running a PMF algorithm on the following adjusted matrix (Zhang et al., 2006):

$$\begin{array}{}\text{(9)}& {\mathbf{\text{A}}}_{\mathrm{1}}=\mathbf{\text{C}}\odot \mathbf{\text{A}}+(\mathbf{\text{1}}-\mathbf{\text{C}})\odot \mathbf{\text{W}}{\mathbf{\text{H}}}^{(t-\mathrm{1})},\end{array}$$

with C_ij=1 if A_ij is known and C_ij=0 if A_ij is unknown, and 1 is a matrix of ones.

Recent work has looked into expanding on this approach to continuous weights, as seen in most PMF problems of real-time atmospheric mass spectrometry data (Yahaya et al., 2019; Yahaya, 2021; Yahaya et al., 2021). To handle the continuous case, a variation of Eq. (8) is maximized (Yahaya et al., 2019):

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\mathrm{argmax}}_{\mathbf{\text{WH}}}\mathbb{E}\left(\log \right(\mathbb{P}(\mathbf{\text{C}}\odot \mathbf{\text{A}},(\mathbf{\text{1}}-\mathbf{\text{C}})\\ & \odot {\mathbf{\text{A}}}_{\mathrm{theo}}\left|\mathbf{\text{W}}\mathbf{\text{H}}\right)\left)\right|\mathbf{\text{C}}\odot \mathbf{\text{A}},\mathbf{\text{W}}{\mathbf{\text{H}}}^{(t-\mathrm{1})}),\end{array}\end{array}$$

where C is now a weight matrix containing estimates of confidence as a value between 0 and 1 in a given data point, and A_theo is the theoretical true data. Maximizing Eq. (10) is equivalent to running any PMF algorithm on the matrix A₁ formed in Eq. (9).

To form the confidence matrix C from the uncertainty matrix Σ, Yahaya (2021) suggests scaling the weights (in this case $\frac{\mathrm{1}}{{\mathit{\sigma }}_{ij}}$) so that the maximum value is 1. However, previous testing has primarily focused on problems with binary weights (Yahaya et al., 2019; Yahaya, 2021; Yahaya et al., 2021).

1.6 Determining the number of factors

The squared weighted residual error, $Q={\sum }_{i=\mathrm{1}}^m{\sum }_{j=\mathrm{1}}^n{\left(\frac{{\mathbf{\text{A}}}_{ij}-{\sum }_{l=\mathrm{1}}^k{\mathbf{\text{W}}}_{il}{\mathbf{\text{H}}}_{lj}}{{\mathit{\sigma }}_{ij}}\right)}^{\mathrm{2}}$, can be used to determine whether a given solution either overfits or underfits the data (Ulbrich et al., 2009). If measurement errors are normally distributed, then Q will follow a χ² distribution with $mn-k(m+n)$ degrees of freedom. Thus, to avoid overfitting, the number of factors in a solution is chosen such that $Q\approx mn-k(m+n)$ (Paatero et al., 2002). When additional error is present, measuring the convergence of Q or the weighted residual error with additional factors can determine if these additional factors add much information to the model.

Another method used in determining the number of factors is the lack of rotational ambiguity of a solution (Ulbrich et al., 2009). Consider a simple case where the data matrix A is the product of two rank-two matrices $\mathbf{\text{W}}=[\mathit{a},\mathit{b}]$ and $\mathbf{\text{H}}=[\mathit{y},\mathit{z}{]}^T$, with $\mathit{a},\mathit{b},\mathit{y},$ and z being column vectors. An exact solution can also be obtained by finding $\mathbf{\text{W}}=[\mathit{c},\mathit{d},\mathit{b}]$ and $\mathbf{\text{H}}=[\mathit{y},\mathit{y},\mathit{z}]$, where $\mathit{c}+\mathit{d}=\mathit{a}$, in a process known as factor splitting (Ulbrich et al., 2009). Thus, a three-factor solution introduces rotational ambiguity, as any two factors c and d can be chosen as long as they add up to a. The same analysis can be seen by analyzing solutions with four or more factors, and we can conclude that a large amount of rotational ambiguity is a potential sign of overfitting.

1.7 Data

In this paper, we use the data from “Ambient measurements of highly oxidized gas-phase molecules during the Southern Oxidant and Aerosol Study (SOAS) 2013”, measuring highly oxidized multifunctional molecules (HOMs) over a forest site in Alabama from 22 June to 7 July 2013 (Massoli et al., 2018). The dataset contains mass spectra concentrations of 1059 different ions over 27 336 different time stamps. Additionally, initial uncertainties associated with each measurement are also included. PMF was applied to the data using the PMF2 algorithm described in Paatero (1997), checking solutions from 2 to 10 factors, where a six-factor solution was obtained (Massoli et al., 2018). We note that some of the uncertainties were artificially increased for this PMF2 analysis. The authors concluded that a significant portion of the secondary organic aerosol (SOA) was the result of interactions between biogenic and anthropogenic emissions (Massoli et al., 2018).

We use this six-factor solution as a reference solution and test whether the RHALS algorithm can recreate formulated factors and those found from PMF2. Analyses of results for different numbers of factors (other than the original six identified in Massoli et al., 2018) were not considered in order to maintain interpretability of the algorithm output. For reference, the PMF2 factor mass spectra and the time trends over all of the data are shown in Fig. 1. The factor time series, as well as the time series of the total mass concentration, is also shown in Fig. 2. Both plots show total concentration amounts over the entire time series and mass spectra, respectively. Discussion of the chemical interpretations of the data are not presented, and the scope of this effort is limited to the mathematical results from implementing the RHALS algorithm.

Figure 1 (a–f) Mass spectra profiles of the six PMF2 factors labeled by order; (g) total mass spectra concentration of the data.

Figure 2 Time series profiles of the six PMF2 factors, overlaid with the total time series concentration of the data. A rolling average is used, with values representing average concentrations over the previous 2 h. The individual time series are stacked on top of each other in order to compare the total PMF2 time series to the overall time series of the data.

2 Methods

In this section, we present the derivation of the basic weighted HALS algorithm in Sect. 2.1, a simple rotation algorithm in Sect. 2.2, our new external weighting algorithm in Sect. 2.3, and inclusion of L1 and L2 regularization in Sect. 2.4.

2.1 HALS algorithm

The derivation of the unweighted algorithm for HALS is detailed in Erichson et al. (2018). Here we present a similar derivation, taking into account uncertainties in the data that act as weights. Another derivation of weighted HALS is given in Ho (2008).

The HALS algorithm applies block coordinate descent methods in order to minimize the cost function Q_j by minimizing a “block”, or an outer product of individual factors, of W and H at a certain time while keeping the other factors fixed (Erichson et al., 2018).

$$\begin{array}{}\text{(11)}& Q_j=\left|\right|({\mathbf{\text{R}}}_j-{\mathbf{\text{W}}}_{(:,j)}{\mathbf{\text{H}}}_{(j,:)})\oslash \mathbf{\Sigma }|{|}_F^{\mathrm{2}}\end{array}$$

In Eq. (11), R_j is the jth residual, with

$$\begin{array}{}\text{(12)}& {\mathbf{\text{R}}}_j=\mathbf{\text{A}}-\sum _{i\ne j}^k{\mathbf{\text{W}}}_{(:,i)}{\mathbf{\text{H}}}_{(i,:)}=\mathbf{\text{A}}-\mathbf{\text{W}}\mathbf{\text{H}}+{\mathbf{\text{W}}}_{(:,j)}{\mathbf{\text{H}}}_{(j,:)}.\end{array}$$

${\mathbf{\text{W}}}_{(:,j)}$ is the jth column (or jth factor) of W, ${\mathbf{\text{H}}}_{(j,:)}$ is the jth row and factor of H, and k is the number of factors that the algorithm is aiming to find. As defined previously, Σ contains the uncertainties associated with each element in the data matrix A.

To derive update rules for HALS, partial derivatives of Eq. (11) are taken with respect to the factors ${\mathbf{\text{W}}}_{(:,j)}$ and ${\mathbf{\text{H}}}_{(j,:)}$. In order to incorporate the uncertainties (Σ), we present a variation on the derivation presented in Erichson et al. (2018) by considering just a row (i) and column (p) of the weighted residual.

$$\begin{array}{}\text{(13)}& {\displaystyle }{Q}_j^i=\left|\right|({\mathbf{\text{R}}}_{j(i,:)}-{\mathbf{\text{W}}}_{ij}{\mathbf{\text{H}}}_{(j,:)}){\mathbf{\Sigma }}_i^{-\mathrm{1}}|{|}_F^{\mathrm{2}}\text{(14)}& {\displaystyle }{Q}_j^p=\left|\right|{\mathbf{\Sigma }}_p^{-\mathrm{1}}({\mathbf{\text{R}}}_{j(:,p)}-{\mathbf{\text{H}}}_{jp}{\mathbf{\text{W}}}_{(:,j)})|{|}_F^{\mathrm{2}}\end{array}$$

In Eqs. (13) and (14), Σ_i are Σ_p are diagonal matrices (of size m and n, respectively), with the diagonal elements corresponding to the elements of the ith row (for Σ_i) and pth column (for Σ_p) of Σ. Expanding Eq. (13) and Eq. (14) using the fact that $\left|\right|\mathbf{\text{X}}|{|}_F^{\mathrm{2}}=\mathrm{Tr}({\mathbf{\text{X}}}^T\mathbf{\text{X}})$, where Tr(A) denotes the trace of the matrix A, we get

$$\begin{array}{}\text{(15)}& {\displaystyle }\begin{array}{rl}{Q}_j^i& =\mathrm{Tr}({\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{R}}}_{j(i,:)}^T{\mathbf{\text{R}}}_{j(i,:)}{\mathbf{\Sigma }}_i^{-\mathrm{1}}\\ & -\mathrm{2}{\mathbf{\text{W}}}_{ij}{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{R}}}_{j(i,:)}^T{\mathbf{\text{H}}}_{(j,:)}{\mathbf{\Sigma }}_i^{-\mathrm{1}}\\ & +{\mathbf{\text{W}}}_{ij}^{\mathrm{2}}{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{H}}}_{(j,:)}^T{\mathbf{\text{H}}}_{(j,:)}{\mathbf{\Sigma }}_i^{-\mathrm{1}}),\end{array}\text{(16)}& {\displaystyle }\begin{array}{rl}{Q}_j^p& =\mathrm{Tr}({\mathbf{\text{R}}}_{j(:,p)}^T{\mathbf{\Sigma }}_p^{-\mathrm{1}}{\mathbf{\Sigma }}_p^{-\mathrm{1}}{\mathbf{\text{R}}}_{j(:,p)}\\ & -\mathrm{2}{\mathbf{\text{H}}}_{jp}{\mathbf{\text{R}}}_{j(:,p)}^T{\mathbf{\Sigma }}_p^{-\mathrm{1}}{\mathbf{\Sigma }}_p^{-\mathrm{1}}{\mathbf{\text{W}}}_{(:,j)}\\ & +{\mathbf{\text{H}}}_{jp}^{\mathrm{2}}{\mathbf{\text{W}}}_{(:,j)}^T{\mathbf{\Sigma }}_p^{-\mathrm{1}}{\mathbf{\Sigma }}_p^{-\mathrm{1}}{\mathbf{\text{W}}}_{(:,j)}).\end{array}\end{array}$$

Differentiating with respect to W_ij and H_jp,

$$\begin{array}{}\text{(17)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial Q_j^i}{\partial {\mathbf{\text{W}}}_{ij}}}& =-\mathrm{2}\mathrm{Tr}\left({\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{R}}}_{j(i,:)}^T{\mathbf{\text{H}}}_{(j,:)}{\mathbf{\Sigma }}_i^{-\mathrm{1}}\right)\\ & +\mathrm{2}{\mathbf{\text{W}}}_{ij}\mathrm{Tr}\left({\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{H}}}_{(j,:)}^T{\mathbf{\text{H}}}_{(j,:)}{\mathbf{\Sigma }}_i^{-\mathrm{1}}\right),\end{array}\text{(18)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial Q_j^p}{\partial {\mathbf{\text{H}}}_{jp}}}& =-\mathrm{2}\mathrm{Tr}\left({\mathbf{\text{R}}}_{j(:,p)}^T{\mathbf{\Sigma }}_p^{-\mathrm{1}}{\mathbf{\Sigma }}_p^{-\mathrm{1}}{\mathbf{\text{W}}}_{(:,j)}\right)\\ & +\mathrm{2}{\mathbf{\text{H}}}_{jp}\mathrm{Tr}\left({\mathbf{\text{W}}}_{(:,j)}^T{\mathbf{\Sigma }}_p^{-\mathrm{1}}{\mathbf{\Sigma }}_p^{-\mathrm{1}}{\mathbf{\text{W}}}_{(:,j)}\right).\end{array}\end{array}$$

To eliminate the matrix traces, it is easy to show that Eqs. (17) and (18) can be rewritten as

$$\begin{array}{}\text{(19)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial Q_j^i}{\partial {\mathbf{\text{W}}}_{ij}}}& =-\mathrm{2}{\mathbf{\text{R}}}_{j(i,:)}({\mathbf{\text{H}}}_{(j,:)}^T\oslash ({\mathbf{\Sigma }}_{(i,:)}^T\odot {\mathbf{\Sigma }}_{(i,:)}^T\left)\right)\\ & +\mathrm{2}{\mathbf{\text{W}}}_{ij}{\mathbf{\text{H}}}_{(j,:)}({\mathbf{\text{H}}}_{(j,:)}^T\oslash ({\mathbf{\Sigma }}_{(i,:)}^T\odot {\mathbf{\Sigma }}_{(i,:)}^T\left)\right),\end{array}\text{(20)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial Q_j^p}{\partial {\mathbf{\text{H}}}_{jp}}}& =-\mathrm{2}{\mathbf{\text{R}}}_{j(:,p)}^T({\mathbf{\text{W}}}_{(:,j)}\oslash ({\mathbf{\Sigma }}_{(:,p)}\odot {\mathbf{\Sigma }}_{(:,p)}\left)\right)\\ & +\mathrm{2}{\mathbf{\text{H}}}_{jp}{\mathbf{\text{W}}}_{(:,j)}^T({\mathbf{\text{W}}}_{(:,j)}\oslash ({\mathbf{\Sigma }}_{(:,p)}\odot {\mathbf{\Sigma }}_{(:,p)}\left)\right).\end{array}\end{array}$$

Setting Eqs. (19) and (20) to zero and solving for the factor values yield

$$\begin{array}{}\text{(21)}& {\displaystyle }{\mathbf{W}}_{ij}={\displaystyle \frac{{\mathbf{\text{R}}}_{j(i,:)}({\mathbf{\text{H}}}_{(j,:)}^T\oslash ({\mathbf{\Sigma }}_{(i,:)}^T\odot {\mathbf{\Sigma }}_{(i,:)}^T\left)\right)}{{\mathbf{\text{H}}}_{(j,:)}({\mathbf{\text{H}}}_{(j,:)}^T\oslash ({\mathbf{\Sigma }}_{(i,:)}^T\odot {\mathbf{\Sigma }}_{(i,:)}^T\left)\right),}}\text{(22)}& {\displaystyle }{\mathbf{H}}_{jp}={\displaystyle \frac{{\mathbf{\text{R}}}_{j(:,p)}^T({\mathbf{\text{W}}}_{(:,j)}\oslash ({\mathbf{\Sigma }}_{(:,p)}\odot {\mathbf{\Sigma }}_{(:,p)}\left)\right)}{{\mathbf{\text{W}}}_{(:,j)}^T({\mathbf{\text{W}}}_{(:,j)}\oslash ({\mathbf{\Sigma }}_{(:,p)}\odot {\mathbf{\Sigma }}_{(:,p)}\left)\right).}}\end{array}$$

Substituting Eq. (12) into Eqs. (19) and (20) yields the following update rules:

$$\begin{array}{}\text{(23)}& {\displaystyle }{\mathbf{\text{W}}}_{ij}\leftarrow {\left[{\mathbf{\text{W}}}_{ij}+{\displaystyle \frac{\begin{array}{c}{\mathbf{\text{A}}}_{(i,:)}({\mathbf{\text{H}}}_{(j,:)}^T\oslash ({\mathbf{\Sigma }}_{(i,:)}^T\odot {\mathbf{\Sigma }}_{(i,:)}^T\left)\right)\\ -{\mathbf{\text{W}}}_{(i,:)}\mathbf{\text{H}}({\mathbf{\text{H}}}_{(j,:)}^T\oslash ({\mathbf{\Sigma }}_{(i,:)}^T\odot {\mathbf{\Sigma }}_{(i,:)}^T\left)\right)\end{array}}{{\mathbf{\text{H}}}_{(j,:)}({\mathbf{\text{H}}}_{(j,:)}^T\oslash ({\mathbf{\Sigma }}_{(i,:)}^T\odot {\mathbf{\Sigma }}_{(i,:)}^T\left)\right)}}\right]}_{+},\text{(24)}& {\displaystyle }{\mathbf{\text{H}}}_{jp}\leftarrow {\left[{\mathbf{\text{H}}}_{jp}+{\displaystyle \frac{\begin{array}{c}{\mathbf{\text{A}}}_{(:,p)}^T({\mathbf{\text{W}}}_{(:,j)}\oslash ({\mathbf{\Sigma }}_{(:,p)}\odot {\mathbf{\Sigma }}_{(:,p)}\left)\right)\\ -{\mathbf{\text{H}}}_{(:,p)}^T{\mathbf{\text{W}}}^T({\mathbf{\text{W}}}_{(:,j)}\oslash ({\mathbf{\Sigma }}_{(:,p)}\odot {\mathbf{\Sigma }}_{(:,p)}\left)\right)\end{array}}{{\mathbf{\text{W}}}_{(:,j)}^T({\mathbf{\text{W}}}_{(:,j)}\oslash ({\mathbf{\Sigma }}_{(:,p)}\odot {\mathbf{\Sigma }}_{(:,p)}\left)\right)}}\right]}_{+},\end{array}$$

where the operator $[\cdot {]}_{+}$ projects all negative update values to 0. In practice, the authors in Erichson et al. (2018) utilize the following simplified form of Eqs. (23) and (24):

$$\begin{array}{}\text{(25)}& {\displaystyle }{\mathbf{W}}_{ij}\leftarrow {\left[{\mathbf{\text{W}}}_{ij}-{\displaystyle \frac{{\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}}{Q}_j^i}{{\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}}^{\mathrm{2}}{Q}_j^i}}\right]}_{+},\text{(26)}& {\displaystyle }{\mathbf{H}}_{jp}\leftarrow {\left[{\mathbf{\text{H}}}_{jp}-{\displaystyle \frac{{\mathrm{\nabla }}_{{\mathbf{\text{H}}}_{jp}}{Q}_j^p}{{\mathrm{\nabla }}_{{\mathbf{\text{H}}}_{jp}}^{\mathrm{2}}{Q}_j^p}}\right]}_{+}.\end{array}$$

Thus, one can add additional auxiliary functions to the cost function Q, such as regularization and rotation terms, and add them to the update rules based on the new gradient and Hessian values. Furthermore, writing the update rules as in Eqs. (25) and (26) includes the calculation of the projected gradient. This can be used as a stopping condition criterion that avoids the computational costs of calculating other convergence statistics (Erichson et al., 2018).

We note that when the uncertainties are equal to 1, Σ can be disregarded in Eqs. (23) and (24), and the update rules are identical to those in Erichson et al. (2018). In Eqs. (23) and (24), the Hessians $\mathbf{\text{H}}({\mathbf{\text{H}}}_{(j,:)}^T\oslash ({\mathbf{\Sigma }}_{(i,:)}^T\odot {\mathbf{\Sigma }}_{(i,:)}^T\left)\right)$ and ${\mathbf{\text{W}}}^T({\mathbf{\text{W}}}_{(:,j)}\oslash ({\mathbf{\Sigma }}_{(:,p)}\odot {\mathbf{\Sigma }}_{(:,p)}\left)\right)$ should be found prior to multiplication by ${\mathbf{\text{W}}}_{(i,:)}$ and ${\mathbf{\text{H}}}_{(:,p)}^T$, respectively, to minimize computational costs. We do not preallocate the products WH and H^TW^T, which is also not done in Erichson et al. (2018).

2.2 Rotational considerations

Rotating solutions to induce slight variations in the output of a PMF algorithm may be necessary, especially in cases where the interpretation of the solutions yields some unrealistic results (e.g., a factor is zero during a period in which it is expected to be present, or two factors appear mixed in their time series and/or spectra). As mentioned in Sect. 1.4, we do not attempt to constrain these “rotations” to be norm preserving. However, it is possible to find approximate rotations. As detailed in Paatero and Hopke (2009), for a specific factor value W_ij or H_jp, an auxiliary term can be added to the cost function to pull the component towards a set value ${\mathbf{\text{W}}}_{ij}^{*}$ or ${\mathbf{\text{H}}}_{ij}^{*}$. Defining

$$\begin{array}{}\text{(27)}& Q_{{\mathrm{aux}}_{\mathit{\text{W}}}}^{ij}=s({\mathbf{\text{W}}}_{ij}-{\mathbf{\text{W}}}_{ij}^{*}{)}^{\mathrm{2}}\text{ or }{Q}_{{\mathrm{aux}}_{\mathit{\text{H}}}}^{pj}=s({\mathbf{\text{H}}}_{jp}-{\mathbf{\text{H}}}_{jp}^{*}{)}^{\mathrm{2}},\end{array}$$

with s determining the strength of the pull, we find ${\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}}{Q}_{{\mathrm{aux}}_W}^{ij}=\mathrm{2}s({\mathbf{\text{W}}}_{ij}-{\mathbf{\text{W}}}_{ij}^{*})$, ${\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}}^{\mathrm{2}}{Q}_{{\mathrm{aux}}_W}^{ij}=\mathrm{2}s$, ${\mathrm{\nabla }}_{{\mathbf{\text{H}}}_{j,p}}{Q}_{{\mathrm{aux}}_H}^{pj}=\mathrm{2}s({\mathbf{\text{H}}}_{jp}-{\mathbf{\text{H}}}_{jp}^{*})$, and ${\mathrm{\nabla }}_{{\mathbf{\text{H}}}_{jp}}^{\mathrm{2}}{Q}_{{\mathrm{aux}}_H}^{jp}=\mathrm{2}s$. Adding these terms to the gradient and Hessian in Eqs. (25) and (26), we can derive new update rules for W_ij and H_jp.

To probe possible rotations of entire factors, we introduce “pulling equations” to the cost function, which pull the elements of W and H to the desired rotations:

$$\begin{array}{}\text{(28)}& {\displaystyle }{Q}_{{\mathrm{aux}}_{\mathit{\text{H}}}}^{\left(j\right)}=a{\left(\sum _{i=\mathrm{1}}^n(\mathrm{1}-{\mathbf{\text{H}}}_{ji})\right)}^{\mathrm{2}},\text{(29)}& {\displaystyle }{Q}_{{\mathrm{aux}}_{\mathit{\text{W}}}}^{\left(j\right)}=b{\left(\sum _{i=\mathrm{1}}^m{\mathbf{\text{W}}}_{ij}\right)}^{\mathrm{2}}.\end{array}$$

These equations correspond to r<0, with both a and b positive. For r>0, W and H would be interchanged. These pulling equations are similar to those introduced in Paatero et al. (2002), although their meanings are slightly changed.

Taking derivatives of Eqs. (28) and (29),

$$\begin{array}{}\text{(30)}& {\displaystyle }{\mathrm{\nabla }}_{{\mathbf{\text{H}}}_{jp}}{Q}_{{\mathrm{aux}}_H}^{jp}=-a(\mathrm{1}-{\mathbf{\text{H}}}_{jp}),\text{(31)}& {\displaystyle }{\mathrm{\nabla }}_{{\mathbf{\text{H}}}_{j}p}^{\mathrm{2}}{Q}_{{\mathrm{aux}}_H}^{jp}=a,\text{(32)}& {\displaystyle }{\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}}{Q}_{{\mathrm{aux}}_W}^{ji}=b{\mathbf{\text{W}}}_{ij},\text{(33)}& {\displaystyle }{\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}}^{\mathrm{2}}{Q}_{{\mathrm{aux}}_W}^{ji}=b,\end{array}$$

we add these terms to the gradient and Hessian in Eqs. (25) and (26) to receive new update rules. Optimally, we would expect these equations to pull the values of H towards more positive values and the values of W towards zero, which mirrors the effect of setting r to be less than zero in Eq. (6). The end result may be different solutions that are more realistic and interpretable.

2.3 External weighting

To perform external weighting, we first find $\widehat{\mathbf{\text{A}}}=\mathbf{\text{A}}\oslash \mathbf{\Sigma }$ (i.e., we divide the data elementwise by the uncertainties). Note that the uncertainty matrix Σ must only contain nonzero real entries, and therefore the algorithm cannot handle problems with binary weights (such as PMF problems with missing entries). One approach to PMF with binary weights is the expectation maximization (EM) approach detailed in the section titled “Expectation maximization” (Yahaya et al., 2021; Zhang et al., 2006).

After elementwise division, RHALS is applied to $\widehat{\mathbf{\text{A}}}$ to form the estimated scaled factor matrices $\widehat{\mathbf{\text{W}}}$ and $\widehat{\mathbf{\text{H}}}$. The unscaled estimates W and H are found by the relation

$$\begin{array}{}\text{(34)}& \mathbf{\text{W}}\mathbf{\text{H}}=\left(\widehat{\mathbf{\text{W}}}\widehat{\mathbf{\text{H}}}\right)\odot \mathbf{\Sigma }.\end{array}$$

W and H can then be found iteratively via alternating least squares:

$$\begin{array}{}\text{(35)}& {\displaystyle }\mathbf{\text{W}}=\left(\right(\widehat{\mathbf{\text{W}}}\widehat{\mathbf{\text{H}}})\odot \mathbf{\Sigma }){\mathbf{\text{H}}}^T(\mathbf{\text{H}}{\mathbf{\text{H}}}^T{)}^{-\mathrm{1}}\text{ or }\mathbf{\text{W}}=(\left(\widehat{\mathbf{\text{W}}}\widehat{\mathbf{\text{H}}}\right)\odot \mathbf{\Sigma }){\mathbf{\text{H}}}^{†},\text{(36)}& {\displaystyle }\mathbf{\text{H}}=\left({\mathbf{\text{W}}}^T\mathbf{\text{W}}{)}^{-\mathrm{1}}{\mathbf{\text{W}}}^T\right(\left(\widehat{\mathbf{\text{W}}}\widehat{\mathbf{\text{H}}}\right)\odot \mathbf{\Sigma })\text{ or }\mathbf{\text{H}}={\mathbf{\text{W}}}^{†}(\left(\widehat{\mathbf{\text{W}}}\widehat{\mathbf{\text{H}}}\right)\odot \mathbf{\Sigma }).\end{array}$$

Here, W^† and H^† denote the pseudo-inverses of W and H, and all negative elements are set to zero after the factor matrix is updated. Mathematically, the two methods for calculating W and H detailed in Eqs. (35) and (36), respectively, are identical, as long as the ranks of the factor matrices are equal to k. Evidence of this is briefly detailed in Appendix B. Of course, if the factor matrices were of a lower rank (or contained very small singular values), then a lower-rank factorization should first be found instead. Since both algorithms are identical mathematically, one could be favorable if it provided a speed advantage, depending on the computational efficiency of the pseudo-inverse algorithm called. Both update rules are 𝒪(mnk) (Feng et al., 2018), and running this code in MATLAB on a single CPU, we found the update rules using the pseudo-inverses were faster.

To begin the iteration, we have to initialize either W or H. The most intuitive way to do this is by assuming the unweighted factors are similar to the weighted factors and setting H (or W) to the weighted values but scaled to the magnitude of the original data. Since each entry of the data matrix is approximated as ${\text{A}}_{ij}={\sum }_{l=\mathrm{1}}^k{\mathbf{\text{W}}}_{il}{\mathbf{\text{H}}}_{lj}$, and assuming relatively equal magnitudes for W and H, the appropriate scale factor should be $\sqrt{\frac{\stackrel{\mathrm{‾}}{\mathbf{\text{A}}}}{k}}$, where $\stackrel{\mathrm{‾}}{\mathbf{\text{A}}}$ denotes the element mean of A. When W and H contain values with greatly differing magnitudes, perhaps a different scaling factor should be used, although this has not been extensively tested. Thus, the algorithm should be initialized by

$$\begin{array}{}\text{(37)}& {\mathbf{\text{H}}}_{\mathrm{0}}=\left(\sqrt{{\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathbf{\text{A}}}}{k{\stackrel{\mathrm{‾}}{\mathbf{\text{H}}}}^{\mathrm{2}}}}}\right)\widehat{\mathbf{\text{H}}}.\end{array}$$

In Eq. (37), $\stackrel{\mathrm{‾}}{\mathbf{\text{H}}}$ denotes the element mean of $\widehat{\mathbf{\text{H}}}$. Measuring the change in W and H between iterations, the algorithm typically converges relatively quickly – within 20 to 40 iterations. In practice, one can use L2 regularization in the external weighting steps, equal to 0.01 for our data, as the least squares method may become increasingly ill-posed as the number of factors increases. This value may need to be altered based on the magnitude of the values in the data and the number of factors. However, we found that adding L2 regularization lowered the similarity of the factors to the given factors from the solution using PMF2.

We used this method for all externally weighted algorithms tested in Sect. 4. Theoretically, any PMF algorithm could be used for the post-processing step. This may become relevant, since as noted in Sect. 1.2.3, the nonnegative ALS method described above may have convergence issues for certain factorization problems (Gillis, 2020).

2.4 Regularization

As in the RHALS algorithm presented in Erichson et al. (2018), we add L1 and L2 regularization. L1 regularization is added to the cost function through the L1 norm, $\left|\right|\cdot |{|}_{\mathrm{1}}$, which is the maximum sum of the absolute value of the components in a given column in the factor matrices. L2 regularization is added through the L2 (or Euclidean) norm, $\left|\right|\cdot |{|}_{\mathrm{2}}^{\mathrm{2}}$. L1 regularization is typically added to reduce sparsity, while L2 regularization is added to control the Euclidean norms of the factors and avoid overfitting ill-posed problems (Erichson et al., 2018). The cost function now becomes

$$\begin{array}{}\text{(38)}& \begin{array}{rl}& \left|\right|({\mathbf{\text{R}}}_j-{\mathbf{\text{W}}}_{(:,j)}{\mathbf{\text{H}}}_{(j,:)})\oslash \mathbf{\Sigma }|{|}_F^{\mathrm{2}}+\mathit{\alpha }|\left|{\mathbf{\text{W}}}_{(:,j)}\right|{|}_{\mathrm{1}}\\ & +\mathit{\beta }\left|\right|{\mathbf{\text{H}}}_{(j,:)}|{|}_{\mathrm{1}}+\mathit{\gamma }|\left|{\mathbf{\text{W}}}_{(:,j)}\right|{|}_{\mathrm{2}}^{\mathrm{2}}+\mathit{\delta }\left|\right|{\mathbf{\text{H}}}_{(j,:)}|{|}_{\mathrm{2}}^{\mathrm{2}}.\end{array}\end{array}$$

In Eq. (38), $\mathit{\alpha },\mathit{\beta },\mathit{\gamma }$, and δ control the amount of regularization that is added. Typically α and β, as well as γ and δ, are set equal. ${\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}}{Q}_j^i=\mathit{\beta }+\mathrm{2}\mathit{\gamma }{\mathbf{\text{W}}}_{ij}$, ${\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}}^{\mathrm{2}}{Q}_j^i=\mathrm{2}\mathit{\gamma }$, and ${\mathrm{\nabla }}_{{\mathbf{\text{H}}}_{jp}}{Q}_j^p=\mathit{\alpha }+\mathrm{2}\mathit{\delta }{\mathbf{\text{H}}}_{jp}$, ${\mathrm{\nabla }}_{{\mathbf{\text{H}}}_{jp}}^{\mathrm{2}}{Q}_j^p=\mathrm{2}\mathit{\delta }$. These terms are added to the update rules in Eqs. (25) and (26). Optimal parameter values for $\mathit{\alpha },\mathit{\beta },\mathit{\gamma }$, and δ in Eq. (38) are often found using an L-curve analysis, which measures the tradeoff between minimizing a norm and minimizing the residual error (Hansen and O'Leary, 1993).

In the SOAS dataset, the average data value is $\mathrm{1.9969}\times {\mathrm{10}}^{-\mathrm{5}}$ for the mass-to-charge ratio. When external weighting is applied, the average value of A⊘Σ is 1.4851. The L curves below are thus plotted using two different scales: one with a scale of 10⁻⁵ and one with a scale of 1. The deterministic, unweighted algorithm is first applied to a small dataset consisting of the first 100 rows and columns of the SOAS dataset, scaled to contain different average values, to produce the L curves shown in Fig. 3. Values of the regularization parameters are listed in the graphs, with “e” representing the exponential function with base 10, and the number to the right is the exponent (e.g., e03=10³). The optimal regularization value is one that is located in the bend of the graph, meaning it minimizes a solution’s norm while having low residual error. In our results, we put emphasis on minimizing the latter, so we choose regularization values towards the left side of the bend. Analyzing the figures, the L1 regularization parameters α and β are chosen between $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ for the scale of 10⁻⁵ and between 0.1 and 1 for the scale of 1. The L2 regularization parameters γ and δ are chosen around 0.001 for the scale of 10⁻⁵ and 100 for the scale of 1. We note that choosing much smaller regularization values does not drastically increase the norms of the solution, suggesting that regularization is not especially necessary for this problem. However, for more ill-posed problems, the L1 and L2 norms may become extremely large as the amount of regularization tends to zero (Hansen and O'Leary, 1993).

Figure 3 L curves for L1 and L2 regularization, applied to data with different average values. (a) L1 regularization for data with a magnitude of around 10⁻⁵, (b) L1 regularization for data with a magnitude of around 1, (c) L2 regularization for data with a magnitude of around 10⁻⁵, and (d) L2 regularization for data with a magnitude of around 1.

As a rule of thumb, the L1 regularization parameters should be chosen to be of the same order of magnitude as the data, and the L2 regularization parameters should be chosen to be around 100 times that. Note that neither the data matrix size nor the number of factors affect these optimal regularization parameter values. External and internal weighting will also not affect these values, as the magnitudes of the gradient and the Hessian will be around the same value.

An increasingly popular alternative to traditional regularization is nuclear norm regularization, which can be applied to matrix factorization without the need for rank constraints (Hu et al., 2013; Sun and Mazumder, 2013; Fornasier et al., 2011). The nuclear norm is defined as

$$\begin{array}{}\text{(39)}& \left|\right|\mathbf{\text{A}}|{|}_{*}=\sum _{i=\mathrm{1}}^n{\mathit{\sigma }}_i(\mathbf{\text{A}}),\end{array}$$

where σ_i(A) is the ith singular value of A. It is the convex envelope to the rank function, and thus finding a PMF solution that minimizes the nuclear norm also has a minimal rank (Hu et al., 2013). Two possible ways to apply nuclear norm regularization to weighted PMF are by the alternating direction method of multipliers (ADMM), as demonstrated in Sun and Mazumder (2013), and by reconstruction of the nuclear norm into a Frobenius norm, which can then be treated as L2 regularization (Fornasier et al., 2011). However, both these approaches involve key computations with the large, low-rank product WH and may be computationally expensive to implement. Furthermore, both algorithms still involve some arbitrary choice of rank by requiring a present amount of nuclear norm regularization. In traditional PMF, factor profiles are optimized from a predetermined number of factors. Our implementation and methodology are consistent with this approach.

3 Computational details

The computer used for these calculations has an 11th Gen Intel^® Core™ i7-1165G7 dual-quad-core CPU with a speed of 2.80 GHz. MATLAB™ R2021a is used for all calculations. While no algorithms detailed in this paper explicitly utilize multicore processing, they may rely on built-in MATLAB functions (such as svd()) that utilize multicore processing for increased efficiency.

4 Results

We measure the accuracy of RHALS using the weighted residual error of the final algorithm, again defined as $\sqrt{{\sum }_{i=\mathrm{1}}^m{\sum }_{j=\mathrm{1}}^n{\left(\frac{{\mathbf{\text{A}}}_{ij}-{\sum }_{l=\mathrm{1}}^k{\mathbf{\text{W}}}_{il}{\mathbf{\text{H}}}_{lj}}{{\mathit{\sigma }}_{ij}}\right)}^{\mathrm{2}}}$; the correlation coefficients between the time series factors produced by RHALS and those of the reference solution; and the cosine similarity between the mass spectra profiles produced by RHALS and the reference solution.

In the MU algorithm, the elements of W and H are initialized by taking the absolute value of random numbers drawn from a standard normal distribution. These values are then scaled by a factor of ${\stackrel{\mathrm{‾}}{\mathbf{\text{A}}}}^{\mathrm{2}}$. The ALS and HALS algorithms are initialized by the nonnegative double SVD (NNDSVD) approach detailed in Boutsidis and Gallopoulos (2008). Note that the latter approach will still vary between random seeds as the NNDSVD is also found from utilizing the randomization technique in Sect. 1.3.

All algorithms cease updates when a stopping criterion is met or the maximum number of iterations is reached. For the MU and ALS algorithms, the value of the cost function is calculated after each iteration, and the algorithm halts when the percent change in this value, divided by the initial cost function value, is less than a set tolerance. The HALS algorithm differs by using the percent change in the projected gradient in its stopping criterion (detailed in Lin, 2007). The projected gradient is defined as (Erichson et al., 2018)

$$\begin{array}{}\text{(40)}& {\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}}^P\left\{\begin{array}{l}{\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}}\text{ if }{\mathbf{\text{W}}}_{ij}>\mathrm{0}\\ min(\mathrm{0},{\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}})\text{ if }{\mathbf{\text{W}}}_{ij}=\mathrm{0}\end{array}\right.\end{array}$$

and similarly for ${\mathrm{\nabla }}_{{\mathbf{\text{H}}}_{jp}}^P$. When the percentage change in ${\sum }_{i=\mathrm{1}}^m{\sum }_{j=\mathrm{1}}^k({\mathrm{\nabla }}_{{\mathbf{\text{W}}}_{ij}}^P{)}^{\mathrm{2}}+{\sum }_{i=\mathrm{1}}^k{\sum }_{j=\mathrm{1}}^n({\mathrm{\nabla }}_{{\mathbf{\text{H}}}_{ij}}^P{)}^{\mathrm{2}}$ (divided by its initial value) is less than a specified tolerance, the algorithm halts. Below, the maximum number of iterations allowed is set to be 100 (which is rarely reached by any algorithm), and the tolerance is set to be $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ for the HALS and ALS algorithms. A tolerance of 10⁻⁵ is used for the MU algorithm to account for a larger initial cost from the random initialization.

When comparing factors from different algorithms, it is important to note that the ratio of values between a time series factor and a mass spectra factor might be slightly different. Consider two algorithms that compute the identical factorizations

$$\begin{array}{rl}& {\mathbf{\text{W}}}_{\mathrm{1}}{\mathbf{\text{H}}}_{\mathrm{1}}=\sum _{l=\mathrm{1}}^k{\mathbf{\text{W}}}_{(:,l)}^{\left(\mathrm{1}\right)}{\mathbf{\text{H}}}_{(l,:)}^{\left(\mathrm{1}\right)}\text{ and }\\ & {\mathbf{\text{W}}}_{\mathrm{2}}{\mathbf{\text{H}}}_{\mathrm{2}}=\sum _{l=\mathrm{1}}^k{\mathbf{\text{W}}}_{(:,l)}^{\left(\mathrm{2}\right)}{\mathbf{\text{H}}}_{(l,:)}^{\left(\mathrm{2}\right)}=\sum _{l=\mathrm{1}}^k\left(c_l{\mathbf{\text{W}}}_{(:,l)}^{\left(\mathrm{1}\right)}\right)\left({\displaystyle \frac{\mathrm{1}}{c_l}}{\mathbf{\text{H}}}_{(l,:)}^{\left(\mathrm{1}\right)}\right),\end{array}$$

with c_l being an arbitrary positive constant. It can be seen above that ${\mathbf{\text{W}}}_{(:,l)}^{\left(\mathrm{2}\right)}$ is identical to ${\mathbf{\text{W}}}_{(:,l)}^{\left(\mathrm{1}\right)}$, scaled in magnitude by c_l. It would appear that one factorization yields a stronger signal for factor l, despite the fact that the factorizations are identical. Thus, when comparing individual factors between algorithms, the time series concentrations are cumulative over the entire mass spectra, and the mass spectra concentrations are cumulative over the entire time series (as also presented with the PMF2 results). Specifically, we plot the time series of factor l as the sum over the columns of the outer product ${\mathbf{\text{W}}}_{(:,l}{\mathbf{\text{H}}}_{(l,:)}$, and we plot the mass spectra of factor l as the sum over the rows. To avoid this issue, one can add scaling coefficients to rescale the factors at each iteration so that they will approximately have the same magnitude (Lu and Wu, 2004). This has not yet been implemented in our code.

4.1 Computational efficiency

Table 1 shows the number of operations required at each step (update of W and H). Preprocessing steps, such as finding A⊘Σ and Σ⊙Σ, are excluded from the table. External weighting eliminates almost all elementwise operations needed, which is a slow memory-bound operation (Jia et al., 2020). However, systems with an abundance of free memory and/or GPUs may find that internal weighting methods are sufficiently quick for small- or medium-sized problems. Furthermore, the performance of matrix operations within these algorithms may vary between programming languages and libraries, such as between the NumPy, NumExpr, and Theano libraries in Python (Bergstra et al., 2010).

Table 1 Comparison of the number of operations per step for each algorithm used. We also track the number of elementwise multiplications and divisions (implemented in MATLAB as .* and ./) in the third column.

Analyzing Table 1 allows us to see why, in practice, it takes longer for internally weighted ALS and HALS to converge than internally weighted MU, as the former will have more matrix multiplication operations and elementwise operations for our target rank, k=6. However, when the algorithms are internally weighted, ALS and HALS become much more computationally efficient, allowing runtime to diminish as well.

To show numerically how the cost of RHALS scales with the size of the system, we consider the size of the data matrix vs. the runtime of the RHALS algorithm, as well as the deterministic internally weighted HALS and MU algorithms. Figure 4 shows how the computational costs of the different algorithms vary with the number of rows (and columns) in a square data matrix. The number of factors is set to five, and the data are formed by multiplying two rank-five matrices filled with sampled values from the PMF2 solution of the SOAS dataset. Uncertainty values are also randomly sampled from the given uncertainty data, and noise is added to the data matrix through random normal values centered around 0 with standard deviations equal to the uncertainties. L1 regularization is set to 1 for all factors, and L2 regularization is set to 50.

Figure 4 Size of data matrix vs. runtime (in seconds) in RHALS algorithm, performed with three different random seeds (containing the absolute value of random normal variables). The x axis shows the base 10 log of n, and the y axis plots the base 10 log of runtime. The median runtime is plotted, with the error bars plotting the maximum and minimum runtimes. Note that runtime initially decreases due to a coincidental reduction in the number of steps to convergence, but problems of environmental interest are generally located on the right of the graph or beyond its right edge.

As seen in Fig. 4, the logarithms of size and runtime are approximately linear for all three algorithms. RHALS clearly outperforms MU and HALS in terms of efficiency, requiring roughly 5 % and 1 % of the runtime of the MU and deterministic HALS algorithms, respectively. It should also be noted that the randomization step, which accounts for about 5 % of the algorithm computational costs, can be parallelized using GPUs (Erichson et al., 2018). The post-algorithm external weighting step can also be parallelized as it only involves matrix multiplication and an SVD (Lahabar and Narayanan, 2009). This post-weighting step accounts for around 60 % of the algorithm's computational costs.

4.2 Simple case

We now present the RHALS algorithm for a small test case. To do this, we form a data matrix by combining underlying “true” factors and random noise. We form a 5952 by 400 matrix using the first three factors from the PMF2 solution to the SOAS dataset. We choose the time stamps between 9385 and 15 336 from the PMF2 time series factor matrix W_PMF2 as the PMF2 factors are clearly distinguishable from each other in this time interval. We choose the first 400 mass spectra profiles from the mass spectra factor matrix H_PMF2 as the bulk of the mass spectra concentrations lies within this range.

A data matrix is then formed by combining the factors and adding random normal noise with a standard deviation equal to that of the uncertainty of the data:

$$\begin{array}{}\text{(41)}& {\mathbf{\text{A}}}_{ij}^{\mathrm{test}}=max(\mathrm{0},{\mathbf{\text{W}}}_{(i,:)}^{\mathrm{test}}{\mathbf{\text{H}}}_{(:,j)}^{\mathrm{test}}+N(\mathrm{0},{\mathbf{\Sigma }}_{ij}\left)\right).\end{array}$$

Next, we run RHALS over 20 trials, each with a different initialization, in order to obtain many different solutions that exist in the solution space. Here, initializations are formed using the randomized SVD method described earlier, which generated sufficiently different solutions. Figure 5 shows the weighted residual error of the RHALS algorithm over 20 trials.

Figure 5 RHALS error over 20 trials. Mean: 1476.60, standard deviation: 20.47.

The average weighted error over 20 trials is 1.4749×10³, and the algorithm took an average time of 0.0561 s to converge with an average of 48.15 block coordinate descent steps per trial. As one can see, solutions can vary over different trials, and a smaller weighted error could be used as justification of one solution over another.

Next, we compare the similarities of the converged factors and the original factors that formed the test case. The averages of the similarities between the three factors are shown in Fig. 6. A cosine similarity is used for the mass spectra, while the correlation coefficient is used for the time series.

Figure 6 (a) Average time series similarity and (b) average mass spectra similarity to formed factors in the small test case.

As one can see, the RHALS algorithm recreates the mass spectra and time series factors almost perfectly for all trials tested, regardless of the weighted residual error. In practice, a cosine similarity over 0.95 or a correlation coefficient over 0.90 between factor profiles represents almost identically interpretable solutions. Any trial is viable to be chosen as a “good” solution. Each solution yields a mass spectra similarity over 0.994 and a time series correlation over 0.974. There is low variance in the similarity metrics among trials, and a higher weighted error corresponds with less important parts of the factor.

4.3 Large dataset

4.3.1 Comparing different algorithms

Next, we analyze the complete dataset and compare the RHALS factors to all PMF2 factors. Table 2 shows a table of diagnostics for the HALS, RHALS, ALS, and MU algorithms applied to the complete SOAS dataset, averaged over three trials. As one can see, algorithms with external weighting demonstrate a dramatic reduction in computational costs, albeit at the cost of a larger error. The ALS algorithms were by far the slowest, and the internally weighted HALS algorithm yielded the lowest average weighted error. The internally weighted MU algorithm produced comparable results to the internally weighted HALS algorithm less than two-thirds of the time; however, both the internally and the externally weighted cases took the most steps to converge. The two fastest algorithms were the externally weighted HALS and RHALS algorithms, taking 1.01 and 0.50 s, respectively.

Table 2 Average statistics of different algorithms over five different SVD initializations, with a tolerance of 10⁻⁴ (for MU, tolerance is 10⁻⁵ and a random initialization), L1 regularization of 1, and L2 regularization of 50.

PMF factor profiles are often visually analyzed and compared to known candidate profiles for identification. Thus, the question arises as to whether or not an algorithm utilizing external weighting still contains interpretable factors in light of the decreased accuracy. Figure 7 shows the cumulative time series from one solution of the externally and internally weighted HALS algorithm, with the factors produced by the externally weighted algorithm and the internally weighted algorithm overlaid. Figure 8 shows the cumulative mass spectra of this solution, with the factors from the externally weighted algorithm and the internally weighted algorithm laid out side by side. Similar graphs for the ALS and MU algorithms are shown in Figs. A4 and A5 respectively in Appendix A.

Figure 7 Comparison of 2 h rolling average of time series for externally weighted and internally weighted HALS factors. Externally weighted error: 7.1321×10³. Internally weighted error: 6.4841×10³.

Figure 8 Comparison of mass spectra for internally weighted (on the left) and externally weighted (on the right) HALS factors.

Upon visual inspection, the time series produced by the externally weighted factors have peaks and troughs at almost identical times, but the magnitudes of these peaks and troughs can vary throughout the factors. Interestingly, there exists a consistent difference in the magnitude of the time series between internal and external weighting. Specifically, the external weighting algorithms seem to consistently overpredict the concentrations of the time series of the second factor compared to the internally weighted algorithms. The difference is small compared to the concentrations of the factors, but this may lead to an over-interpretation of the importance of the second factor with an externally weighted algorithm. It should be noted that there exists a (similarly) large variation in relative factor signals between trials of internally weighted algorithms, although caution should be taken with regard to the possibility that external weighting introduces extra error in this analysis. More research is needed to understand the overestimation or underestimation of factor magnitudes in external weighting.

Comparing the mass spectra yields a similar analysis. Most factors of externally weighted algorithms share the same spikes of key ions as in the internally weighted ions, although sometimes at different magnitudes. Occasionally, some externally weighted mass spectra factors will look quite dissimilar to the corresponding internally weighted factor or contain noticeable spikes and divots in key ions, as seen in the second factor for the MU and HALS algorithm. It should also be noted that the solutions of the externally weighted algorithms could be rotated away from each other. These differences may be extreme enough to encourage users to utilize multiple trials to search for multiple solutions.

4.3.2 Comparison between expectation maximization and external weighting

To test the EM approach to uncertainty weighting, as mentioned in the section titled “Expectation maximization”, the weights σ_ij in the uncertainty matrix Σ are scaled so that ${max}_{ij}\left(\frac{\mathrm{1}}{{\mathit{\sigma }}_{ij}}\right)=\mathrm{1}$. Since the bulk computational component of the algorithm constructs the matrix A₁ in Eq. (9), the authors in Yahaya et al. (2019) and Yahaya et al. (2021) recommend updating A₁ only after convergence or a maximum number of iterations of 20 or 50. They also note that applying the expectation step too early in the algorithm led to poorer performance due to the number of errors in the estimates of W and H. In order to apply the projected gradient stopping criterion discussed in Sect. 4, we chose to reconstruct A₁ at fixed iterations – after 10 and 20 PMF steps for different experiments. The first construction of A₁ is evaluated at the 1st, 5th, and 10th step. We used the NNDSVD initialization in Boutsidis and Gallopoulos (2008) and also varied the tolerance of the stopping condition between 10⁻⁴ and 10⁻⁶. We summarize these results by presenting the range of average values across the different variations.

We present a comparison of external weighting and the EM algorithm in Table 3, using the ranges of values from the different experiments listed above. Each value is an average over 20 trials. We compare the convergence times of the algorithms, the number of steps, the weighted errors, and the similarity to the PMF2 solution for both W (correlation) and H (cosine similarity).

Table 3 Average statistics of external weighting (EW) versus expectation maximization (EM) algorithms over 20 trials. Internally weighted (IW) HALS is provided as a reference. Ranges of the values of the algorithms run using different variations of expectation maximization (EM) steps are presented. The correlation of the columns of W and similarity of the rows of H to the PMF2 solution are also listed.

As demonstrated in Table 3, some variations of the EM algorithm were able to outperform externally weighted HALS and RHALS in total time, as well as in weighted error. For instance, running the expectation maximization algorithm with the first calculation of A₁ taking place at the 10th step and recalculating A₁ after 20 additional steps until the convergence criterion with a tolerance of 10⁻⁵ was reached yielded an average weighted error of 6.92×10³ in 0.6703 s. This computational speed compares to RHALS, while the accuracy bests both RHALS and externally weighted HALS. However, no expectation maximization algorithm was as successful at recreating the PMF2 time series factors, as seen in column five of Table 3. Externally weighted RHALS and HALS also provided mass spectra factors with a higher similarity to the PMF2 factors, with the exception of two runs of the expectation maximization algorithm, one recalculating A₁ after 20 steps, with the first calculation at the fifth step, and the other after 10 steps, with the first calculation at the first step. These yielded average similarities of 0.8856 and 0.8861, respectively, although they both took over 2.90 s to run, much slower than any other algorithm tested besides internally weighted HALS.

We also tested how well each algorithm produced factors that were within a rotation of the PMF2 factors. As detailed in Sect. 1.4, the factor profiles of $\widehat{\mathbf{\text{W}}}\widehat{\mathbf{\text{H}}}={\mathbf{\text{WT}}}^{-\mathrm{1}}\mathbf{\text{TH}}$ for a square matrix T may be closer to the desired solution than the original factors W and H. To see the extent that W can be rotated towards W_PMF2, we find the matrix T that minimizes the total squared differences between H_PMF2 and TH and then find the average correlation between the rows of W_PMF2 and WT⁻¹. A symmetrical approach can be made to find the cosine similarity for H. For internally weighted HALS, externally weighted HALS, and RHALS, we found average post-rotation time series correlations to be 0.9391, 0.9021, and 0.8902, respectively, and post-rotation mass spectra similarities to be 0.9602, 0.9433, and 0.9518, respectively. When this approach was tested on HALS and RHALS using the EM algorithm, average post-rotation time series correlations varied within 0.8528–0.8819 and 0.8314–0.8544, respectively, and post-rotation mass spectra similarities within 0.9230–0.9363 and 0.9147–0.9287, respectively.

Ultimately, we found that external weighting recreated the PMF2 factors more consistently than EM, both before and after rotation. This may be due to the fact that the scaling of the weights used for the EM step is not perfectly analogous to creating a set of weights that represent the confidence of each data point. Thus, the EM method using this scaling may not capture the key error-weighted patterns in the data as well as external weighting.

4.3.3 Complete analysis of RHALS algorithm

Figure 9 shows the weighted residual error for a six-factor solution (equivalent to the number of factors in PMF2) over 20 different nonnegative SVD initializations. The average weighted residual error is 7.2743×10³, with a convergence time of 0.5163 s over 38.45 steps per trial. If the solutions differ significantly, weighted residual error is a useful metric for choosing a solution from the 20 trials. The 11th and 13th trials appear promising, with the lowest weighted residual errors of 7.1474×10³ and 7.1026×10³, respectively. The average similarities between RHALS factors and PMF2 factors are presented in Fig. 10, with the same similarity metrics as listed above.

Figure 9 RHALS error over 20 trials. Mean error: 7.2743×10³, standard deviation of error: 120.8. Mean time: 0.5163, standard deviation of time: 0.0292. Mean steps: 38.45, standard deviation of steps: 6.88.

Figure 10 (a) Average time series and (b) average mass spectra similarities between RHALS and PMF2.

We see that different initializations can lead to different solutions, in terms of both similarity to the given PMF2 solution and weighted error, suggesting that convergence to a global minimum is not always achieved. This further emphasizes the importance of using multiple initializations in order to find an optimal solution.

Only the 2nd, 11th, 13th, and 17th trials have an average cosine similarity between mass spectra factors over 0.95, along with a correlation over 0.90 between time series factors. The 13th trial holds the highest similarity scores, with a time series correlation of 0.9468 and a cosine similarity of 0.9787. As seen in Fig. 9, these solutions have a small weighted residual error, further justifying picking a solution with low weighted error.

Figure 11 (a) Average time series correlation, (b) average mass spectra similarity, and (c) average of time series and mass spectra similarities between the solutions of the different trials for RHALS.

Figure 12 RHALS mass spectra series factors. (a) First factor similarity with PMF2=0.9631, (b) second factor similarity with PMF2=0.9582, (c) third factor similarity with PMF2=0.9912, (d) fourth factor similarity with PMF2=0.9970, (e) fifth factor similarity with PMF2=0.9795, and (f) sixth factor similarity with PMF2=0.9835.

Figure 13 The 2 h rolling average of RHALS time series factors overlaid with PMF2 factors.

As we saw with the small test case, every RHALS solution contained factors similar to the true factors. However, only 4 out of 20 trials produced solutions surpassing 0.95 in mass spectra similarity and 0.90 in time series correlation for the entire dataset. Running 100 trials, it was found that RHALS found a solution matching these criteria 27 % of the time. Thus the probability of not finding a good solution in 10 trials would be around (1–0.27)¹⁰≈4.3 % and ≈0.2 % in 20 trials. This rate could vary depending on the number of factors and the rotational ambiguity of the problem, as in the “bad” trials RHALS may simply find a rotated version of the “true” solution.

In the potential case that interest lies in factors that occur most frequently, an alternative approach to picking a solution given multiple trials would be to focus on solutions in which factors are found repeatedly. In Fig. 11, we compare the average similarities of a given trial to the other 19 trials. Again, cosine similarity is used for the mass spectra, and the correlation coefficient is used for the time series. The bottom plot shows the two graphs averaged to give a total metric of the similarity of a solution to the other solutions.

Analyzing these figures, it appears that most solutions have about the same similarity to each other, with a time series correlation of around 0.87 and a mass spectra similarity of around 0.93. A few solutions, such as those in the trials between 6 and 9, can be ruled out as outliers due to a low similarity to the other solutions. While the 12th trial yields the highest similarity between it and the other solutions (0.9148 averaged between the time series correlation and the mass spectra similarity), no solution drastically outperforms any of the other solutions. Note that the solutions from the 2nd, 11th, 13th, and 17th trials all perform well in this analysis, with average similarities of 0.9028, 0.8993, 0.9016, and 0.9036, respectively.

From visual inspection of Figs. 9 and 11, the solution from the 13th trial is barely outperformed in the above analysis, while holding a clear lead in accuracy. Thus, it seems to be a rational choice to decide that this solution is the “best”. The solution from the 13th trial is plotted in Figs. 12 and 13. Figure 12 shows the cumulative mass spectra profiles from the RHALS solution, while Fig. 13 shows the cumulative time series profiles from the RHALS solution compared with those from the PMF2 solution. Again, the profiles are plotted as a sum over the entire time series and mass spectra, respectively.

The similarities of the RHALS factors to the “true” PMF2 factors are very promising; however, the magnitudes of some factors differ. Specifically, RHALS underestimates the significance of the first factor, while overestimating the significance of the second and fifth factors. Again, this additional error should be considered when interpreting RHALS factors.

Another potential concern from the RHALS algorithm is that the factors have a general bias towards higher or lower magnitudes. Specifically, since RHALS enforces nonnegativity merely by setting negative elements to zero, in addition to regularization, it could be hypothesized that RHALS would produce factors of lower magnitude than the actual data. In Fig. A1 in Appendix A, the time series of the six RHALS factors, stacked on top of each other, are plotted against the total time series of the data matrix. We sometimes see larger magnitudes from the RHALS algorithm, along with seeing smaller magnitudes other times, but there seems to be no general pattern of bias towards solutions of greater or lesser magnitude generated by RHALS.

Finally, we determine if a solution with a different number of factors is more optimal with the RHALS algorithm for this dataset. Figure A6 shows the convergence of weighted residual error of RHALS as the number of factors increases, with the random initialization from the 13th trial. Error dramatically decreases as the factors increase from one to six, while barely improving for solutions with a larger number of factors. Therefore, a six-factor solution is justified by the RHALS algorithm.

4.4 Testing rotations

For rotations of entire factors, we first tested the rotational algorithm detailed in Paatero and Hopke (2009) to determine whether or not it was applicable to RHALS. Unfortunately the pulling equations used in the multilinear engine (ME) do not seem to work for RHALS as well.

Implementing the approach laid out in Sect. 2.2 to some of the solutions, we find that this method can potentially find better factorizations in the solution space. Ranging the values of a and b in Eqs. (28) and (29) between 0 and 500 (which is of a similar magnitude of the optimal amount of L2 regularization), we test the rotational method on the RHALS solutions from the 11th, 13th, and 17th trials. The graphs are plotted, in which a positive pulling value corresponds to a “pull-up” of W and a “pull-down” of H and vice versa for negative pulling values. The values of a and b are the magnitude of the pulling parameter. We present the rotation of the 13th solution in Figs. 14 and 15, with the rotations of the 11th solution in Figs. A7 and A8 and rotations of the 17th solution in Figs. A9 and A10 in Appendix A.

Figure 14 Weighted residual error of rotated solutions from the 13th trial of RHALS.

Figure 15 (a) W, (b) H, and (c) total average similarities of the rotated solutions to the PMF2 solution from the 13th trial of RHALS. The total similarity refers to an average of the similarity metrics for the time series and mass spectra factors.

As can be seen by analyzing the rotation of the solution from the three different trials, solutions that have lower weighted residual error and are closer to the target solution can be found through this simple pulling method, regardless of the direction of the pull. Interestingly, weighted residual error was able to be decreased to under 7000 when rotating the solution from the 13th trial towards larger values of W with a pulling parameter of 400. Note that the matrix norms do not always increase with larger positive values of the pulling parameter (and vice versa). This is most clearly seen in Fig. A11 in Appendix A with the rotated factors from the 11th trial. As the factors are pulled in one direction, the matrix norms may respond in the opposite direction.

The magnitude of change in the similarity to the PMF2 factors is small for all solutions, and rotations may only be worthwhile to look at when no good solution exists from the onset. The benefit of the rotations is more apparent when rotating poor solutions.

5 Conclusions

As the size of datasets has grown, computational costs have become increasingly expensive for traditional PMF algorithms. Thus, randomized and hierarchical algorithms are attractive alternative methods. Specifically, the RHALS algorithm was shown to provide a reduction in runtime compared to the multiplicative update algorithm, as well as the deterministic HALS algorithm. Furthermore, we proposed a novel approach to handling uncertainties in a weighted factorization problem. While this approach, coupled with randomization, slightly reduced the accuracy of the algorithm, it dramatically decreased the computational cost. Ultimately, we showed that our weighted RHALS algorithm was able to almost completely recreate the factors in both a formed test matrix and a real dataset and is a useful tool for finding nonnegative factors of large datasets, particularly in the context of real-time atmospheric mass spectrometry.

Appendix A: ALS derivation

The cost function Q is defined as

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}Q& =\left|\right|(\mathbf{\text{A}}-\mathbf{\text{W}}\mathbf{\text{H}})\oslash \mathbf{\Sigma }|{|}_F^{\mathrm{2}}+\sum _{i=\mathrm{1}}^m\mathit{\alpha }|\left|{\mathbf{\text{W}}}_{(i,:)}\right|{|}_{\mathrm{1}}\\ & +\sum _{j=\mathrm{1}}^n\mathit{\beta }\left|\right|{\mathbf{\text{H}}}_{(:,j)}|{|}_{\mathrm{1}}+\mathit{\gamma }|\left|\mathbf{\text{W}}\right|{|}_F^{\mathrm{2}}+\mathit{\delta }\left|\right|\mathbf{\text{H}}|{|}_F^{\mathrm{2}},\end{array}\end{array}$$

where $\mathit{\alpha },\mathit{\beta },\mathit{\gamma }$, and δ are regularization parameters. As with the derivation of HALS, the elementwise division of Σ is eliminated by considering a row or column of the residual at a time. Thus we find

$$\begin{array}{}\text{(A2)}& {\displaystyle }\begin{array}{rl}& \mathrm{0}\leftarrow {\displaystyle \frac{\partial Q_i}{\partial {\mathbf{\text{W}}}_{(i,:)}}}={\displaystyle \frac{\partial }{\partial {\mathbf{\text{W}}}_{(i,:)}}}\left(\right|\left|\right({\mathbf{\text{A}}}_{(i,:)}-{\mathbf{\text{W}}}_{(i,:)}\mathbf{\text{H}}\left){\mathbf{\Sigma }}_i^{-\mathrm{1}}\right|{|}_F^{\mathrm{2}}\\ & +\mathit{\alpha }\left|\right|{\mathbf{\text{W}}}_{(i,:)}|{|}_{\mathrm{1}}+\sum _{j=\mathrm{1}}^n\mathit{\beta }|\left|{\mathbf{\text{H}}}_{(:,j)}\right|{|}_{\mathrm{1}}\\ & +\mathit{\gamma }\left|\right|{\mathbf{\text{W}}}_{(i,:)}|{|}_{\mathrm{2}}^{\mathrm{2}}+\mathit{\delta }|\left|\mathbf{\text{H}}\right|{|}_F^{\mathrm{2}}),\end{array}\text{(A3)}& {\displaystyle }\begin{array}{rl}& \mathrm{0}\leftarrow {\displaystyle \frac{\partial Q_j}{\partial {\mathbf{\text{H}}}_{(:,j)}}}={\displaystyle \frac{\partial }{\partial {\mathbf{\text{H}}}_{(:,j)}}}\left(\right|\left|{\mathbf{\Sigma }}_j^{-\mathrm{1}}\right({\mathbf{\text{A}}}_{(:,j)}-\mathbf{\text{W}}{\mathbf{\text{H}}}_{(:,j)}\left)\right|{|}_F^{\mathrm{2}}\\ & +\sum _{i=\mathrm{1}}^m\mathit{\alpha }\left|\right|{\mathbf{\text{W}}}_{(i,:)}|{|}_{\mathrm{1}}+\mathit{\beta }|\left|{\mathbf{\text{H}}}_{(:,j)}\right|{|}_{\mathrm{1}}+\mathit{\gamma }\left|\right|\mathbf{\text{W}}|{|}_F^{\mathrm{2}}\\ & +\mathit{\delta }\left|\right|{\mathbf{\text{H}}}_{(:,j)}\left|{|}_{\mathrm{2}}^{\mathrm{2}}\right),\end{array}\end{array}$$

where Σ_i is a diagonal n×n matrix with the diagonal values equal to the ith row of Σ, and Σ_j is a diagonal m×m matrix with the diagonal values equal to the jth column of Σ. Using the fact that $\left|\right|\mathbf{\text{X}}|{|}_F^{\mathrm{2}}=\mathrm{Tr}({\mathbf{\text{X}}}^T\mathbf{\text{X}})$, with Tr being the trace of the matrix, Eqs. (A2) and (A3) can be rewritten as

$$\begin{array}{}\text{(A4)}& {\displaystyle }\begin{array}{rl}\mathrm{0}& ={\displaystyle \frac{\partial }{\partial {\mathbf{\text{W}}}_i}}\left(\mathrm{Tr}\right({\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{A}}}_i^T{\mathbf{\text{A}}}_i{\mathbf{\Sigma }}_i^{-\mathrm{1}}-\mathrm{2}{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{A}}}_i^T{\mathbf{\text{W}}}_i\mathbf{\text{H}}{\mathbf{\Sigma }}_i^{-\mathrm{1}}\\ & +{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{H}}}^T{\mathbf{\text{W}}}_i^T{\mathbf{\text{W}}}_i\mathbf{\text{H}}{\mathbf{\Sigma }}_i^{-\mathrm{1}})+\mathit{\alpha }|\left|{\mathbf{\text{W}}}_{(i,:)}\right|{|}_{\mathrm{1}}\\ & +\sum _{j=\mathrm{1}}^n\mathit{\beta }\left|\right|{\mathbf{\text{H}}}_{(:,j)}|{|}_{\mathrm{1}}+\mathit{\gamma }|\left|{\mathbf{\text{W}}}_{(i,:)}\right|{|}_{\mathrm{2}}^{\mathrm{2}}+\mathit{\delta }\left|\right|\mathbf{\text{H}}\left|{|}_F^{\mathrm{2}}\right),\end{array}\text{(A5)}& {\displaystyle }\begin{array}{rl}\mathrm{0}& ={\displaystyle \frac{\partial }{\partial {\mathbf{\text{H}}}_j}}\left(\mathrm{Tr}\right({\mathbf{\text{A}}}_j^T{\mathbf{\Sigma }}_j^{-\mathrm{1}}{\mathbf{\Sigma }}_j^{-\mathrm{1}}{\mathbf{\text{A}}}_j-\mathrm{2}{\mathbf{\text{H}}}_j^T{\mathbf{\text{W}}}^T{\mathbf{\Sigma }}_j^{-\mathrm{1}}{\mathbf{\Sigma }}_j^{-\mathrm{1}}{\mathbf{\text{A}}}_j\\ & +{\mathbf{\text{H}}}_j^T{\mathbf{\text{W}}}^T{\mathbf{\Sigma }}_j^{-\mathrm{1}}{\mathbf{\Sigma }}_j^{-\mathrm{1}}\mathbf{\text{W}}{\mathbf{\text{H}}}_j)+\sum _{i=\mathrm{1}}^m\mathit{\alpha }|\left|{\mathbf{\text{W}}}_{(i,:)}\right|{|}_{\mathrm{1}}\\ & +\mathit{\beta }\left|\right|{\mathbf{\text{H}}}_{(:,j)}|{|}_{\mathrm{1}}+\mathit{\gamma }|\left|\mathbf{\text{W}}\right|{|}_F^{\mathrm{2}}+\mathit{\delta }\left|\right|{\mathbf{\text{H}}}_{(:,j)}\left|{|}_{\mathrm{2}}^{\mathrm{2}}\right),\end{array}\end{array}$$

where ${\mathbf{\text{A}}}_{(i,:)}={\mathbf{\text{A}}}_i$, ${\mathbf{\text{A}}}_{(:,j)}={\mathbf{\text{A}}}_j$, ${\mathbf{\text{W}}}_{(i,:)}={\mathbf{\text{W}}}_i$, and ${\mathbf{\text{H}}}_{(:,j)}={\mathbf{\text{H}}}_j$. Using vector derivative rules and the fact that $\mathrm{Tr}\left(\mathbf{\text{A}}\mathbf{\text{B}}\mathbf{\text{C}}\right)=\mathrm{Tr}\left(\mathbf{\text{C}}\mathbf{\text{A}}\mathbf{\text{B}}\right)=\mathrm{Tr}\left(\mathbf{\text{B}}\mathbf{\text{C}}\mathbf{\text{A}}\right)$ when ABC is square, the following gradients are found.

$$\begin{array}{}\text{(A6)}& {\displaystyle }\begin{array}{rl}& \mathrm{0}\leftarrow -\mathrm{2}{\mathbf{\text{A}}}_i{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{H}}}^T+\mathrm{2}{\mathbf{\text{W}}}_i\mathbf{\text{H}}{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{H}}}^T\\ & +\mathit{\alpha }{\mathbf{1}}_{\mathrm{1}\times k}+\mathrm{2}\mathit{\gamma }{\mathbf{\text{W}}}_i\end{array}\text{(A7)}& {\displaystyle }\begin{array}{rl}& \mathrm{0}\leftarrow -\mathrm{2}{\mathbf{\text{W}}}^T{\mathbf{\Sigma }}_j^{-\mathrm{1}}{\mathbf{\Sigma }}_j^{-\mathrm{1}}{\mathbf{\text{A}}}_j+\mathrm{2}{\mathbf{\text{W}}}^T{\mathbf{\Sigma }}_j^{-\mathrm{1}}{\mathbf{\Sigma }}_j^{-\mathrm{1}}\mathbf{\text{W}}{\mathbf{\text{H}}}_j\\ & +\mathit{\beta }{\mathbf{1}}_{k\times \mathrm{1}}+\mathrm{2}\mathit{\delta }{\mathbf{\text{H}}}_j\end{array}\end{array}$$

In Eqs. (A6) and (A7), 1 is a vector of ones. Finally, the following update rules are found.

$$\begin{array}{}\text{(A8)}& {\displaystyle }\begin{array}{rl}& {\mathbf{\text{W}}}_i\leftarrow \left[\right({\mathbf{\text{A}}}_i{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{H}}}^T-{\displaystyle \frac{\mathit{\alpha }}{\mathrm{2}}}{\mathbf{1}}_{\mathrm{1}\times k})\\ & (\mathbf{\text{H}}{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\Sigma }}_i^{-\mathrm{1}}{\mathbf{\text{H}}}^T+\mathit{\gamma }{\mathbf{\text{I}}}_{k\times k}{)}^{-\mathrm{1}}{]}_{+}\end{array}\text{(A9)}& {\displaystyle }\begin{array}{rl}& {\mathbf{\text{H}}}_j\leftarrow \left[\right({\mathbf{\text{W}}}^T{\mathbf{\Sigma }}_j^{-\mathrm{1}}{\mathbf{\Sigma }}_j^{-\mathrm{1}}\mathbf{\text{W}}+\mathit{\delta }{\mathbf{\text{I}}}_{k\times k}{)}^{-\mathrm{1}}\\ & ({\mathbf{\text{W}}}^T{\mathbf{\Sigma }}_j^{-\mathrm{1}}{\mathbf{\Sigma }}_j^{-\mathrm{1}}{\mathbf{\text{A}}}_j-{\displaystyle \frac{\mathit{\beta }}{\mathrm{2}}}{\mathbf{1}}_{k\times \mathrm{1}}){]}_{+}\end{array}\end{array}$$

Here, I is the identity matrix.

Figure A1 RHALS time series factors laid out with total time series from data, plotted as a 2 h rolling average.

Figure A2 Comparison of 2 h rolling average of time series for externally weighted and internally weighted ALS factors. Externally weighted error: 7.3666×10³. Internally weighted error: 6.4768×10³.

Figure A3 Comparison of 2 h rolling average of time series for externally weighted and internally weighted MU factors. Externally weighted error: 6.9457×10³. Internally weighted error: 6.4292×10³.

Figure A4 Comparison of mass spectra for internally weighted (on the left) and externally weighted (on the right) ALS factors.

Figure A5 Comparison of mass spectra for internally weighted (on the left) and externally weighted (on the right) MU factors.

Figure A6 Number of factors versus weighted residual error for the RHALS algorithm.

Figure A7 Weighted residual error of rotated solutions from the 11th trial of RHALS.

Figure A8 (a) W, (b) H, and (c) total average similarities of the rotated solutions to the PMF2 solution from the 11th trial of RHALS.

Figure A9 Weighted residual error of rotated solutions from the 17th trial of RHALS.

Figure A10 (a) W, (b) H, and (c) total average similarities of the rotated solutions to the PMF2 solution from the 17th trial of RHALS.

Figure A11 Frobenius norms of rotated factor matrices for the (a–b) 11th, (c–d) 13th, and (e–f) 17th trials.

Appendix B: Mathematical equivalence of pseudo-inverse update and ordinary least squares

For simplicity, we label $\left(\stackrel{\mathrm{̃}}{\mathbf{\text{W}}}\stackrel{\mathrm{̃}}{\mathbf{\text{H}}}\right)\odot \mathbf{\Sigma }$ as $\stackrel{\mathrm{̃}}{\mathbf{\text{y}}}$. Let us consider the ordinary least squares update for W:

$$\mathbf{\text{W}}=\stackrel{\mathrm{̃}}{\mathbf{\text{y}}}{\mathbf{\text{H}}}^T(\mathbf{\text{H}}{\mathbf{\text{H}}}^T{)}^{-\mathrm{1}}.$$

H is a k×n matrix that is assumed to be of a full-row rank and can thus be decomposed into the exact SVD factorization H=USV^T, where U is a k×k orthogonal matrix, S is a k×k diagonal matrix, and V is an n×k matrix with orthonormal columns. The pseudo-inverse of H is defined as VS⁻¹U^T. V^TV=I and ${\mathbf{\text{U}}}^T={\mathbf{\text{U}}}^{-\mathrm{1}}$. The update for W can be rewritten as

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}\mathbf{\text{W}}& =\stackrel{\mathrm{̃}}{\mathbf{\text{y}}}\mathbf{\text{V}}\mathbf{\text{S}}{\mathbf{\text{U}}}^T(\mathbf{\text{U}}\mathbf{\text{S}}\mathbf{\text{I}}\mathbf{\text{S}}{\mathbf{\text{U}}}^T{)}^{-\mathrm{1}}\\ & =\stackrel{\mathrm{̃}}{\mathbf{\text{y}}}\mathbf{\text{V}}\mathbf{\text{S}}{\mathbf{\text{U}}}^T(\mathbf{\text{U}}{\mathbf{\text{S}}}^{\mathrm{2}}{\mathbf{\text{U}}}^T{)}^{-\mathrm{1}}\\ & =\stackrel{\mathrm{̃}}{\mathbf{\text{y}}}\mathbf{\text{V}}\mathbf{\text{S}}{\mathbf{\text{U}}}^T\left(\mathbf{\text{U}}{\mathbf{\text{S}}}^{-\mathrm{2}}{\mathbf{\text{U}}}^T\right)\\ & =\stackrel{\mathrm{̃}}{\mathbf{\text{y}}}\mathbf{\text{V}}\mathbf{\text{S}}\mathbf{\text{I}}{\mathbf{\text{S}}}^{-\mathrm{2}}{\mathbf{\text{U}}}^T\\ & =\stackrel{\mathrm{̃}}{\mathbf{\text{y}}}\left(\mathbf{\text{V}}{\mathbf{\text{S}}}^{-\mathrm{1}}{\mathbf{\text{U}}}^T\right)\\ & =\stackrel{\mathrm{̃}}{\mathbf{\text{y}}}{\mathbf{\text{H}}}^{†},\end{array}\end{array}$$

which is equal to the pseudo-inverse update described in Sect. 2.3. A similar argument can be made for the update for H.