An NMF-Based Building Block for Interpretable Neural Networks With Continual Learning

The main point of this work is to make some progress on developing learning methods with a better balance of interpretability and predictive performance compared to the existing MLP-based approaches. In particular, we make the following contributions:

• •

  In Section 2.4 we introduce a new neural network building block called the Predictive Factorized Coupling (PFC) block as a more interpretable alternative to the MLP. Since its declarative model is still matrix factorization, it potentially retains the interpretable parts-based nature of NMF, but extends it to support its use as a general differentiable predictive module.

• •

  In Section 5.2 we demonstrate that the PFC block has competitive accuracy with the (fully-connected) MLP on MNIST, Fashion MNIST, and CIFAR10 classification. In Section 5.7 we demonstrate that a (fully-connected) residual network consisting of two PFC blocks is also able to be trained without optimization difficulties and that it performs similarly. Section 5.6 show an example of the increased interpretability by visualizing in-domain vs out-of-domain input examples.

• •

  In Section 3 we develop a factorized RNN by starting from the well-known vanilla RNN and then replacing its MLP building block with our PFC blocks. This results in an RNN modeled as a single matrix factorization, effectively extending the parts-based nature of NMF to the modeling of sequential data. In Section 5.8 we demonstrate its interpretability advantages compared to the standard vanilla RNN on a simple sequential learning task with a known interpretable solution and observe that it is consistently able to learn the minimal transition model of the solution, even when the network is heavily over parameterized. In all of our sequence learning tasks we find that the factorized RNN performs either competitively or better compared with the vanilla RNN when both are trained with the usual BPTT. This includes learning a repeating sequence (Section 5.8), Copy Task (Section 5.9), Sequential MNIST (Section 5.10), and audio source separation (Section 5.11). We also perform some ablations and find two unexpected results of particular interest: The factorized RNN is able to solve simple tasks such as learning a repeating sequence and the Copy Task using only alternating NMF update rules so that backpropagation is not used at all. We also observe that in both the factorized and conventional RNNs, using backpropagation but disabling BPTT often results in a surprisingly minimal degradation in accuracy. The models do tend to become significantly less parameter efficient without BPTT, however.

• •

  In Section 5.3 we show that our non-replay-based continual learning method is competitive with approaches that rely on replay on the split MNIST Class-IL scenario [13]. Our PFC-based model performs better than the MLP on this task and we introduce a sliding window optimizer in Section 4 for PFC-based models that results in further accuracy improvements.

• •

  In Section 5.4 we show the superiority of a PFC-based model over the MLP on a non-i.i.d. training task.

• •

  In Section 5.5 we show how a PFC-based model can support knowledge removal after training by leveraging the sliding window optimizer that we introduce in Section 4.

We also mention the following limitations:

• •

  The PFC block is slower and consumes more memory during training than the MLP block due to the use of an iterative inference algorithm that effectively turns each replaced MLP into a corresponding RNN. It is therefore not intended to be a general MLP replacement, but rather it is intended to be used in modeling tasks where its better interpretability, parts-based modeling prior, and/or suitability to continual learning is needed.

• •

  Since the PFC block is modeled as a matrix factorization, it similarly cannot discriminate between two inputs that differ only by a scale factor, since the corresponding outputs would also only differ by the same scale factor. The NMF modeling assumption also constrains the input features to being non-negative, although this can potentially be relaxed to semi-NMF to support negative data values at the expensive of possibly reduced interpretabiity.

• •

  This work is preliminary. Our experiments only involve relatively small datasets. In terms of the suitability of the PFC block as an MLP replacement, we only consider two multi-block architectures in this paper: a 2-block residual network and a factorized RNN. Evaluating on larger datasets and/or more complex multi-block architectures is left as future research.

Non-negative matrix factorization (NMF) is a matrix decomposition method where a data matrix $V\in\mathbb{R}^{M\text{x}N}$ is factorized into two matrices $W\in\mathbb{R}^{M\text{x}R}$ and $H\in\mathbb{R}^{R\text{x}N}$, with the property that all three matrices have no negative elements. The objective is to find two non-negative matrices whose product closely approximates the initial matrix, such that:

Here, $R$ is a tunable hyperparameter that leads to a more compressed approximation as it is reduced. NMF was originally proposed by [14] as positive matrix factorization and later popularized by [15]. NMF’s strength lies in its ability to provide interpretable decompositions. In particular, the non-negativity constraint is found to lead to sparse and parts-based representations because only additive, not subtractive, combinations are allowed [15]. In many real-world applications, such as image representation or document analysis, this aligns well with the nature of the data where the measured quantities (like pixel intensity or word counts) are inherently non-negative.

It is also possible to relax the non-negativity constraint so that only the factor matrix $H$ is constrained to be non-negative, which is called semi-NMF [16]. This allows for more flexibility in representing data but can also reduce the interpretabiity.

The factorization is not unique and is typically computed using iterative algorithms, such as gradient descent or multiplicative updates, for example [17]. This involves first specifying a differentiable reconstruction loss function such as mean-squared error (MSE), for example, and then jointly optimizing $W$ and $H$ to minimize the loss. These algorithms alternately fix one of the factor matrices and update the other, until approximate convergence to a fixed point. In this paper we will use the same terminology from [18] to discriminate between these two updates, so that the NMF left update step refers to applying the update rule once to update the $W$ factor matrix, which we also refer to as the NMF learning update step. Likewise, the NMF right update step refers to applying the update rule once to update the $H$ factor matrix, which we also refer to as the NMF inference update step.

The usual modeling convention is that the columns of $V$ correspond to $N$ input feature vectors, which are $M$-dimensional. For example, the neural network interpretation of NMF shown in Figure 3 of [15] corresponds to:

where $v_{n}$ and $h_{n}$ are now individual column vectors (at columns index $n$) in $V$ and $H$, respectively. Here, $v_{n}$ represents the observed input features (visible variables), while $W$ represents the neural network weights and $h_{n}$ represent the hidden variables which are inferred from $v_{n}$ via the repeated application of the NMF right update rule until convergence. The columns of $W$ are often referred to as the learned basis vectors in the literature and we will also sometimes refer to them as (parts-based) prototypes in this paper.

We can see from the basic NMF formulation in Eq. 2.1 that it attempts to learn an approximation to a supplied data matrix. While this can be suitable for many tasks, it also tends to limit its use to unsupervised learning and tasks where NMF is able to provide sufficient modeling expressiveness. For these reasons, NMF is not typically used for supervised learning or for tasks requiring more expressiveness, such as modeling sequential data, for example.

The k-nearest-neighbors (k-NN) algorithm [10] is one of the simplest and most interpretable machine learning methods used for classification and regression. It is considered an instance or prototype-based method since the model stores the entire training dataset and then makes predictions using a similarity measure.

Suppose that we have a set of training example pairs $\{(x_{1},y_{1}),(x_{2},y_{2}),\ldots,(x_{N},y_{N})\}$ where $x_{n}$ is the input feature vector and $y_{n}$ is the corresponding target output value or vector. In the classification setting, $y_{n}$ would contain the ground truth class label and could be represented either as an integer label or a 1-hot vector. In the regression setting, $y_{n}$ could more generally be an arbitrary real-valued vector. Since the following formulation could apply to either, we will refer to it as nearest neighbor prediction.

In typical machine learning algorithms, the model adjusts its weights (i.e., parameters) to minimize a loss function that measures the difference between the model’s predictions and the actual target values. However, there is no such loss function in nearest neighbors. Rather, the learning procedure is extremely simple as it only consists of storing the training examples in a suitable data structure as they become available. We will still refer to this data structure as the “weights” since it contains the learned knowledge extracted from the training set, which in this case is just the training examples themselves. For easier connection to the matrix-factorization-based models that follow in later sections, we will formulate nearest neighbor prediction as a special case of matrix factorization. For each training example $(x_{n},y_{n})$, we create a corresponding column vector $w_{n}$ in the weights matrix $W$ as the vertical concatenation of the target $y_{n}$ on top of input feature $x_{n}$. Specifically, let $y_{n}\in\mathbb{R}^{C}$ (i.e., $C$ class labels) and let $x_{n}\in\mathbb{R}^{M}$. We then construct $w_{n}\in\mathbb{R}^{C+M}$ as:

As a result, $W$ will be a $(C+M)$ x $N$ matrix containing all of the training examples as its columns. We will find it useful to split $W$ into an upper “prediction” sub-matrix $W_{y}$ consisting of the first $C$ rows (containing the $y_{n}$ targets), and a lower “recognition” sub-matrix $W_{x}$ consisting of the last $M$ rows (containing the $x_{n}$ input features):

The nearest neighbor algorithm also requires us to define a suitable distance or similarity metric. For this, we can specify a similarity function $sim(x_{q},x_{k})$ which computes a similarity score between two supplied vectors based on their Euclidean distance or cosine similarity, for example.

Given a new input example $x$, we can use the model to perform inference and predict the output (either classification or regression) as follows. In the first step, we use $sim(x,x_{k})$ to compute the similarity score of $x$ against each of the columns $x_{k}$ in $W_{x}$. The general form of the algorithm finds the k-nearest neighbors (k-NN) but we consider only the 1-nearest neighbor (1-NN) method here for simplicity. Let $m$ denote the index of the best-matching column in $W_{x}$ (which contains the training example $x_{m}$). We refer to this step as “recognition” since the input $x$ was recognized as its nearest neighbor $x_{m}$, which represent the model’s reconstruction of $x$. The second step then involves selecting the same column of $W_{y}$ so that the model will output $y_{m}$ as its prediction. We therefore refer to this as the “prediction” step.

We can express the inference solution to the 1-NN in the form of a matrix-vector product that will make the reconstructive-predictive aspect of the model explicit. We introduce a 1-hot vector $h\in\mathbb{R}^{N}$ having the same dimensionality as the column count in $W$, where only index $m$ is set to 1. This allows us to express the model’s output prediction in the form of the following matrix-vector product:

We can likewise express the model’s input prediction as:

Using Eq. 2.4, these can be combined so that both the input and output predictions are given by a single product:

Note that given the (1-hot) inference solution in $h$, it picks out a single column in the weights to use as the predictions for both the output and input. For the more general k-NN algorithm, we could consider computing $h$ by setting all elements to 0 except those for indices corresponding to the k-nearest neighbors, for which we would use the corresponding value of the similarity function and then normalize these $k$ values to sum to 1 in $h$. Several interpretable aspects of the 1-NN prediction algorithm are worth mentioning. For example, it provides confidence estimation by making use of the similarity score. When the model makes a wrong prediction (or when an out-of-distribution example is supplied), we can inspect the reconstructed input $x_{pred}$ since it shows what the input was recognized as. Continual learning and knowledge removal are easily supported as well, since the training examples comprise the columns of $W$ and can therefore easily be added or removed as necessary.

A drawback of the 1-NN predictor is that the output $y_{pred}$ is not a differentiable function of $W$ and the input $x$. This prevents us from connecting it to a loss function such as classification loss in order to learn the weights using backpropagation, and also prevents its use as a building block for more complex architectures. The reliance on a manually specified similarity function can limit the prediction accuracy. In addition, inference can be expensive for large datasets since the number of columns in $W$ grows as the number of training examples. Despite these drawbacks, 1-NN and the more general k-NN can still sometimes produce a good balance of predictive performance and interpretability depending on the dataset. It is also worth mentioning that learnable versions of nearest neighbor prediction exist and are referred to as Learnable Vector Quantization (LVQ) [11] [12]. However, we are not aware of a fully differentiable version that would enable its use as a building block for more expressive architectures.

In contrast to nearest neighbor methods, which are non-parametric, most machine learning algorithms use a set of learnable parameters. In modern neural networks, also known as deep neural networks (DNNs), most of the learnable parameters tend to be contained in the MLP blocks as discussed in the Section 1. The versatility of DNNs lies in their ability to use backpropagation [19] [20] [21] for learning parameters that minimize a chosen loss function—ranging from classification to regression losses. This flexibility allows us to take a foundational component like the MLP and scale it to create architectures of varying complexity such as convolutional networks, deep residual networks, recurrent neural networks (RNNs), transformer architectures, etc.. It is this ability to optimize the parameters so as to minimize a desired loss function that appears to contribute to DNNs’ superior predictive performance over other more interpretable methods.

The basic MLP corresponds to the sequential connection of two affine transformations (linear layers) with a non-linear activation function in between. The hidden layer activation vector $h$ is expressed as the following (differentiable) function of the parameters and input vector $x$:

where the weights matrix $W_{xh}$ and bias vector $b_{h}$ are the parameters of the first linear layer and $\sigma()$ can be an arbitrary differentiable activation function. Common choices for $\sigma()$ include ReLU, GELU [22], and tanh, for example. The MLP’s predicted output, $y_{\text{pred}}$, is then obtained by applying a second affine transformation to $h$:

where the weights matrix $W_{hy}$ and bias vector $b_{y}$ are the parameters of the second linear layer. Note that $y_{\text{pred}}$ is differentiable with respect to the parameters of both linear layers, hidden activations, and input. This allows is to be used as a building block to create expressive neural architectures. Also note that unlike a single linear layer or Single Layer Perceptron (SLP), an MLP with at least one hidden layer can approximate complex non-linear functions [5] [6]. This ability comes from the non-linear activation function $\sigma()$ used in the hidden layer. As discussed in Section 1, a drawback of MLP-based models is that they are considered black-box models lacking interpretabiity and have challenges dealing with continual learning and training on non-i.i.d. examples.

Although the above presentation of MLPs might make them seem completely unrelated to nearest neighbor methods, we can make a connection between them. We do this by showing that a particular choice of weights initialization and activation function choice for the MLP makes it equivalent to the k-NN predictor. This is the reason why we overloaded $h$ to refer to both the MLP hidden activations in this section, while also using it to refer to the 1-hot vector (or k-nonzero in the case of k-NN) $h$ in Eqs 2.5 2.6 2.7. To see the connection, we first ignore the MLP bias terms. Instead of using the usual backpropagation procedure to learn the weights, we instead initialize MLP weights $W_{xh}$ and $W_{hy}$ to the corresponding nearest neighbor weights $W_{x}$ and $W_{y}$ that were used in Eqs 2.5 2.6 2.7 and also normalize the columns of $W_{xh}$ to have unit L2 norm. Recall that $W_{x}$ contained all of the training input vectors $x_{n}$ as its columns while $W_{y}$ contained the corresponding target output vectors $y_{n}$ as its columns, for a total of $N$ columns in each weight matrix. For the activation function $\sigma()$, we use a k-max activation followed by softmax, which has the effect of only passing the k-largest values and then normalizing them to have unit sum. For example, with $k=1$, the output $h$ of $\sigma()$ will be a 1-hot vector similar to the $h$ vector that we used for nearest neighbors. The final step is that we need to ensure that the input $x$ to the MLP is normalized to have unit L2 norm. Note that with this setup, the matrix product in Eq. 2.8 computes the cosine similarity between each of the training inputs in $W_{xh}$ and $x$, corresponding to this choice of similarity measure in nearest neighbors. The activation function then has non-zero outputs only for the indices corresponding to the k-nearest neighbors in $W_{xh}$, where these values must sum to 1 in $h$. As a result, we see from Eq. 2.9 that the MLP then outputs a linear combination of the corresponding $k$ columns in $W_{hy}$, (which contain the target training vectors $y_{n}$). Without using backpropagation, this “nearest-neighbor” MLP remains interpretable. However, if we attempt to train it, the interpretation of the weights as containing prototypes is potentially lost. For example, with the usual backpropagation training, there is no input reconstruction loss used by default, and so the input prediction $x_{pred}=W_{xh}h$ (analogous to Eq. 2.6) may not necessarily result in a good reconstruction.

In this section we introduce the Predictive Factorized Coupling (PFC) block and discuss how it relates to the models in the previous background sections. Having covered the related models, we can now easily introduce the PFC block simply by re-interpreting the matrix product shown in Eq. 2.7 with the NMF modeling interpretation of Eq. 2.2 instead of k-NN, so that it represents the predicted input and output vectors after using an NMF algorithm to solve for $h$:

Interpreted as NMF, it shows the prediction for a single column vector $v$ of data matrix $V$ in the factorization $V\approx WH$, where $V$ contains the input vectors $x_{n}$ to be recognized as the columns of its lower sub-matrix $X$ and the corresponding target output vectors ${y_{target}}_{n}$ to be predicted as the columns of its upper sub-matrix $Y_{targets}$:

When used as a neural building block, we consider $X$ to contain the input vectors to the PFC block. The block then infers (solves for) $H$ and predicts $Y_{pred}$ for its output. Similar to neural network training, the corresponding targets $Y_{targets}$ are not available during the prediction (i.e., inference) process. $V$ is therefore partially observed during inference, since only its $X$ sub-matrix is available to the block as input. $Y_{targets}$ is then only used for the purpose of computing the prediction loss when it is available.

Recall that under the NMF modeling constraints, the three matrices of $V\approx WH$ are only required to be non-negative (the non-negativity constraints for $V$ and $W$ can also potentially be removed if we allow semi-NMF, but we will assume NMF here for the purpose of describing the model). We see that $H$ is now less constrained compared to the nearest neighbor model that required $h$ to be 1-hot for 1-NN or contain only $k$ non-zero elements for the k-NN case. Since $W$ is learnable under NMF, we no longer need to use the same number of column vectors (which are often called basis vectors in the NMF literature) as training examples. Similar to Section 2.1, we let the hyperparameter $R$ specify the number of learnable basis vectors in $W$, and these can now be initialized to random non-negative values as an alternative to initializing with training examples. $R$ also specifies the dimension of the hidden vector $h$. By keeping $h$ internal to the block, we are free to later add or remove basis vectors from $W$ without changing the external interface of the block. With this interpretation, $h$ corresponds to the inferred hidden activations that represent an encoding of the input in terms of the (parts-based) basis vectors in $W$. From Eq. 2.10 we see that $W$ is also composed of two sub-matrices, $W_{y}$ and $W_{x}$, which represent the learned basis vectors as coupled “input-output” parts or prototypes. These could also be interpreted as learned key-value factors, where the input vector $x$ serves as the query, and the block is seen to perform a kind of factorized attention over parameters in predicting its output.

The factorization expression in Eq. 2.10 corresponding to the PFC block was first proposed in our earlier work [18] where we referred to it as a “coupling module” or “coupling factorization” and only considered its use in sequential data models. This contrasts with its more general usage as a building block for both sequential and non-sequential models in the current work. Refer to Section 6 for a more detailed discussion of related work.

To motivate the idea, we first need to review the vanilla RNN. In the following, we initially assume that the weights ($W$ matrices) have already been learned, so that we only need to consider the inference or forward pass through the network to compute its output predictions from its inputs. Given an input sequence of length $T$ containing the feature vectors $x_{0},x_{1},\dots,x_{T-1}$, the RNN maps them in order, one at a time, into a corresponding output sequence of vectors $y_{0},y_{1},\dots,y_{T-1}$. The subscript denoting the position in the sequence is often called the “time slice”, even when there is no notion of time involved. The reason they must be mapped one at a time is because the RNN maintains an internal hidden state vector $h_{k}$ which is consumed as an additional (hidden) input during each time slice, modified, and produced as a (hidden) output for use in the next time slice. So, in time slice $k$, the RNN will consume input $x_{k}$ and previous hidden state input $h_{k-1}$. It then produces a new hidden state $h_{k}$ and output $y_{k}$. The vanilla RNN does this in two stages. In the first stage, we first update the hidden state

where $\sigma()$ denotes an arbitrary activation function and $b_{h}$ is the bias vector. In the second stage, we compute the output from the updated hidden state:

Note that Eq. 3.1 can be rewritten as a single linear layer followed by nonlinear activation function:

where we let $W_{z}$ refer to the combined weights:

and let $z_{k}$ refer to the combined inputs:

Combining Eq. 3.2 and Eq. 3.3, we can express a single time slice of the RNN computation as:

This shows that each time slice $k$ of the RNN can be interpreted as an MLP that takes input $z_{k}$ and produces outputs $y_{k}$. From Eq. 3.5, we see that the previous hidden state $h_{k-1}$ appear together with $x_{k}$ in the input $z_{k}$, and the updated hidden state $h_{k}$ corresponds to the hidden layer of the MLP after the $\sigma()$ activation as shown in Eq. 3.3.

Section 3.1 showed that the each time slice of the vanilla RNN can be interpreted as an MLP. With this interpretation, it is now interesting to consider the model that results from replacing each of these MLP blocks with a corresponding PFC block. This corresponds to keeping the output linear layer in Eq. 3.2 (although with the bias term removed). We replace the input linear layer and activation function from Eq. 3.3 with the following vector factorization for the $k$’th time slice:

With this change, we have now reversed the direction of the linear mapping compared to the MLP so that the input $z_{k}$ is approximately a linear function of the hidden state $h_{k}$ (contrast this to the MLP in Eq. 3.3 where the pre-activation $h_{k}$ is a linear function of $z_{k}$). Our factorized representation is also simplified compared to Eq. 3.3 since we have removed the need for an activation function and bias term. The tradeoff is that now when we are given an input $z_{k}$, we will require an iterative NMF update algorithm to solve for $h_{k}$, which could be more computationally costly compared to the MLP.

With the computed $h_{k}$, we then compute the output as in the vanilla RNN using Eq. 3.2. Applying Eq. 3.2 and Eq. 3.7 to all time slices $k\in 0\dots T-1$, we finally arrive at the factorized vanilla RNN expressed as a single matrix factorization of the form $V\approx WH$:

Using Eq. 3.5, we can also express the left matrix $V$ in terms of $y_{k}$, $h_{k-1}$, and $x_{k}$. This brings us to the key result which lets us express the the factorized RNN with all of the model inputs, outputs, hidden states, and weights together in a single matrix factorization:

Note that for a single time slice, our model corresponds to the following vector factorization:

If we use $W$ to denote the three stacked weights sub-matrices, the notation simplifies even further to the following:

Now contrast the simplicity of the factorized RNN for a single time slice in Eq. 3.11 with the corresponding vanilla RNN expression for a single time slice in Eq. 3.6. Note that they differ in that Eq. 3.11 is a declarative representation while Eq. 3.6 is imperative. That is, for the factorized RNN we will still need to find suitable algorithms to actually solve the factorization, whereas the vanilla RNN expression explicitly tells us the steps needed to produce the output.

We can further simplify the notation by replacing each of the sequences in Eq. 3.9 with their respective sub-matrices. Letting $Y=\begin{bmatrix}y_{0}&y_{1}&y_{2}&\dots&y_{T-1}\end{bmatrix}$, $H_{prev}=\begin{bmatrix}h_{-1}&h_{0}&h_{1}&\dots&h_{T-2}\end{bmatrix}$, and $X=\begin{bmatrix}x_{0}&x_{1}&x_{2}&\dots&x_{T-1}\end{bmatrix}$ results in the following:

Regarding the hidden states, we see that the initial “previous” $h_{-1}$ only appears in the left $V$ matrix while the final state $h_{T-1}$ only appears in the right $H$ matrix. The other hidden states are duplicated since $h_{k}$ for $k\in[0\dots T-2]$ appear in both $H_{prev}$ and $H$. Also note that $W_{h}$ must be an $R$ x $R$ square sub-matrix of $W$ since if the $h_{k}$ are $R$-dimensional then each of $W_{y}$, $W_{h}$, and $W_{x}$ must also have $R$ columns.

A simple method of training the factorized RNN consists of performing alternating NMF updates to $W$ and $H$, while also copying the inferred states from $H$ to the corresponding duplicated positions in the data matrix after each inference update step to satisfy the consistency constraints. This is similar to the approach that we used for the sequential models in [18]. The details of this are as follows. We use the simple SGD-based algorithm which is reviewed in Appendix A for our experiments, but a variety algorithms can potentially be used.. Starting from the factorization in Eq. 3.9, we begin by initializing the weights $W={W_{y},W_{h},W_{x}}$ to small random values and initializing the hidden states to either zeros or small random values. This corresponds to setting sub-matrices $H_{prev}$ and $H$ to zeros or small random values in Eq. 3.12. Let $X=\begin{bmatrix}x_{0}&x_{1}&x_{2}&\dots&x_{T-1}\end{bmatrix}$ represent the training inputs and $Y=\begin{bmatrix}y_{0}&y_{1}&y_{2}&\dots&y_{T-1}\end{bmatrix}$ represent the corresponding target output values that we want to predict. Since we want to predict the outputs, using only the provided inputs, we first need to infer the hidden states. For that we use the subpart of Eq. 3.9 corresponding to Eq. 3.7:

The task is then to solve for the $h_{k}$, which we can do by alternating between matrix factorization updates of the right $H$ matrix followed by enforcing the constraint that the duplicated $h_{k}$ must have equal values. We can do this by simply copying the $h_{k}$ from the $H$ (right factor matrix) to $H_{prev}$ (left data matrix) after each NMF update to $H$. Once the updates converge, we can then use the top-most sub-factorization of Eq. 3.9 to compute the predicted outputs $\hat{y_{k}}$ as a linear function of the $h_{k}$:

Now that the inference part of the algorithm has complete, we perform the learning updates. We replace the predicted outputs with the target outputs above and perform a NMF update on $W_{y}$. Similarly, we perform an NMF update on $W_{h}$ and $W_{x}$ in Eq. 3.13. We then repeat the above procedure until convergence.

Since the hidden state updates propagate one slice forward for each NMF update, it will take the same number of updates as the sequence length for information from the first slice to potentially reach the last slice. As a result, if the sequence of long, the initial NMF updates late in the sequence could be considered wasted computation since they will be operating on hidden states that only contain information from nearby slices. Whether or not this actually becomes an issue in practice would seem to depend on how far information can actually propagate in an RNN, which is outside the scope of this paper.

As an alternative to performing the “full batch” updates as above, we can consider performing the inference “one time slice at a time”. Specifically, we start at the first slice $k=0$ and wait for the inference procedure to converge on that slice before containing with the next. Since the (output) inferred hidden state $h_{0}$ has now converged, we then copy it into the duplicated (input) location in the next ($k=1$) slice of $H_{prev}$. We can now increment the current slice to $k=1$ and continue in the same way so that the input $h_{k}$ states in the current slice of $V$ can always be considered to have already converged. This is similar to how inference is carried out in the vanilla RNN, since both RNNs share the same temporal dependency ordering between the hidden states. Although the slice-at-a-time inference option seems less able to take advantage of parallel hardware, it also seems potentially more efficient when the sequence length is extremely long since we perform the minimum number of iterations needed for convergence on each state before moving on to the next. Both options seem interesting but a detailed empirical comparison of their relative efficiency is outside the scope of this paper. Perhaps a batch-wise version could also be interesting as a topic for future research. We only consider the latter slice-at-a-time option in the experiments in Section 5.8.

The training algorithm in Section 3.3 used standard MF update rules to learn the model weights. These updates optimize the local reconstruction loss of the data matrix. Depending on the task, we empirically observed that such algorithms can sometimes be sufficient and we demonstrate example of this in the experiments Sections 5.8 5.9. However, we generally found this method to fail to perform well on more complex tasks. For this reason we use the algorithm unrolling approach as discussed in Section 2.4 for all other experiments.

It is straightforward to apply unrolling to the factorized RNN since the inference algorithm remains unchanged: we evaluate the computation graph and backpropagate through it. Since we replaced each MLP of the vanilla RNN with a corresponding PFC block, the computation graph dependencies result in the inference progressing one time slice at a time, similar to the vanilla RNN. However, note that in performing the inference in a given time slice, the PFC block’s unrolled NMF update steps form another RNN (internal to each PFC block) with length corresponding to the number of unrolled iterations. We therefore have a computation graph corresponding to an RNN within an RNN.

We compute the loss using the MSE loss between the targets and the predicted values. Since the factorized RNN has three predicted sub-matrices in $V$ (Eq. 3.12), this means we will have three loss terms. Once the inference procedure converges so that the inferred hidden states are available in $H$, we can predict $\hat{V}$ as follows:

The total loss is the the sum of the three loss terms with an arbitrary non-negative scale factor to adjust the relative strength of each term:

We then backpropagate through the loss to compute the error gradients and use an existing SGD-based optimizer such as RMSprop etc. to update the weights. Since the inference procedure operated one time slice at a time, we can see that the gradients then flow backward through the entire sequence, effectively making it a instance of backpropagation through time (BPTT).

We first evaluate a simple 1-hidden-layer MLP-based classifier as a baseline model on the MNIST [27], Fashion MNIST [28], and CIFAR10 [29] datasets. For each dataset, we train models for hidden dimensions sizes of 300, 2000, and 5000. The input feature size is equal to the number of image pixels when the image is flattened into a vector. This is $28x28=784$ for MNIST and Fashion MNIST since they contain 28x28 grayscale images and $32x32x3=3072$ for CIFAR10 since color images are used. The output layer dimension is 10 since all three datasets have 10 class labels. Table 1 shows the accuracy results.

We train a 1-block PFC-based network on the same datasets and with the same parameter sizes as the MLP from Section 5.1. Recall that the PFC basis vector count corresponds to the MLP hidden dimension. We also train with and without enforcing non-negative parameters (i.e., NMF vs semi-NMF) and also evaluate the effect of including vs disabling the input reconstruction loss term.

We train with early stopping after 20 epochs with no validation loss improvement. Table 2 shows the accuracy results, averaged over 3 training runs. We see that using semi-NMF generally leads to slightly better accuracy. The impact of input reconstruction loss on accuracy is less clear and seems to vary depending on the specific dataset and parameter configuration. Since both the NMF constraint and enabling reconstruction loss can potentially lead to better interpretability, we will enable the reconstruction loss in all remaining experiments. We will also use the NMF parameter constraint in the remaining experiments unless otherwise mentioned.

Comparing the PFC results in Table 2 with the MLP results in Table 1 shows the PFC-based models to perform competitively. We see that the PFC-based model performs significantly better compared to the MLP on CIFAR10, slightly better on Fashion MNIST, and similarly on MNIST, although the MLP does perform slightly better on MNIST for the case of 300-dimensional hidden dimension when the NMF constraint is used on the PFC-based network.

In this experiment we evaluate and compare the performance of the PFC and MLP-based models on the Split MNIST task [30] under the Class-IL scenario as described in [13]. In Split MNIST, the MNIST dataset is split into 5 task partitions, so that each task only contains two digits: Split 0 contains digits 0 and 1, Split 1 contains digits 2 and 3, and so on up to Split 4. The Class-IL scenario is the most difficult of the three considered in [13], as it requires the model solve the tasks that have appeared so far as well as inferring the task ID. Since the model is not told which of the 5 tasks it needs to solve, it only receives the image pixels as input and needs predict the correct digit label. The model will therefore have 28x28 = 784 inputs corresponding to the image pixels flattened into a vector and it will have 10 outputs corresponding to the digit labels. We train the model on one split at a time, ordered by the split number so that Split 0 will be the first task. Within each task the examples are presented i.i.d. to the model. The validation loss is computed over all tasks seen so far and early stopping is used to end the current task and move on to the next once the best validation loss is achieved.

Neural networks are most effectively trained when the stream of training examples is presented i.i.d. and the class labels are balanced. Training difficulties being to arise when distribution shifts are introduced. In this experiment, we train MLP and PFC-based models on the MNIST classification task. However, rather than the usual i.i.d. training in which the examples are presented in shuffled order, we instead present the examples in a fixed label-sorted order. Specifically, we sort the training examples ascending by class label so that all of the digit 0 examples are presented first, followed by all of the digit 1 examples, and so on, until finally presenting all of the digit 9 examples. We also train the networks with the usual i.i.d. (shuffled) ordering for comparison to provide an upper bound on the achievable accuracy.

With the PFC-based model, we also have the option of using the sliding learnable window optimizer from Section 4, which is implemented in RMSpropSLW. Since the examples are presented in the same order each epoch, we simply reset the learnable window to the starting position (i.e., the left-most column of each weight matrix) at the beginning of each epoch. This can potentially improve performance when training on non-i.i.d. data since optimizer updates for a particular example (or its batch) are constrained to a small learnable window of the weights. As a result, only nearby training batches (in terms of the number of training iterations between them) are capable of overwriting previous learned knowledge. More widely spaced batches cannot interfere with each other.

It is sometimes desireable to remove learned knowledge from a model. For example, if it is discovered that certain training batches contained errors, or if certain knowledge needs to be removed for legal reasons. A large model might take a long time to train, and so the ability to quickly remove specific knowledge could be preferable to retaining from scratch.

In this experiment, we train a network on the MNIST classification task in which some of the training examples have incorrect class labels. We show that this reduces the classification performance, as is expected. Provided that the bad examples/batches can be identified after training has completed, we show that it is possible to perform unlearning by removing the subset of model weights that were influenced by these bad examples during training, restoring much of the lost classification accuracy.

We use a model with 1 PFC block and 3000 available basis columns for each of the two weight matrices. For simplicity, we constrain the corrupted training examples so that they appear a contiguous range of batches. We use the same fixed shuffling of examples each epoch so that the i’th batch of each epoch always contains the same examples. There are then a total of 1020 batches per epoch (with 50 examples per batch). Batches 500 through 800 are corrupted by changing the class label to different (and incorrect) class. This is accomplished by incrementing the label modulo the number of classes.

We use the same ‘RMSpropOptimizerSlidingWindow‘ optimizer from the continual learning experiments. We reset the sliding window to the left-most position at the beginning of each epoch and use a deterministic dataset loader so that the i’th batch always contains the same examples in each epoch. This will cause the learning updates corresponding to the i’th batch to be stored within a knowable range of columns in the weights matrices, corresponding to the position of the sliding window while the said batch was active.

After the model is trained, 2600 of the 3000 basis columns are in use, as Figure 3(a) shows. The right-most 400 columns remain at their randomly initialized values, but this is not a problem and does not affect the results since these columns are not activated during the NMF inference process. Since the training data included a significant fraction of corrupted examples, the classification accuracy is a somewhat low 83.81% on the MNIST test set.

Next, we perform the unlearning operation, which will attempt to remove the knowledge obtained from the corrupted training examples from the network. This works as follows. First, recall that training batches 500 through 800 were identified as containing the corrupted labels. Since our SLW optimizer constrains the optimizer update for each batch to a small learnable window of the weights, we need to find the corresponding union of all positions of the learnable window during learning updates for these batches. We find that batch index 500 corresponds to the left-most column of the learnable window being at index 1250, and batch index 800 corresponds to the right-most column of the learnable window being at index 2050. That is, during this range of (corrupted example) batch updates, the sliding learnable window covered columns in the weight matrices ranging from 1250 through 2050, so that any knowledge learned from these batches must be in this subset of the weights. It is then straightforward to remove the knowledge, such as by deleting these weights, setting them to zero, or reinitializing to random values. For this experiment, we set these weights to zero, as shown in Figure 3(b). With the corrupted knowledge removed, we then re-evaluate the network on the test set and see that the accuracy has improved to 92.77%. For reference, if we train the model from scratch excluding the corrupted examples, we get 95.35% on the test set, which sets an upper bound on the unlearning accuracy. These results are summarized in Table 5. This shows that our method is effective in restoring accuracy without retraining, provided that we are able to identify which training batches were bad. Note that if any training examples are identified as bad, then the entire batch and corresponding region of the weights must be thrown out. This method is therefore most effective when the subset of batches to remove corresponds to a contiguous sequence of batches in the training data.

In this section, we train a network containing 1 PFC block on an image classification task and then evaluate it on both in-domain and out-of-domain (OOD) inputs. Since the network produces reconstructed input features during the recognition process, it is interesting to visualize and compare these reconstructions on in-domain vs OOD inputs. We might expect that the learned weights would correspond to the parts of the in-domain images, so that the reconstruction quality should be better on in-domain vs OOD inputs. This is because when the network is given an OOD image, it is forced to attempt to reconstruct it using on the “in domain” parts, potentially reducing the reconstruction quality compared to in-domain inputs.

We now empirically investigate these effects using MNIST as the in-domain images and Fashion MNIST as the OOD images. We train a small network on the MNIST classification task using 100 weight templates to enable the easy visualization of all weights in a single plot figure. We use the combined input reconstruction and classification loss as usual, except that here we use an adjustable trade-off between the two in the form of a hyperparameter $\lambda\in[0,1]$:

Figure 4 shows the weights, reshaped into images, for models trained using different values of $\lambda$. As $\lambda$ is increased, the strength of the classification loss increases and the reconstruction loss decreases. Thus, sub-figure 4(a) corresponds to a loss that emphasizes classification accuracy since only a small reconstruction loss is used. This is reflected in the good classification accuracy. We see that some of the weights resemble MNIST digits, but the images appear somewhat “noisy”. The middle sub-figure 4(b) shows the weights using the default blend of classification and reconstruction loss used in most of the other experiments. Here we see that the weights resemble MNIST digits or parts thereof. Finally, sub-Figure 4(c) shows the effect of a strong reconstruction loss. We see that the weights now appears as more localized image parts and that the classification accuracy is significantly lower.

We now compare the input reconstructions produce by the model for in-distribution vs OOD examples. For the following, keep $\lambda$ set to the default value of 0.5. For the in-distribution visualization, we train on MNIST and also evaluate on a batch of MNIST test images, as shown in Figure 5. Sub-figure 5(a) shows a batch of input MNIST test images and sub-figure 5(b) shows the corresponding input reconstructions generated by the model. We see that many of reconstructed images resemble the corresponding inputs, although some are less recognizable.

We now move on to the visualization of the model’s reconstructed images when given OOD inputs. For this we will use a batch of images from the Fashion MNIST test set, which contains grayscale images of fashion products of the same size as MNIST images. Recall that the model is constrained to reconstruct an input image by solving for the best additive combination of its weights (i.e., using the 100 images in 4(b)). However, since the available weights correspond to MNIST images and/or their parts, we might expect the reconstructed Fashion MNIST inputs to have less resemblance to the actual images. Indeed, this is what we observe in Figure 6, which shows the OOD inputs and their reconstructions.

In this experiment we use more than one PFC block to build a more complex architecture. The main purpose of this is to verify that if there are no optimization difficulties then the network should produce similar accuracy as that 1-block network in Section 5.2. We construct a simple fully-connected residual network containing 2 PFC blocks in which the first block uses a skip connection. That is, we let $x$ represent the input to the first PFC block, which then produces output $y_{1}$. The input to the second PFC block is then given by $x_{2}=relu(x-y_{1})$. The relu is optional when using the semi-NMF assumption but is needed when using the NMF constraint to prevent the input $x_{2}$ from becoming negative. The second PFC block then outputs the final prediction $y_{pred}$.

We train on the same datasets as in Section 5.2 using the same basis vector sizes (equivalent to MLP hidden dimension). Since there are two blocks, the total number of parameters is doubled from the 1-block model. Here we only train using the NMF modeling constraint (non-negative parameters) and we use both prediction and input reconstruction MSE loss terms for each PFC block. Table 6 shows the accuracy results on the test set. We see that the accuracy results appear in a similar range compared to those of the 1-block model in Table 2.

We developed the factorized RNN in Section 3.2 and modeled it as the matrix factorization of the form $V\approx WH$ shown in Eq. 3.9 which we repeat here:

We also provided a more compact notation in which the sequences are replaced by their respective sub-matrices in Eq. 3.12, which we repeat here:

We then provided two basic training procedures. Here we consider the first method which uses matrix factorization update rules for both inference and learning as described in Sections 3.3. It is initially unclear whether such simple update rules could even learn a task spanning several time slices, since there is no backward flow of error gradient-like information through the sequence. We therefore think it seems appropriate to start with a relatively simple temporal learning task.

We now test the factorized RNN and compare to the vanilla RNN on the more difficult Copy Task [32] using the setup from [33]. This involves generating a random sequence of tokens which are supplied one at a time in each time slice to the network. After this we supply a padding token for some fixed number of time slices. We then supply a “remember” token for a number of time slices equal to the length of input sequence that was supplied earlier, during which the network must attempt to output the remembered tokens of the input sequence. This task thus tests the network’s ability to recall information that was presented earlier and it gets more difficult as the task spans longer time intervals. We use the same parameters as [33] in which the input sequence is of length 10 and there are 10 distinct token values, which we supply to the network as 1-hot vectors. This combined with the pad and remember tokens lead to 12 distinct token values in total, so that the input 1-hot vectors will be 12-dimensional. We can then still adjust the difficulty by controlling the padding length $T_{pad}$.

The Sequential MNIST classification task is a common task used for evaluating the performance of RNNs. It adapts the MNIST dataset [27], which consists of 28x28 pixel grayscale images of handwritten digits (0 to 9), for a sequential data processing context. In the standard MNIST task, the entire image is presented to the model at once, as in Section 5.2. However, in the Sequential MNIST task, the image is presented as a sequence of pixels, typically in a row or column-wise manner. We evaluate the column-wise version in this experiment. The image is unrolled column by column, resulting in a 28-time-slice sequence, where each slice is a 28-dimensional vector representing one column of pixels. Each time slice of the RNN outputs a 10-dimensional class prediction vector, although we only compute the loss on the final time slice of the sequence, ignoring the model predictions from earlier time slices.

Audio tends to be well modeled as a mixture of source components. Much of the music we listen to is explicitly mixed together from isolated recordings (often called stems). In the source separation problem, we are interested in separating multiple sources that have been mixed together. For example, suppose we have a recording by a small band in which the vocals, guitar, and drums have all been captured together in a single audio channel. A source separation system would then be tasked with taking this recording as input and producing an output containing only a single desired source such as the vocals.

A useful inductive bias in this case would be to have the model $f()$ satisfy the linear superposition property: $alpha1*f(x_{1})+alpha2*f(x_{2})=f(alpha1*x_{1}+alpha2*x_{2})$. Any scaled mixture of the two sources should produce the corresponding scaled output. This is a desirable property since it implies that if the model is capable of recognizing the sources in isolation, it automatically generalizes to handling the mixture case. Such a property could potentially allow the model to generalize better from limited training data.

NMF can potentially satisfy this property, and there are several existing works in which it has been applied to related audio problems [34]. Since our factorized RNN is essentially a sequential extension of matrix factorization, it seems interesting to apply it to the source separation task and compare its generalization performance against a standard vanilla RNN (which does not have the superposition property due to its use of non-linear activation functions). We train both factorized and vanilla RNNs on the source separation task using the MUSDB18 dataset. MUSDB18 contains 150 music tracks of various genres corresponding to approximately 10 hours of music, split into a training dataset containing 100 tracks and a test dataset containing 50 tracks. It includes the isolated tracks (stems) for drums, bass, other, and vocals. This makes it straightforward to use for training and evaluating source separation models.

Perhaps the single most related existing work is our previous work on Positive Factor Networks [18], which similarly proposed NMF-based building blocks that can be composed to build expressive and interpretable architectures. Although the factorized equivalent of the vanilla RNN that we introduce in this paper was not considered, various other factorized-RNN-like architectures (which were referred to as dynamic positive factor networks) were considered, including some more sophisticated architectures such as a sequential data model for target tracking, and a model employing dilated deconvolutional layers in a Wavenet-like [35] architecture. A block with the same factorized model as our PFC block appears in Eq 3.1 of [18], which was referred to a “coupling module” or “coupling factorization” in that paper. Similar iterative NMF algorithms were used to perform inference in these models. However, a key distinction between that work and our present work is that it did not use backpropagation for learning the parameters, and instead relied on applying the NMF left-update rules to learn them. We also experimented with NMF left-update rules in the current work for learning the parameters as an ablation in Sections 5.8, 5.9, 5.10, but found them to underperform backpropagation on the more challenging learning tasks. As a result of not using backpropagation for parameter learning, the models presented in [18] failed to produce results competitive with other approaches on supervised learning tasks. The PFC block and corresponding neural networks constructed from them that we consider in the present work can therefore be considered as unrolled positive factor networks, or positive factor networks employing backpropagation-based training.

NMF was originally proposed by [14] as positive matrix factorization and later popularized by [15] [17]. NMF is related to other methods such as sparse coding [36] that use a similar dictionary learning model but with differences in the modeling constraints, regularizers, and/or algorithms used. NMF is often observed to provide parts-based decompositions of data, which can be useful when interpretable decompositions are desired. It also potentially satisfies the (non-negative) linear superposition property and as a result has been applied to audio processing tasks such as source separation and music transcription in which the audio features in the input data matrix are assumed to be well modeled as an additive combination of audio sources (i.e., individual instruments and/or notes) [34].

Our PFC block is related in that its declarative model corresponds to a masked predictive NMF in which the data matrix is partitioned along the row dimension into “input” and “output” vectors of the block, with the right factor matrix $H$ corresponding to inferred hidden activations. The output vectors as masked during recognition (during inference of $H$) and the resulting inferred $H$ is then used to predict the outputs. This allows our block to retain the additive superpositional NMF model, while also supporting the construction of more expressive differentiable architectures such as factorized RNNs. The resulting neural networks (or positive factor networks) can then be interpreted as an extension of NMF that increases its modeling expressiveness and suitability for supervised learning tasks.

The key enabler of the increased predictive performance of our present work compared to the Positive Factor Networks of [18] is the modification of the learning algorithm to use backpropagation instead of relying on NMF left-update steps. As discussed in Section 2.4, backpropagation-based training is enabled by unrolling the iterative NMF inference steps into an RNN-like structure in the computation graph, making the PFC block differentiable and therefore compatible with backpropagation training.

This process of unrolling an iterative optimization algorithm into a neural network and training with backpropagation is referred to as algorithm unrolling or unrolled neural networks in the literature [25]. An existing example of unrolled NMF appears in [37]. The idea that improved results could be obtained on supervised tasks by adding an additional task-specific loss function to an iterative and differentiable optimization algorithm and using it to optimize the parameters was proposed in [38] and [39], with earlier related ideas being proposed in [40] and [41]. More recent works have explored the use of unrolled convolutional sparse coding for improved robustness to corrupt and/or adversarial input images [42], [43]. Of these, [44], [43] also use FISTA to accelerate the convergence of the unrolled inference. As far as we are aware, ours is the first work to explore the use of a modular unrolled block (PFC block) supporting the construction of arbitrary neural architectures, unrolled factorized RNNs, and the first to explore the interpretability advantages of unrolled NMF-based networks for continual learning.

As discussed in Section 2.2, our PFC block shares some similarities with classification and regression based on the k-nearest neighbors (k-NN) algorithm [10]. We show that the PFC block can be derived by first formulating k-NN prediction as a matrix-vector product, and then generalizing and interpreting it as a matrix factorization.

Learning Vector Quantization (LVQ) [11] [12] is a prototype-based method that extends the k-nearest neighbors algorithm by introducing the concept of learnable prototypes. An advantage of this approach over matrix-factorization-based methods is its faster recognition since an iterative algorithm is not used. We are not aware of a fully differentiable version of LVQ that could match MLP predictive performance and make it possible to use as a building block for more complex architectures such as RNNs, however.

This work is preliminary and we leave several unanswered questions and possibilities for future research directions. We list some of them here:

• •

  As discussed in [45], methods such as ours that perform recognition by iterating to a fixed point can potentially make use of adaptive computation. This could potentially be used to improve recognition efficiency, as well as providing an additional form of confidence estimation, since “difficult” inputs could require more iterations to converge.

• •

  When unrolling the NMF inference updates, we observed that memory usage can potentially be reduced by computing the initial several iterations “without gradients” and only unrolling the last several iterations “with gradients”. In some cases, this resulted in improved efficiency, but in others it resulted in reduced accuracy, and so a better understanding is still needed.

• •

  Although we found FISTA to be one effective method in accelerating the NMF inference step of the PFC block, it could be interesting to also consider other approaches to further accelerate the inference.

• •

  Our sliding learnable window optimizer introduced in Section 4 enabled improved continual learning in PFC-based models. However, more sophisticated approaches based on basis-vector-specific learning rates could potentially lead to improved continual learning performance and/or parameter efficiency. For example, it could be interesting to consider methods that identify the important or in-use basis vectors and reduce their learning rates accordingly so as to make them more resistant to being overwritten or modified later in training.

• •

  We did not consider any form of sparsity regularization (e.g., L1 penalty) in our experiments, leaving regularization methods other than the non-negativity constraint and weight decay unexplored in the current work.

• •

  As Table 7 shows, disabling BPTT in both factorized and vanilla RNNs often resulted in only a small drop in accuracy, although more parameters were needed to recover accuracy. It is also unclear why the semi-NMF parameter constraint resulted in better accuracy compared to NMF when BPTT was disabled. It could be interesting to conduct a more detailed investigation as future research.

• •

  The declarative model of the PFC block in Eq. 2.10 is symmetric with respect to the inputs and outputs. This allows us to potentially reverse the prediction direction of the block so that we can consider swapping their roles and make the block compute the “inputs” given the “outputs”. In addition to reversible blocks/layers, it could be interesting to explore the case where only some subset of the input and output are observed and the block then jointly predicts all unobserved elements.

• •

  As discussed in Section 2.4, the PFC block can be interpreted as performing factorized attention over parameters. It could be interesting to also consider its use as an attention block that performs factorized attention over key and value activations (output from other upstream blocks) instead of over parameters, making it more like a factorized version of the attention block used in the transformer. Specifically, in Eq. 2.10 the $W_{x}$ matrix would be replaced by a keys matrix $K_{x}$ and $W_{y}$ would be replaced with a corresponding values matrix $V_{y}$ which would then represent output activations from other layers/blocks instead of learnable parameters.

• •

  Note that the PFC block has the inductive bias of NMF, which is different than the MLP. In our experiments, we observed it to perform similar to and sometimes better compared to the MLP. However, it is possible that the NMF inductive bias could make PFC-based models particularly well suited or poorly suited depending on the modeling assumptions of the datasets used. Also note that relaxing the non-negativity constraint to semi-NMF could support the use of datasets containing negative values, at the expense of possibly reduced interpretabiity, however. We leave such an exploration as future research.