# Alleviating Sequence Information Loss with Data Overlapping and Prime Batch
Sizes

Noémien Kocher  Christian Scuito  Lorenzo Tarantino  Alexandros Lazaridis
Andreas Fischer  Claudiu Musat

###### Abstract

In sequence modeling tasks the token order matters, but this information can
be partially lost due to the discretization of the sequence into data points.
In this paper, we study the imbalance between the way certain token pairs are
included in data points and others are not. We denote this a token order
imbalance (TOI) and we link the partial sequence information loss to a
diminished performance of the system as a whole, both in text and speech
processing tasks. We then provide a mechanism to leverage the full token order
information---Alleviated TOI---by iteratively overlapping the token
composition of data points. For recurrent networks, we use prime numbers for
the batch size to avoid redundancies when building batches from overlapped
data points. The proposed method achieved state of the art performance in both
text and speech related tasks.

##  1 Introduction

Modeling sequences is a necessity. From time series Connor et al. (1994); Lane
and Brodley (1999) to text Sutskever et al. (2011) and voice Robinson (1994);
Vinyals et al. (2012), ordered sequences account for a large part of the data
we process and learn from. The data are discretized and become, in this
paradigm, a list of tokens.

The key to processing these token sequences is to model the interactions
between them. Traditionally Rosenfeld (2000) this has been achieved with
statistical methods, like N-grams.

With the advances in computing power and the rebirth of neural networks, the
dominant paradigm has become the use of recurrent neural networks (RNNs)
Mikolov et al. (2010).

The dominance of RNNs has been recently challenged with great success by self-
attention based models Vaswani et al. (2017). Instead of modeling the sequence
linearly, Transformer-based models use learned correlations within the input
to weight each element of the input sequence based on their relevance for the
given task.

Series discretization. Both RNNs and self-attention models take as input data
points\---token sequences of a maximum predefined length---and then create
outputs for each of them. These tend to be much shorter in size, compared to
the size of the full dataset. While for humans time seems to pass
continuously, this discretization step is important for the machine
understanding of the sequence.

Figure 1: The common way of building data points given a dataset of contiguous
tokens. Here we illustrate a dataset with a contiguous list of 13 tokens, from
which we build 3 data points of 4 tokens each. This process keeps the order of
the tokens inside the data points, but loses the order information from token
pairs that happen to fall between adjacent data points.

A side effect of this step is a partial loss of the token order information.
As portrayed in Figure 1, we notice that the token order information within a
data point are kept. On the other hand, the knowledge about the token order at
the boundaries of data points is lost. We name the situation Token Order
Imbalance (TOI).

As the discretization in Figure 1 is the current standard of sequence
processing, we denote this as standard Token Order Imbalance (TOI). We
hypothesize that this loss of information unnecessarily affects the output of
the neural network models.

Alleviated Token Order Imbalance. A first contribution in this work is a
mechanism to ensure that all token sequences are taken into account, i.e.
every token pair is included in a data point and does not always fall between
two data point boundaries. Thus, all sequence information is available for
subsequent processing. The proposed method, denoted Alleviated TOI, employs a
token offset in the data point creation to create overlapped data point
sequences in order to achieve this effect.

Batch Creation with Alleviated TOI. A second contribution is a strategy for
batch creation when using the proposed Alleviated TOI method. We have observed
an unintended data redundancy within batches introduced by the overlapped data
point sequences. A strategy for avoiding this data redundancy is surprisingly
simple but effective: Always use a prime number for the batch size. The
intuition behind the prime batch size is that it ensures a good distribution
of the batches over the entire dataset. If used naively, the Alleviated TOI
policy leads to very similar data points being selected in a batch, which
hinders learning. By decoupling the batch size and the token offset used in
the token creation, this negative effect is effectively removed.

We then compare the Alleviated TOI with the Standard TOI and show that, on the
same dataset and with the same computation allocated, the Alleviated TOI
yields better results. The novel TOI reduction method is applicable to a
multitude of sequence modeling tasks. We show its benefits in both text and
voice processing. We employ several basic and state of the art RNNs as well as
Transformers and the results are consistent---the additional information
provided by the Alleviated TOI improves the final results in the studied
tasks.

For text processing we focus on a well-studied task---language modeling---
where capturing the sequence information is crucial. Using Alleviated TOI (P)
with the Maximum Over Softmax (MoS) technique on top of a recurrent cell Yang
et al. (2017) we get the new state of the art on the Penn-Tree-Bank dataset
without fine-tuning with 54.58 perplexity on the test set. We also obtain
results comparable to the state of the art on speech emotion recognition on
the IEMOCAP Busso et al. (2008) dataset111To make our results reproducible,
all relevant source codes are publicly available at
https://github.com/nkcr/overlap-ml.

The paper continues with an overview of the related work in Section 2, a
description of the alleviated TOI mechanism in Section 3 and a detailed
description of the batch generation in Section 4. The experimental design
follows in Section 5 and the results are detailed and interpreted in Section
6.

##  2 Related work

At the core of our work is the idea that the way that data samples are
provided for training a model can affect speed or capabilities of the model.
This field is broad and there are several distinct approaches to achieve it.
Notable examples include curriculum learning Bengio et al. (2009) and self-
paced learning Kumar et al. (2010), where data points for training are
selected based on a metric of easiness or hardness. In Bayesian approaches
Klein et al. (2016), the goal is to create sub-samples of data points, whose
traits can be extrapolated as the full dataset.

Our work thus differs from the aforementioned methods in the fact that we
focus on exploiting valuable but overlooked information from sequences of
tokens. We change the way data points are generated from token sequences and
extend the expressivity of a model by providing an augmented, and well sorted,
sequence of data points. This method has a related effect of a randomized-
length backpropagation through time (BPTT) Merity et al. (2017), which yields
different data points between epochs. It also resembles classical text data-
augmentation methods, such as data-augmentation using thesaurus Zhang and
LeCun (2015).

Our method takes a step forward and proposes a systematic and deterministic
approach on building data points that provides the needed variety of data
points without the need of randomized-length backpropagation through time
(BPTT). This has the effect of producing a text-augmentation without the need
of using external resources such as a thesaurus, but only requires the dataset
itself. Our method uses a concept of overlapped data points, which can be
found in many areas such as data-mining Dong and Pei (2007), DNA sequencing Ng
(2017), spectral analysis Ding et al. (2000), or temporal data Lane and
Brodley (1999). In language modeling however, this approach of overlapped data
points has not yet been fully exploited. On the other hand, extracting frame-
based acoustic features such as mel-fequency cepstral coefficients (MFCCs)
using overlapping windows is a common technique in speech processing and more
specifically in automatic speech recognition (ASR) Chiu et al. (2018); Xiong
et al. (2016); Kim and Stern (2016). We hypothesize that extending the current
overlapping technique to a higher level, that is using a sliding overlapping
window over the already extracted features, will be proven beneficial. We
believe this to have a positive impact on speech processing tasks such as
speech emotion recognition (SER). This is because the emotional load in an
spoken utterance expands over larger windows than frame-, phoneme- or
syllable-based ones Frijda (1986).

We investigate the proposed method using a simple LSTM model and a small-size
Transformer model on the IEMOCAP dataset Busso et al. (2008), composed of five
acted sessions, for a four-class emotions classification and we compare to the
state of the art Mirsamadi et al. (2017) model, a local attention based
BiLSTM. Ramet et al. (2018) showed in their work a new model that is
competitive to the one previously cited, following a cross-valiadation
evaluation schema. For a fair comparison, in this paper we focus on a non-
cross-valiation schema and thus compare our results to the work of Mirsamadi
et al. (2017), where a similar schema is followed using as evaluation set the
fifth session of IEMOCAP database. It is noteworthy that with a much simpler
method than presented in Ramet et al. (2018), we achieve comparable results,
underscoring the importance of the proposed method for this task as well.

##  3 Alleviated Token Order Imbalance

Let a token pair denote an ordered pair of tokens---for instance token A
followed by token B, as in the sequence $"ABCDEFG..."$. When splitting a token
sequence into data points $"D1,D2,.."$, if the split is fixed, as in $D1$
always being equal to $"ABC"$, $D2$ always being equal to $"DEF"$, etc., then
the information contained in the order of tokens $C$ and $D$ for instance is
partially lost. This occurs as there is no data point that contains this token
pair explicitly. We call the $"CD"$ token pair a split token pair and its
tokens, $C$ and $D$, are denoted as split tokens.

In its most extreme form, split token pair order information is lost
completely. In other cases, it is partially taken into account implicitly. In
recurrent cells, for instance, the internal state of the cell allows for the
order information of split tokens pairs to be used. This is due to the serial
processing of the data points containing the split tokens.

As some token pairs are taken into account fully, others partially and others
not at all, we denote this situation as token order imbalance (TOI).

In this paper, we propose to alleviate the TOI by means of overlapping
sequences of data points. The aim is to avoid the loss of information between
the last token of a data point and the first token of its subsequent data
point. Instead of splitting the sequence of tokens only once, we repeat this
process multiple times using different offsets. Each time we subdivide the
sequence of tokens with a new offset, we include the links that were missing
in the previous step. Finally, the overlapping sequences of data points are
concatenated into a single sequence, forming the final dataset.

Figure 2: Illustration of an Alleviated TOI (3) made from a single contiguous
list of 13 tokens. With a Standard TOI and $N$=3 (ie. 3 tokens per data
point), a contiguous list of 13 tokens would produce 4 data points, which is
illustrated by the first overlapped sequence. Here, an Alleviated TOI (3)
splits the contiguous list of tokens 3 times with each time a different offset
($0,1,2$ respectively). This finally leads to a list of 11 data points coming
from the 3 appended overlapped sequences.

Figure 2 illustrates an Alleviated TOI (3), which means the sequence of data
points is split three times instead of only once, producing 3 overlapped
sequences that will then be concatenated.

Our Alleviated TOI (P) method is detailed in the pseudo-code below, where
`olp_sequence` holds an overlapped sequence and $P$ is the number of times we
subdivide the sequence of tokens with a different offset:

    Let N = Number of tokens per data point
        P = Number of overlapped sequences
        Step = N / P
    DataPoints = empty list
    FOR i = 0..P-1
      olp_sequence = create data points
                     from Dataset with
                     offset (i * Step)
      Add olp_sequence to DataPoints
    RETURN DataPoints

When we apply an Alleviated TOI (P), this means that we are going to create
$P$ times a sequence of data points with different offsets. Therefore, the
final dataset will be the concatenation of $P$ repetitions of the original
dataset, with data points shifted by a specific and increasing offset at token
level for each repetition.

For example, given a sequence $S_{1}$ with $N=70$ tokens per data point and an
Alleviated TOI (P) with $P=10$, the step size will be $\frac{N}{P}=7$ tokens.
Therefore, starting from the sequence $S_{1}$, nine additional sequences of
data points will be created: $S_{2}$ starting from token 7, $S_{3}$ starting
from token 14, $S_{4}$ starting from token 21 and so on until $S_{10}$.

When using Alleviated TOI (P), with $P$ smaller than the data point size,
within an epoch, a split token pair---that is a token pair that is split in
the original data point splitting---becomes part of a data point $P-1$ times.
A token pair that is never split will be part of the data point $P$ times.

We can thus define a token order imbalance ratio that describes the imbalance
between the number of times we include split token pairs and the number of
times we include pairs that are not split:

$(P-1)/P$

We notice that the higher $P$, the closer the ratio becomes to 1. We
hypothesize that the closer the ratio becomes to 1, the better we leverage the
information in the dataset. We thus expect that for higher values of P the
Alleviated TOI (P) method will outperform versions with lower values, with
Alleviated TOI (1) being the Standard TOI, which is now prevalent.

We quantify the additional computational cost of Alleviated TOI (P). Since our
method only results in $P$ (shifted) repetitions of the dataset, each epoch
using the augmented dataset would take $\sim P$ times longer than an epoch
over the original dataset. Therefore, we ensure fair comparison by allowing
baseline models to train for $P$ times more epochs than a model using
Alleviated TOI (P).

##  4 Batch Creation with Alleviated TOI

Series discretization may also occur at higher levels than data points, in
particular when building batches for mini-batch training of neural networks.
We can distinguish two types of batches, i.e. sequential and distributed
batches. The former keep the data point sequences intact, thus creating split
token pairs only between two consecutive batches. The latter distribute data
points from different parts of the dataset to approximate the global
distribution, thus creating split token pairs between all data points in
batches.

In principle, our proposed method alleviates the TOI in both cases, since
multiple overlapping sequences of data points are generated. However, we have
observed an unintended interference with the batch creation in the case of
distributed batches. In this section we explain the problem in detail and
propose a simple but effective solution---choosing a prime batch size.

Figure 3 illustrates the three levels of data representation in the case of
distributed batches. Data points are built from $N$ consecutive tokens to
capture the sequential information. Batches are then built from $K$ parts of
the data point sequence to capture the global distribution. An example of this
approach is the batching procedure used in Zoph and Le (2016); Merity et al.
(2017); Yang et al. (2017); Zołna et al. (2017) for word language modeling,
where the basic token is a word.

Figure 3: Three levels of data representation used to create distributed
batches. The dataset is a sequence of tokens on which data points are built by
splitting the sequence into subsequences of N tokens. Batches of K data points
are then built by subdividing the sequence of data points into K equal parts.
Here, the first part contains the first two data points, the second part the
following two, and the last data point is dropped. Each batch then uses one
element of each part.

The batching mechanism can be seen as building a 2-dimensional matrix, where
each row contains a batch. Consider a sequence of $M$ data points and a batch
size of $K$. In order to build batches, the data points are split into $K$
parts, represented as $\frac{M}{K}\times 1$ column vectors. They are
concatenated to form a $\frac{M}{K}\times K$ matrix, such that the rows
correspond to batches.

(a) Matrix of batches with batch size of 20  (b) Matrix of batches with batch
size of 19

Figure 4: Illustrations of the 2D matrix of batches with different $P$-values of Alleviated TOI (P). On the left we used a batch size of 20 and on the right we used a prime batch size of 19. Each data point is a pixel and each row is a batch. The grayscale value models the proximity of the data points with respect to the dataset. Therefore, two pixels with similar color represents two data points that are close in the dataset. The illustrations demonstrate how different values of $P$ affect the content of the batches, which can lack a good distribution over the dataset. Ideally, each row should contain a gradient of different grayscale values. We can observe how using a prime batch size affects the distribution of data points within the batches, where the matrices on the right offer a better distribution. This effect is especially well visible for the Alleviated TOI 10.  Figure 5: Data point repetition with period $q$ for $M$ data points, $K$ batches, and Alleviated TOI (P). Data point 1’ is the same as data point 1 with a token offset. $P$ | Period $q$ | Repetitions
---|---|---
2 | 10 | 2
5 | 4 | 5
7 | 20 | 1
10 | 2 | 10
Table 1: Data point repetition with period $q$ for batch size $K=20$ and
Alleviated TOI (P).

When applying the proposed Alleviated TOI (P) method (see Section 3), we
augment the original dataset to a total of $P\cdot M$ data points, adding
additional data points with token offsets. Therefore, the $\frac{P\cdot
M}{K}\times K$ matrix used for batch creation may contain repeated data points
within the same batch as illustrated in Figure 5. A repeated data point
differs from the previous data point only marginally due to the token offset.
This redundancy can be problematic, as the batches are not well-distributed
over the entire dataset anymore.

With respect to the batch matrix, a repeated data point occurs iff
$\frac{P\cdot M}{K}\cdot q=n\cdot M$ with period $q<K$ and $q,n\in\mathbb{N}$.
This is equivalent to

$\frac{P}{K}\cdot q=n,\ \ \ q<K,q,n\in\mathbb{N}$

independent of the number of data points $M$. A repetition thus occurs iff the
greatest common divisor (GCD) of $P$ and $K$ is larger than 1. Otherwise, for
GCD$(P,K)=1$ a data point repeats only after period $q=K$, i.e. there is no
repetition within the same batch.

Table 1 lists exemplary periods for a batch size of $K=20$ and different
values of $P$ for the Alleviated TOI (P). The worst case is $P=10$ with $10$
repetitions of the same data point within the same batch and the best case is
$P=7$, which avoids any redundancy because the GCD of $P$ and $K$ is 1. Figure
4 illustrates the repetition with grayscale values, where similar grayscale
values indicate that two data points are close within the original data points
sequence.

In general, while we aim for large values of $P$ for reducing the TOI, a
simple solution for avoiding redundancy within batches is to choose a prime
number for the batch size $K$.

##  5 Experimental Setup

To validate the generalization capability of the proposed technique, we apply
it on both text and speech related tasks. We thus run the Alleviated TOI (P)
with language modeling (text) and emotion recognition (speech). The text
datasets used are Penn-Tree-Bank (PTB) Marcus et al. (1993) as preprocessed in
Mikolov et al. (2011), Wikitext-2 (WT2), and Wikitext-103 (WT103) Merity et
al. (2016). The speech dataset is the IEMOCAP database Busso et al. (2008), a
collection of more than 12 hours of recorded emotional speech of 10 native-
English speakers, men and women. The audio data is filtered down to 5.5 hours
containing only angry, happy, neutral and sad utterances.

###  5.1 TOI in Language Modelling

For language modeling, we use three different methods:

  * •

A simple LSTM that does not benefit from extensive hyper-parameter
optimization.

  * •

An Average Stochastic Gradient Descent Weight-Dropped LSTM (AWD-LSTM) as
described in Merity et al. (2017), with the same hyper-parameters.

  * •

The latest State-of-the-Art model: Mixture of Softmaxes (MoS) Yang et al.
(2017).

We compare our results against the original process of building data points,
i.e. Standard TOI, and use the same computation load allocated for each
experiment. We use the same set of hyper-parameters as described in the base
papers, except for the batch size with Alleviated TOI (P), where we use a
prime batch size in order to prevent any repetitions in batches, as described
in Section 4. That is, on the PTB dataset, we use a sequence length of 70 for
all the models. For the Simple LSTM and AWD-LSTM, we use a batch size of 20
and a hidden size of 400. AWD-LSTM and MoS are trained on 1000 epochs, and the
Simple LSTM on 100 epochs. For the MoS model, embedding size used is 280,
batch size 12, and hidden size 980. All the models use SGD as the optimizer.

We set up experiments to compare 4 different token order imbalance setups:
Extreme TOI, Inter-batch TOI, Standard TOI, and Alleviated TOI (P).

#### Extreme TOI

The Extreme TOI setup builds batches using a random sequence of data points.
This removes any order inside the batches (i.e. among data points within a
batch), and among batches.

#### Inter-batch TOI

In the Inter-batch TOI setup, batches are built using an ordered sequence of
data points, but the sequence of batches is shuffled. This keeps the order
inside batches, but removes it among batches. Looking at the 2D matrix of
batches, in Figure 4, this results in shuffling the rows of the matrix.

#### Standard TOI

In the Standard TOI setup, the process of building batches is untouched, as
described in section 3. This keeps the order inside, and among batches.

#### Alleviated TOI (P)

In the Alleviated TOI (P) setup, we apply our proposed TOI reduction by
creating $P$ overlapped data point sequences (see Sections 3 and 4). This
strategy not only keeps the order inside and among batches, but it also
restores the full token order information in the dataset.

###  5.2 TOI in Speech Emotion Recognition

For Speech Emotion Recognition (SER) we use two different models: the encoder
of the Transformer Vaswani et al. (2017) followed by convolutional layers, and
the simple LSTM used in text domain case. Since the Transformer is stateless
and uses self-attention instead, we are able to investigate the effect of
Alleviated TOI (P) independently of LSTM cells.

As with language modeling, we set up experiments to compare the 4 different
token order imbalance strategies: Extreme TOI, Inter-batch TOI, Standard TOI,
and Alleviated TOI (P).

We apply the methodology used in text on the SER task, using the simple LSTM
and a window size of 300 frames. In this case, a data point, instead of being
a sequence of words, is a sequence of frames coming from the same utterance.
Each frame is described by a 384-dimensional features vector. OpenSMILE Eyben
et al. (2013) is used for extracting the features. We opt for the IS09
features set Schuller et al. (2009) as proposed by Ramet et al. (2018) and
commonly used for SER.

Finally, to investigate the effect of the Alleviated TOI (P) strategy
independently of LSTM cells, we design a final experiment in the SER task. We
investigate whether or not we have improved results as we increase $P$, the
number of overlapped data point sequences in a stateless scenario. For this
reason, we use the Transformer model described above.

##  6 Experimental Results

###  6.1 Language Modelling

Table 2 compares the 4 token order imbalance strategies using the AWD model
and three text datasets. We use the test perplexity after the same equivalent
number of epochs. The different Alleviated TOI (P) experiments use a different
number of overlapped sequence: An Alleviated TOI (P) means building and
concatenating $P$ overlapped sequences. Our results indicate that an
Alleviated TOI (P) is better than the Standard TOI, which is better than an
Extreme or Inter-batch TOI. We note a tendency that higher values of $P$ lead
to better results, which is in accordance with our hypothesis that a higher
TOI ratio $(P-1)/P$ improves the results.

Experiment | PTB | WT2 | WT103
---|---|---|---
Extreme TOI | 63.49 | 73.52 | 36.19
Inter-batch TOI | 64.20 | 72.61 | 36.39
Standard TOI | 58.94 | 65.86 | 32.94
Alleviated TOI 2 | 57.97 | 65.14 | 32.98
Alleviated TOI 5 | 57.14 | 65.11 | 33.07
Alleviated TOI 7 | 57.16 | 64.79 | 32.89
Alleviated TOI 10 | 56.46 | 64.73 | 32.85
Table 2:  Perplexity score (PPL) comparison of the AWD model, on the three
datasets, with batch sizes $K=20$ (PTB), $K=80$ (WT2) and $K=60$ (WT103), with
different levels of Token Order Imbalance (TOI). With Alleviated TOI (P), we
use a prime batch size of $K=19$ (PTB), $K=79$ (WT2) and $K=59$ (WT103).

#### Comparison with State of the Art and Simple LSTM.

With the MoS model and an Alleviated TOI, we improve the current state of the
art without fine tuning for the PTB dataset with 54.58 perplexity on the test
set. Table 3 demonstrates how models can be improved by applying our
Alleviated TOI method on 2 latest state-of-the-art models: AWD-LSTM Merity et
al. (2017) and AWD-LSTM-MoS Yang et al. (2017), and the Simple LSTM model. We
compare the results with the same hyper-parameters used on the original papers
with the only exception of the batch size, that must be prime. To ensure
fairness, we allocate the same computational resources for the base model as
well the model with Alleviated TOI, i.e. we train with the equivalent number
of epochs.

Model | test ppl
---|---
AWD-LSTM Merity et al. (2017) | 58.8
AWD-LSTM + Alleviated TOI | 56.46
AWD-LSTM-MoS Yang et al. (2017) | 55.97
AWD-LSTM-MoS + Alleviated TOI | 54.58
Simple-LSTM | 75.36
Simple-LSTM + Alleviated TOI | 74.44
Table 3:  Comparison between state-of-the-art models Merity et al. (2017);
Yang et al. (2017) and a Simple LSTM, and the same models with Alleviated TOI.
The comparison highlights how the addition of Alleviated TOI is able to
improve state-of-the-art models, as well as a simple model that does not
benefit from extensive hyper-parameter optimization.

#### Comparison without prime batch size.

In Table 4 we demonstrate how using a prime batch size with Alleviated TOI (P)
actually impacts the scores. We compare the scores of a prime batch size
$K=19$ with the scores of the original batch size $K=20$ for the AWD model
with Alleviated TOI (P). When using a prime batch size, we observe consistent
and increasing results as $P$ increases. This is due to the good distribution
of data points in the batches regardless of the value of $P$, which is visible
in Figure 4(b) where each row contains a high diversity of grayscale values.
With the original batch size $K=20$, we observe a strong performance for
$P=7$, but a low performance for $P=10$. Again, this effect is related to the
distribution of data points in the batches, which is visible in Figure 4(a).
The matrix with $P=7$ shows a good distribution---corresponding to the strong
performance---and the matrix with $P=10$ shows that each row contains a low
diversity of data points.

Experiment | K=20 | K=19
---|---|---
Alleviated TOI 2 | 59.37 | 57.97
Alleviated TOI 5 | 60.50 | 57.14
Alleviated TOI 7 | 56.70 | 57.16
Alleviated TOI 10 | 65.88 | 56.46
Table 4:  Perplexity score (PPL) comparison on the PTB dataset and the AWD
model. We use two different values for the batch size $K$ \--- the original
one with $K=20$, and a prime one with $K=19$. The results directly corroborate
the observation portrayed in Figure 4, where the obtained score is related to
the diversity of grayscale values in each row.

###  6.2 Speech Emotion Recognition Results

The results on the IEMOCAP database are evaluated in terms of weighted (WA)
and unweighted accuracy (UA). The first metric is the accuracy on the entire
evaluation dataset, while the second is the average of the accuracies of each
class of the evaluation set. UA is often used when the database is unbalanced,
which is true in our case, since the happy class has a total duration that is
half of the second smallest class in speech duration.

Table 5 shows that our proposed method brings value in the speech related task
as well. When choosing the Extreme TOI instead of the Standard TOI approach we
observe a smaller effect than in text related task: this is due to the
different nature of the text datasets (large "continuous" corpuses) and the
IEMOCAP one (composed of shorter utterances). The fact that we can still
observe improvements on a dataset with short utterances is a proof of the
robustness of the method.

Experiment | WA | UA
---|---|---
Extreme TOI (15k steps) | 0.475 | 0.377
Inter-batch TOI (15k steps) | 0.478 | 0.386
Standard TOI (15k steps) | 0.486 | 0.404
Alleviated TOI (15k steps) | 0.553 | 0.489
Alleviated TOI (60 epochs) | 0.591 | 0.523
Table 5:  Token order imbalance (TOI) comparison for the IEMOCAP dataset on a
SER task using angry, happy, neutral and sad classes with a simple LSTM model.

A greater effect is obtained when we increase the size of the dataset with the
proposed Alleviated TOI (P) approach: Due to the increasing offset at each
overlapped sequence, the data fed into the model contains utterances where the
emotions are expressed in slightly different ways. For this reason, the
performance notably increases.

Table 6 reports the result of a final experiment that aims to investigate the
effect of Alleviated TOI (P) independently of LSTM cells. For each Alleviated
TOI (P) setup and Standard TOI described in Table 6, we repeat the training
and evaluation for each of the following window sizes: 100, 200, 300, 400 and
500 frames. The previously described Transformer model is used in these
experiments. The results reported in Table 6 are the mean $\pm$ the standard
deviation computed for different P-values of Alleviated TOI (P).

Experiment | WA (60 epochs) | UA (60 epochs)
---|---|---
Alleviated TOI 1 |  0.591$\pm$0.012 |  0.543$\pm$0.021
Alleviated TOI 2 |  0.594$\pm$0.007 |  0.549$\pm$0.016
Alleviated TOI 3 |  0.605$\pm$0.018 |  0.563$\pm$0.024
Alleviated TOI 5 |  0.608$\pm$0.015 |  0.562$\pm$0.028
Alleviated TOI 10 | 0.617$\pm$0.015 | 0.571$\pm$0.024
Local attention | 0.635 | 0.588

Table 6:  Token order imbalance (TOI) comparison for the IEMOCAP dataset on a
SER task using angry, happy, neutral and sad classes for 60 epochs using the
Transformer model.

The last line of Table 6 refers to Mirsamadi et al. (2017) results. We want to
highlight the fact that the goal of these experiments is to show the direct
contribution of the Alleviated TOI technique for a different model. For this
reason we use a smaller version of the Transformer in order to reduce the
computational cost. We believe that with a more expressive model and more
repetitions, the proposed method may further improve the results.

The results from Table 6 demonstrate that, as we increase the value of $P$,
more significant improvements are achieved. This is in accordance with our
hypothesis that a higher TOI ratio $(P-1)/P$ improves the results.

##  7 Conclusions

In this work, the importance of overlapping and token order in sequence
modelling tasks were investigated. Series discretization is an essential step
in machine learning processes which nonetheless can be responsible for the
loss of the continuation of the tokens, through the token order imbalance
(TOI) phenomenon. The proposed method, Alleviated TOI, has managed to overcome
this drawback and ensures that all token sequences are taken into account. The
proposed method was validated in sequence modelling tasks both in the text and
speech domain outperforming the state of the art techniques.

## References

  * Bengio et al. (2009) Yoshua Bengio, Jérôme Louradour, Ronan Collobert, and Jason Weston. 2009\.  Curriculum learning.  In _Proceedings of the 26th annual international conference on machine learning_ , pages 41--48. ACM.
  * Busso et al. (2008) Carlos Busso, Murtaza Bulut, Chi-Chun Lee, Abe Kazemzadeh, Emily Mower Provost, Samuel Kim, Jeannette N. Chang, Sungbok Lee, and Shrikanth Narayanan. 2008.  Iemocap: interactive emotional dyadic motion capture database.  _Language Resources and Evaluation_ , 42:335--359.
  * Chiu et al. (2018) Chung-Cheng Chiu, Tara N Sainath, Yonghui Wu, Rohit Prabhavalkar, Patrick Nguyen, Zhifeng Chen, Anjuli Kannan, Ron J Weiss, Kanishka Rao, Ekaterina Gonina, et al. 2018.  State-of-the-art speech recognition with sequence-to-sequence models.  In _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 4774--4778. IEEE.
  * Connor et al. (1994) Jerome T Connor, R Douglas Martin, and Les E Atlas. 1994.  Recurrent neural networks and robust time series prediction.  _IEEE transactions on neural networks_ , 5(2):240--254.
  * Ding et al. (2000) Mingzhou Ding, Steven L Bressler, Weiming Yang, and Hualou Liang. 2000.  Short-window spectral analysis of cortical event-related potentials by adaptive multivariate autoregressive modeling: data preprocessing, model validation, and variability assessment.  _Biological cybernetics_ , 83(1):35--45.
  * Dong and Pei (2007) Guozhu Dong and Jian Pei. 2007.  _Sequence data mining_ , volume 33.  Springer Science & Business Media.
  * Eyben et al. (2013) Florian Eyben, Felix Weninger, Florian Gross, and Björn Schuller. 2013.  Recent developments in opensmile, the munich open-source multimedia feature extractor.  In _Proceedings of the 21st ACM International Conference on Multimedia_ , MM ’13, pages 835--838, New York, NY, USA. ACM.
  * Frijda (1986) Nico H. Frijda. 1986.  _The Emotions_.  Cambridge University Press.
  * Kim and Stern (2016) Chanwoo Kim and Richard M Stern. 2016.  Power-normalized cepstral coefficients (pncc) for robust speech recognition.  _IEEE/ACM Transactions on Audio, Speech and Language Processing (TASLP)_ , 24(7):1315--1329.
  * Klein et al. (2016) Aaron Klein, Stefan Falkner, Simon Bartels, Philipp Hennig, and Frank Hutter. 2016\.  Fast bayesian optimization of machine learning hyperparameters on large datasets.  _CoRR_ , abs/1605.07079.
  * Kumar et al. (2010) M Pawan Kumar, Benjamin Packer, and Daphne Koller. 2010.  Self-paced learning for latent variable models.  In _Advances in Neural Information Processing Systems_ , pages 1189--1197.
  * Lane and Brodley (1999) Terran Lane and Carla E Brodley. 1999.  Temporal sequence learning and data reduction for anomaly detection.  _ACM Transactions on Information and System Security (TISSEC)_ , 2(3):295--331.
  * Marcus et al. (1993) Mitchell P. Marcus, Mary Ann Marcinkiewicz, and Beatrice Santorini. 1993.  Building a large annotated corpus of english: The penn treebank.  _Comput. Linguist._ , 19(2):313--330.
  * Merity et al. (2017) Stephen Merity, Nitish Shirish Keskar, and Richard Socher. 2017.  Regularizing and optimizing LSTM language models.  _CoRR_ , abs/1708.02182.
  * Merity et al. (2016) Stephen Merity, Caiming Xiong, James Bradbury, and Richard Socher. 2016.  Pointer sentinel mixture models.  _CoRR_ , abs/1609.07843.
  * Mikolov et al. (2011) Tomáš Mikolov, Anoop Deoras, Stefan Kombrink, Lukáš Burget, and Jan Černockỳ. 2011.  Empirical evaluation and combination of advanced language modeling techniques.  In _Twelfth Annual Conference of the International Speech Communication Association_.
  * Mikolov et al. (2010) Tomáš Mikolov, Martin Karafiát, Lukáš Burget, Jan Černockỳ, and Sanjeev Khudanpur. 2010.  Recurrent neural network based language model.  In _Eleventh annual conference of the international speech communication association_.
  * Mirsamadi et al. (2017) Seyedmahdad Mirsamadi, Emad Barsoum, and Cha Zhang. 2017.  Automatic speech emotion recognition using recurrent neural networks with local attention.  _2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 2227--2231.
  * Ng (2017) Patrick Ng. 2017.  dna2vec: Consistent vector representations of variable-length k-mers.  _arXiv preprint arXiv:1701.06279_.
  * Ramet et al. (2018) Gaetan Ramet, Philip N. Garner, Michael Baeriswyl, and Alexandros Lazaridis. 2018\.  Context-aware attention mechanism for speech emotion recognition.  _2018 IEEE Spoken Language Technology Workshop (SLT)_ , pages 126--131.
  * Robinson (1994) Anthony J. Robinson. 1994.  An application of recurrent nets to phone probability estimation.  _Trans. Neur. Netw._ , 5(2):298--305.
  * Rosenfeld (2000) Ronald Rosenfeld. 2000.  Two decades of statistical language modeling: Where do we go from here?  _Proceedings of the IEEE_ , 88(8):1270--1278.
  * Schuller et al. (2009) Björn Schuller, Stefan Steidl, and Anton Batliner. 2009.  The interspeech 2009 emotion challenge.  In _Tenth Annual Conference of the International Speech Communication Association_.
  * Sutskever et al. (2011) Ilya Sutskever, James Martens, and Geoffrey E Hinton. 2011.  Generating text with recurrent neural networks.  In _Proceedings of the 28th International Conference on Machine Learning (ICML-11)_ , pages 1017--1024.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  In _Advances in Neural Information Processing Systems_ , pages 5998--6008.
  * Vinyals et al. (2012) Oriol Vinyals, Suman V. Ravuri, and Daniel Povey. 2012.  Revisiting recurrent neural networks for robust asr.  In _2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 4085--4088.
  * Xiong et al. (2016) Wayne Xiong, Jasha Droppo, Xuedong Huang, Frank Seide, Mike Seltzer, Andreas Stolcke, Dong Yu, and Geoffrey Zweig. 2016.  Achieving human parity in conversational speech recognition.  _arXiv preprint arXiv:1610.05256_.
  * Yang et al. (2017) Zhilin Yang, Zihang Dai, Ruslan Salakhutdinov, and William W. Cohen. 2017.  Breaking the softmax bottleneck: A high-rank RNN language model.  _CoRR_ , abs/1711.03953.
  * Zhang and LeCun (2015) Xiang Zhang and Yann LeCun. 2015.  Text understanding from scratch.  _arXiv preprint arXiv:1502.01710_.
  * Zołna et al. (2017) Konrad Zołna, Devansh Arpit, Dendi Suhubdy, and Yoshua Bengio. 2017.  Fraternal dropout.  _stat_ , 1050:31.
  * Zoph and Le (2016) Barret Zoph and Quoc V Le. 2016.  Neural architecture search with reinforcement learning.  _arXiv preprint arXiv:1611.01578_.


# Identifying Population Movements with Non-Negative Matrix Factorization from
Wi-Fi User Counts in Smart & Connected Cities

Michael C. Huffman, Armen Davis, Joshua Park, James Curry

###### Abstract

Non-Negative Matrix Factorization (NMF) is a valuable matrix factorization
technique which produces a “parts-based” decomposition of data sets. Wi-Fi
usercounts are a privacy-preserving indicator of population movements in smart
and connected urban environments. In this paper, we apply NMF with a novel
matrix embedding to Wi-Fi usercount data from the University of Colorado at
Boulder Campus for the purpose of automatically identifying patterns of human
movement in a Smart & Connected infrastructure environment.

###### Index Terms:

Wi-Fi, NMF, Non Negative Matrix Factorization, Unsupervised, Smart Cities,
Smart & Connected Cities, IEEE

##  I Introduction

###  I-A Urbanization

The 20th century saw a monumental transition of human populations from rural
settings into urban ones. Beginning in Europe and the United States in the
19th century, industrialization spread across the world and drove an enormous
migration of people into city centers, forever altering the way we live. In
1950, 30% of the world population lived in an urban setting—today it is over
half [1]. Rapidly developing nations saw even more dramatic urban growth over
a comparatively shorter duration.

The 21st century will see more of the same. By 2030, the global population is
expected to pass 8.5 billion, on its way to 9.8 billion by 2050. Driven in
equal measures by this overall global population growth and the persistent
migration into urban areas, urban centers will continue to swell. India, China
and Nigeria alone are projected to add 416 million, 255 million and 189
million people respectively to their cities by 2050. By then, it is expected
that Earth’s urban corridors will have added 2.5 billion people and that 7 out
of 10 people on this planet will live in a city. [1]

Today, many livability hazards characteristic of urban population centers stem
from poor planning and little consideration to scalability: crime[2],
congestion[3], pollution[4], and inequity in access to services[5] are common
features of urban life and without careful consideration, will only become
worse.

###  I-B Smart Cities

Ensuring that city residents live safe, healthy, equitable and productive
lives is the central focus of the Smart & Connected Cities movement[6]. This
is a multidisciplinary push to combine big data analytics, civil engineering
and city planning in order to make urban living more human-friendly. Many
identified multiple areas of city life that could be improved through the
union of big data and city planning: decentralized energy production and
distribution[7], integrated business/economic models that anticipate a demand
for goods and services[8], optimized transportation minimizing pollution and
congestion[9], and real-time warnings with enhanced emergency-response[10]. At
its core, a central premise of Smart Cities is an intelligent response to
accurate and real-time mobility data.

The promises of Smart Cities to improve the human condition can only be
realized through widespread data collection and sophisticated understanding of
the movements of human population through their environments. This is a
critical and active area of research which capitalizes on the growth of
personal connected devices that are so common today. A 2016 study of
population movement in the industrial cities of northern Italy demonstrated
the insights possible from telecommunications data [11]. Another 2015 study of
telecommunications data demonstrates the power to reveal patterns of human
mobility [12]. However, such approaches have significant drawbacks.
Telecommunication data may only be available to carriers and not the
municipalities. Above all, such surveillance techniques are considered
intrusive and come with major privacy concerns[13].

###  I-C Wi-Fi

Wi-Fi describes a family of wireless communication protocols which operate
over unlicensed radio frequencies and allow devices to be networked in small
scale environments (homes, offices, and shops). Introduced in 1999 and
governed by the family of standards IEEE 802.11[14], today Wi-Fi access is a
common feature of urban environments. With the growing popularity of public
hotspots, Wi-Fi expected to see exponential growth in the coming years.
According to a 2018 industry white paper by telecommunications giant Cisco
Systems, Inc., public Wi-Fi hot-spots are expected to increase to 628 million
globally by 2023[15]–a fourfold increase over 5 years.

Wi-Fi has been previously identified as a valuable analogue for population
data. Researchers have shown that Wi-Fi can produce high fidelity population
counts of buildings[16], even more accurately than the HVAC industry standard
of CO2 sensors[17]. Building on this success, Wi-Fi has been identified as a
method to monitor and forecast crowds in large-scale public events[18]. Other
teams have validated Wi-Fi as a method of assessing city behavior with real-
world data[19] and there’s even been work done to develop a real-time city
census from Wi-Fi data[20].

In addition to its accurate representation of human populations, by only
collecting active user counts at Wi-Fi access points the technology can be
rendered anonymous and deployed as a privacy-preserving mode of monitoring
population movement. Because of these attractive features–accuracy, ubiquity,
growth, and anonymity–it is the belief of the authors that Wi-Fi is a safe
tool to monitor population movements and will play a critical role in
realizing the lofty goals of Smart Cities.

###  I-D Matrix Decomposition of Data

Unsupervised, data-driven approaches to automatically extract a low number of
latent components from complex data are of high demand in signal separation
problems such as single-channel audio source separation or hyperspectral
unmixing. Indeed, there is a similar problem present in the Wi-Fi user counts
of Smart & Connected infrastructure. The Smart & Connected Cities movement
hinges on our ability to understand and predict human movement through an
urban environment and respond intelligently. With that in mind, a natural
question is how can bulk movements of people and their devices into and out of
buildings be identified from single-channel signal of population data?

The canonical way to do such a reduction is the singular value decomposition
(SVD). With similar intentions, we stray from this path in favor of the Non-
negative Matrix Factorization (NMF). Unlike SVD, NMF yields strictly non-
negative components of the data. While both factorizations seek similar
results in the form of, for example, dimensionality reduction [24], NMF is
thought to return more interpretable results since many physical quantities
are constrained to be positive or non-negative (e.g. atmospheric pressure,
user counts). For these data, standard SVD will, in general, return
nonphysical components of the original data matrix which leads to
interpretability issues. Due to the non-negativity constraint imposed by NMF,
the components can be thought of as additive parts used to reconstruct the
data [23]. Since Wi-Fi user counts are strictly non-negative, NMF is certainly
a compelling alternative to SVD.

NMF has a rich history dating back to the early 1960s in the fields of
analytical chemistry and remote sensing [24]. The work done in the 60s was not
directly referred to as NMF, although it was equivalent. In fact, NMF did not
appear in its present form until 1994 when a paper was published on positive
matrix factorization [25]. The true explosion of the topic can be directly
traced to the Nature article by Lee and Seung in 1999 when data compression
and feature extraction were shown to be intrinsic properties of NMF [23].

##  II Methods and Materials

###  II-A University of Colorado Boulder Wi-Fi Data

Data for this study was provided by the Office of Information Technologies
(OIT). We generate total active user counts on a per-building basis by
consolidating multiple Wi-Fi access points (APs) that exist in each building.
The access points of the school are designed to be used in a mesh-like mode
where many access points provide coverage for a large area, while devices ebb
and flow between them. As a consequence, a device is only ever counted by one
AP and duplicate counts for a building are avoided.

Unfortunately the University of Colorado Boulder (UCB) system’s sampling
frequency is almost never constant and we have no control over when the system
opts to report device counts. User count sampling intervals are therefore
irregular, and vary across time and space. Sampling intervals generally range
from 3-12 minutes, and irregular outages in data collection can be hours or
days.

By taking the connected counts throughout the day we obtain time series data
which will be analyzed by NMF. The data sets analyzed here have been stripped
of all identifying information leaving only active user counts. As mentioned
earlier, this is one of the most attractive features of this dataset: complete
anonymity with no privacy rights breached. In the Results section of this
paper, we introduce and study the Norlin Library.

###  II-B Data Pre-processing and Matrix Structure

Interpolation
As with many time series techniques, it is helpful to have a standard sampling
frequency of the data for interpretation reasons. Non-constant sampling
frequencies would throw a wrench into the construction and interpretation of
the data matrix $X$. We settled upon ten-minute interpolation of the data
although, admittedly, this was a bit arbitrary. Due to the objective of
understanding general, aggregate behavior, we decided to perform a linear
interpolation as opposed to quadratic, cubic, etc.. We left all fractional
interpolated counts unchanged for the analysis since the NMF components will
not, in general, be whole numbers.

Matrix Embedding
Although factorization techniques act on matrices, most time series data is
typically one-dimensional. Therefore, a critical choice for this type of
analysis is how to impart 2-dimensional structure to the time series. Finding
ways to embed this data into a 2-dimensional matrix is crucial in terms of the
final product. Namely, how $W$ and $H$ are interpreted and how well the
product $WH$ approximates $X$ vary with the design of $X$.

In other time series approaches such as Singular Spectrum Analysis [30], data
is formed into a Hankel matrix, where each column consists of lagged vectors.
In our experience, NMF via a Hankel matrix interpretation has applications in
forecasting for streaming data (though there is an argument to be made about
information loss during Hankelization procedures). Instead, in a further
attempt to extract ”real” meaning, e.g. features from our factorization, we
seek to mirror the daily cycle of building use by choosing to represent each
24 hour period of data collection as a column of the design matrix, $X$. This
embedding is not without precedent: for large parts of the procedure we follow
the approach in this paper by Lee and Staneva [32] by treating each day of
data as a sample. This is justified by the clear daily oscillation in Wi-Fi
patterns coming from general group movements around campus. This leads to a
knowledge based factorization.

Using the ten-minute interpolation, we obtain 144 data points per day. Each
day is treated as a sample, and they appear chronologically as a columns in
$X$. Thus, the factors $W$ and $H$ have simple interpretations. Namely, the
columns of $W$ represent daily patterns and the rows of $H$ represent the
weights of the daily patterns for each day in the original data. This is where
the importance of constant sampling frequencies comes into play. Having the
columns of our data matrix referring the same times of the day is crucial for
a consistent comparison. Further, the columns would likely be of different
dimension without interpolation.

This embedding helps constrain the factorization and optimizes NMF’s ability
to produce results with physical interpretability. Given any day in the
training data, we can reconstruct its pattern via a linear combination of the
columns of $W$ where the coefficients (weights) are represented by a column of
$H$ corresponding to the desired day. It is also worth noting that we can
examine the prominence of each column of $W$ on any given day by examining the
weights in the $H$ matrix.

###  II-C Choice of Inner Dimension

Determining optimal choices for the dimension of the truncated SVD, $r$, has
been an outstanding issue for full-rank matrices. The same can be said for the
choice of inner dimension, $k$, in NMF since both parameters seek similar
results in terms of dimension reduction. For certain matrices, optimal hard-
thresholding techniques have been developed for SVD, but more-informal
strategies such as examining bends/elbows in a plot of decreasing singular
values are still used today [28]. We use a similar strategy for the
determination of $k$ via the graphical examination of mean-squared-error (MSE)
vs. $k$.

We look for the presence of sharp bends/elbows in this plot of MSE vs. $k$
using a sufficiently wide range of $k$ values. The entire NMF procedure is re-
computed for each successive value of $k$ and MSE is computed in the obvious,
element-wise manner between the resulting product $WH$ and $X$.

When in doubt, we advocate for potentially underestimating $k$ rather than
risking an overestimate. No matter what $k$ is specified to be, the NMF
algorithms will make factors $W$ and $H$ that fit the specified dimensions. If
$k$ is overestimated, we risk extracting components that are numerical
artifacts of the algorithm rather than true components of the data. If $k$ is
underestimated, then it is more likely that the components extracted are truly
parts of the data, or contain several hidden parts of the data. We would
rather see a few, true components of the data rather than many components that
might be a mix of truth and artifact.

###  II-D NMF Details

The problem setup for NMF is as follows. Given a data matrix,
$X\in\mathbb{R}^{n\times m}$, we seek factors $W\in\mathbb{R}^{n\times k}$ and
$H\in\mathbb{R}^{k\times m}$ (where the so called inner dimension $k\ll
min\\{n,m\\}$) having the properties

$X\approx WH\;\;\;\;\;\;s.t.\;W,H\succeq 0$ (1)

where $W,H\succeq 0$ denotes element-wise non-negativity. Similar to its
relative, SVD, dimension reduction is usually an assumed goal. The alternative
$X\approx WH^{T}$ is commonly used to describe the same factorization, but we
follow the notation of [24].

This factorization, although easy to describe, is quite difficult to solve in
general. Not only is it fundamentally ill-posed and non-convex, but it is NP-
hard in all but a very few specific cases [24]. With real data, the NP-
hardness is essentially guaranteed as exact low-rank structure rarely exists.

We note that many variants exist in terms of how $W$ and/or $H$ are
constructed. Twenty variants are listed and discussed in-depth in [24]. For
example, deep-NMF seeks $X\approx WH_{1}H_{2}H_{3}...H_{t}$ where
$W,H_{i}\succeq 0$ $\forall i$ and projective-NMF adds the additional
constraint $W=XH^{T}$, $H\succeq 0$. The style of NMF chosen is generally
application-dependent. In this paper, we stick to the traditional form of NMF
given by (1).

In general, no perfect factors are guaranteed to exist with NMF (in contrast
to SVD), which is a consequence of the non-negativity constraint (hence
the“$\approx$”). Therefore, NMF is set up as an optimization problem. For NMF,
there are at least four critical choices that are associated with any
optimization problem: Objective function, regularization terms, descent
algorithm and algorithm initialization. To give a more comprehensive look at
NMF, we will cover the objective function and regularization terms first and
follow with the algorithm and initialization procedure.

####  II-D1 Objective Function

Traditionally, NMF seeks

$\operatorname*{argmin}_{W\succeq 0,H\succeq 0}\left||X-WH\right||_{F}^{2}$
(2)

where $||\cdot||_{F}$ denotes the Frobenius norm. There is the choice of many
different objective functions that differ from the Frobenius norm, however.
The general objective function can be seen as the $\beta$-divergence, with
some choice of $\beta$ [24]. Although, we refer to reader to $l_{1}$-norm,
$l_{\infty}$-norm and weighted versions of (2) seen in [24]. We will not cover
these alternatives. The $\beta$-divergence between two non-negative scalars
$z$, $y$ is given by

$d_{\beta}(z,y)=\begin{cases}\frac{z}{y}-log(\frac{z}{y})-1&\beta=0\\\
zlog(\frac{z}{y})-z+y&\beta=1\\\
\frac{1}{\beta(\beta-1)}\left(z^{\beta}+(\beta-1)y^{\beta}-\beta
zy^{\beta-1}\right)&\beta\neq 0,1\end{cases}$ (3)

The three most common objective functions correspond to $\beta=0,1,2$. In
other words, the Itakura-Saito divergence (IS-divergence), the Kullback-
Leibler divergence (KL-divergence) and the Frobenius norm, respectively. For
our purposes, we chose the Frobenius norm ($\beta=2$) due to it successes in
Blind Source Separation (BSS) [29].

####  II-D2 Regularizations

It is also possible to regularize the problem given in (2) no matter the
objective function. Regularization induces certain solution types (e.g.
sparsity) and is often used to make the solution to the optimization problem
unique. For our purposes, we use the regularization option that is included in
the NMF function in Python and is part of the Scikit-Learn package [31]. Thus,
with $\beta=2$, the optimization problem becomes

$\begin{split}\operatorname*{argmin}_{W\succeq 0,H\succeq
0}\frac{1}{2}\left||X-WH\right||_{F}^{2}+\\\
\rho\alpha\left(\left||\text{vec}(W)\right||_{1}+\left||\text{vec}(H)\right||_{1}\right)+\\\
\frac{1}{2}\alpha(1-\rho)\left(\left||W\right||_{F}^{2}+\left||H\right||_{F}^{2}\right)\end{split}$
(4)

where $\alpha$ controls the intensity of the regularization and $\rho$
controls the relative importance of element-wise $l_{1}$-norms and the
Frobenius norms. As a potentially useful comparison, we note that (4) has many
similar properties to the elastic-net regression problem.

####  II-D3 Algorithm

There are two approaches to solve the problem given in (2) and (4). The first
is the classical gradient descent method. With linear problems gradient
descent is usually straightforward to implement since partial derivatives
exist and can easily be found. For example ([27]) showed that the gradient
descent algorithm is equivalent to the more widely-used multiplicative update
strategy, where each matrix is updated geometrically as opposed to the
arithmetic update of gradient descent.

A second widely-used strategy stems from recognizing that the problem given in
(2) is actually convex once either $W$ or $H$ are known. Given one above,
classical least squares can be used to solve for the other unknown.
Initializing $W$ or $H$ then implies a procedure where classical least squares
optimizes the other matrix during each iteration. This approach is called
alternating least squares and it is the basis for the coordinate descent
strategy.

####  II-D4 Initialization Strategy

Regardless of algorithm choice (multiplicative update or coordinate descent),
we must specify an initial $W$ and/or $H$ to start the procedure. Built-in NMF
solvers in most scientific computing languages (e.g. Scikit-Learn’s NMF
function for Python, MATLAB’s standard nnmf() function) initialize with random
matrices at default. So, each entry in $W$ and/or $H$ is drawn from a
$Uniform(0,1)$ distribution. Due to the non-convex nature of the problem and
the large number of unknowns (every entry in $W$ and $H$), random
initialization often results in slightly poorer outcomes and reproducibility
issues. We advocate for limiting the stochasticity of the algorithm, in favor
of more reproducible results via a SVD initialization strategy.

Namely, we use the $nndsvdar$ call within the Scikit-Learn NMF function. This
initialization comes from the non-negative double SVD (NNDSVD) strategy
proposed in Boutsidis et al.[26]. We remind the reader that, for any matrix
$A\in\mathbb{R}^{m\times n}$, there exists a decomposition

$A=U\Sigma V^{T}$ (5)

such that $U\in\mathbb{R}^{m\times m}$, $\Sigma\in\mathbb{R}^{m\times n}$ and
$V\in\mathbb{R}^{n\times n}$. The columns of $U$ are the eigenvectors of
$AA^{T}$, the columns of V are the eigenvectors of $A^{T}A$, and $\Sigma$ is
rectangular diagonal matrix filled with the singular values. For many
dimension reduction problems, truncated SVD is used to form a rank-$r$
approximation to the original matrix, $A$, via the Eckart-Young Theorem [28].

Standard non-negative double SVD (NNDSVD) examines a rank-$r$ approximation of
the original data matrix. The choice of the inner dimension, $k$, for NMF is
the same as $r$ for the truncated SVD. It then examines the positive and
negative components of each column of $U$ and row of $V^{T}$ and determines
which has a larger contribution (by comparing their norms) to each of the
singular vectors. The branch (positive or negative) with the larger
contribution is then chosen to replace the column of $U$ or row of $V^{T}$ in
its absolute value. The singular values are square-rooted and equally
distributed between the the $U$ and $V$ matrices. These scaled, non-negative
matrices $U$ and $V^{T}$ are then treated as the initialization to $W$ and
$H$, respectively. See Boutsidis and Gallopoulos [26] for the full algorithm
and discussion. The difference between this strategy and the $nndsvdar$ call
is that any values in the NNDSVD initialization equal to zero are replaced
with small, random values. The reason behind this is to allow zero entries to
be updated in the descent algorithm. In their basic forms, both coordinate
descent and multiplicative updates cannot alter zero entries.

##  III Results

Figure 1: Norlin Library Timeseries: raw data from the Spring of 2020\- A
large, unabridged data set was studied with the NMF technique described above.
Prominent features of this timeseries include strong daily cycles, weekly
variation, a large subset of corrupted data, and the beginning of the “COVID
effect” in Colorado.

Figure 2: $X$: All Days of the Study \- The same timeseries data as above
shown in a daily view. Alternatively, this could be considered columns of the
data matrix, X. Note the same types of variation as above: Weekends, Weekdays,
Error in Data and COVID.

###  III-A COVID Effect: Norlin Library

To illustrate this method’s ability to identify structure in a highly varied
dataset, it was applied to UCB’s main library: Norlin Library. This study was
carried out over an 11-week period from the Spring of 2020 spanning the
(reported) beginning of the COVID-19 epidemic in Colorado[33]. With regard to
the setup discussed in the previous section, our data matrix
$X\in\mathbb{R}^{n\times m}$, where $n=144$ (the number of sampling intervals
in a day) and $m=77$ (77 unique days in this study).

Figure 3: Norlin Library: Columns of W or “Typical Days”. After NMF with $k=4$
(see justification in Results), these are the extracted basis vectors by which
the algorithm represents every day in the data set. Also in this figure are
the “Average geometry of Weekdays” and “Average Geometry of Weekends”. Keep in
mind all days are represent by some linear combination of these 4 signals.
We’ve named the columns, $W_{i}$, based on when they represent users in the
library. NOTE: There is no scale on the y-axis scale in this figure, as the
magnitudes of these curves is somewhat arbitrary.

Figure 4: $H$: Activation Contribution to each day’s profile \- Here are the
activations for each day multiplied by the L1 Norm of their corresponding
Column of W. In this way, we can compare the relative signal strengths of
columns of $W$. Do these results make physical sense?

Norlin Library is a major building located north and central on CU’s Boulder
main campus. It houses the Norlin “Stacks”—-a major collection of humanities,
social science and life science documents, and rare volumes-— and includes
diffuse study areas and a coffee shop.

Figure 5: What inner dimension to choose? Relative error vs. Inner dimension
for the NMF factorization of the Norlin Library Study Data strongly supports
the choice of $k=4$ in this study.

Because very few if any classes are taught at Norlin, the building behaves
more as a typical commercial building: there are fewer coordinated movements
of large crowds into and out of the building. Since there are no residents, a
stronger daily structure is observed in the user count data than at the
residential halls. At night, the building is empty aside from the baseline
count of latent Wi-Fi users (networked hardware such as vending machines, web
cameras, and permanent computers). In this way, the flow of humans (and their
devices) might better approximate the surge typical of commercial buildings in
Smart & Connected Cities.

This study began on 01/19/2020 and ran for 11 weeks through 04/04/2020. In
this study, the data was not segregated into “like” days nor were outliers or
corrupted data removed. Rather, the entire unabridged interval of this time
series was passed wholesale to the algorithm.

Before discussing the results of the factorization, let’s take a closer a look
at the original time series data. This interval demonstrates some intriguing
features—some regular, some irregular—that make it an insightful case study of
the techniques described in this paper. Some major ones to note, which are
also highlighted in Figure 2 and Figure 2:

  * •

Weekdays vs. Weekends: User count populations demonstrate much different
structure from weekdays to weekends.

  * •

Corrupted Data: for a roughly 3 day period from 3/01 to 3/03, the data is
obviously corrupted and does not reflect the true behavior of devices in the
library. Typically, many techniques would require operator intervention to
throw out this data. Here it is included to demonstrate the robustness of this
technique.

  * •

COVID-19: Perhaps most conspicuous is the obvious “paradigm” shift from the
pre-COVID era to the COVID era in Colorado. On March 5th, the Governor of the
state announced the first confirmed case of COVID in Colorado. Shortly
thereafter the user count data reflects the steady de-peopling of campus.

Per the procedure prescribed in the previous section, this time series data
was vectorized by day, loaded into a data matrix and factorized. As before, a
variant of Hierarchical Alternating Least Squares (HALS) algorithm was used to
solve the factorization, and initialized by a non-zero SVD method. For this
application, a relatively low inner dimension of $k=4$ was chosen and this
rank-selection is supported by the results of an MSE analysis (Figure 5).

The resulting factorization successfully found a representation of the data by
parts that seems to reflect plausible structures in the movement of device
populations. The 4 columns of $W$ representing the constituent signals of this
system are shown in Figure 4. They have been assigned descriptive names based
on when they represent devices in the Library: $W_{1}$ (AM), $W_{2}$ (PM),
$W_{3}$ (Late) and $W_{4}$ (Baseline). Keep in mind that the NMF has “learned”
a representation of every day in our data set as a linear combination of these
4 vectors. Examining the weight matrix, $H$ gives some insight into the
temporal structure of user counts in the Library. In Figure 4, we show each of
the 4 weights for every day in the times series. These weights are scaled by
the Manhattan distance of their corresponding columns of W, or $W_{j}$:

$h_{i,j,weighted}=h_{i,j}\times\left\|W_{j}\right\|_{1}\times\text{10 min}$
(6)

This scaling is performed to normalize the weights for comparison and impart a
physical intuition to them. $\left\|W_{j}\right\|_{L_{1}}$ is an effective
measure of cumulative device time each day or when multiplied by 10 minutes
can be considered more directly a “device minutes per day” value (due to the
imposed 10 minute sampling interval) represented by each column of $W$
($W_{j}$). Numerical integration may provide a more refined value, but either
is sufficient for understanding the relative contributions of each column of
$W$ for each date. This new value $h_{i,j,weighted}$ therefore allows us to
directly compare the “signal strengths” of each $W_{j}$ relative to one
another, and more critically, how those relationships change temporally.
Several patterns can be identified:

  * •

$W_{1}(AM)$ is strongly active on weekdays and essentially “off” on weekends.
This signal could correspond to the “traditional crowd” that only comes to
campus during weekdays.

  * •

$W_{3}(Late)$ shows a similar pattern, though curiously is offset a day (“off”
on Friday and Saturday). This signal could correspond to the “homework crowd”,
i.e. those individuals who take a break from class and homework at the end of
the week but return Sunday night to do work before the coming week.

  * •

The entire period of error in data collection has a drastically different
profile and stands out clearly from anywhere else on the plot. $W_{4}$ sees a
40-fold increase in expression compared to anywhere else in the sampling
period, while $W_{2}$ sees a unique reduction. It is an interesting finding
that NMF can detect these types of errors unsupervised.

  * •

As COVID begins, the signals for $W_{1}$ (AM), $W_{2}$ (PM), $W_{3}$ (Late)
decay as the campus empties, but the magnitude of $W_{4}$ actually stays
constant around 450 device hours per day. This corresponds to approximately 20
devices on the network day, which indeed appears to be the background device
count.

##  IV Conclusion

NMF was chosen as a decomposition for this data due to its high
interpretability. As discussed earlier though, there are no guarantees made on
the uniqueness or indeed suitability of the results of NMF factorization.
Rather, results must be analyzed and interpreted to confirm physical meaning.
Separating the data from the time series into daily vectors is a natural way
to frame our factorization and imparts a physical intuition to the problem.

In the Norlin Library study, a complex, real-world data set was passed through
the algorithm described in this paper. Despite being highly varied (including
corrupted data), the resulting factorization found structure in the population
time series that reflected the real world use of the library.

Important next steps would be to study how increasing the choice of inner
dimension of the factorization may identify finer features in the population
data or using a separate data source (such as transportation data) to identify
specific instances of coordinated arrivals or departures into the Wi-Fi
environment.

We are excited by the techniques described and believe they may be deployed to
similar use cases in order to quickly, automatically and anonymously find
patterns in population movements within Smart & Connected environments.

## Acknowledgments

Thanks to OIT and the Campus Chief Security Officer for making Wi-Fi data
available. Many of the ideas and framing were developed through interactions
with colleagues in the CU Boulder Technology, Cybersecurity and Policy (TCP)
program, e.g Dan Massey, Kevin Gifford, Keith Gremban, and Lloyd Thrall.
Seminar presentations and speakers had a formative impact on this work. Thanks
to Peter Borrell and Jake McGrath for many valuable formative discussions.
Finally, special thanks to Sandia National Laboratory for the support of this
research project.

## References

  * [1] United Nations, Department of Economic and Social Affairs, Population Division (2019). World Urbanization Prospects: The 2018 Revision (ST/ESA/SER.A/420). New York: United Nations.
  * [2] Edward L. Glaeser ’I&’ Bruce Sacerdote, 1996.“Why Is There More Crime in Cities?,” NBER Working Papers 5430, National Bureau of Economic Research, Inc.
  * [3] Louf, R., Barthelemy, M. “How congestion shapes cities: from mobility patterns to scaling”. Sci Rep 4, 5561 (2014). https://doi.org/10.1038/srep05561
  * [4] Schwela,, Dietrich. “Air Pollution and Health in Urban Areas” Reviews on Environmental Health, vol. 15, no. 1-2, 2000, pp. 13-42. https://doi.org/10.1515/REVEH.2000.15.1-2.13
  * [5] Kristian Behrens, Frédéric Robert‐Nicoud, “Survival of the Fittest in Cities: Urbanisation and Inequality, The Economic Journal”, Volume 124, Issue 581, 1 December 2014, Pages 1371–1400, https://doi.org/10.1111/ecoj.12099
  * [6] Smart and Connected Communities (S & amp;CC). Beta site for NSF - National Science Foundation. (n.d.). Retrieved November 18, 2021, from https://beta.nsf.gov/funding/opportunities/smart-and-connected-communities-scc.
  * [7] F. Imbault, M. Swiatek, R. de Beaufort and R. Plana, “The green blockchain: Managing decentralized energy production and consumption,” 2017 IEEE International Conference on Environment and Electrical Engineering and 2017 IEEE Industrial and Commercial Power Systems Europe (EEEIC / I & CPS Europe), 2017, pp. 1-5, doi: 10.1109/EEEIC.2017.7977613.
  * [8] M. Lom, O. Pribyl and M. Svitek, “Industry 4.0 as a part of smart cities,” 2016 Smart Cities Symposium Prague (SCSP), 2016, pp. 1-6, doi: 10.1109/SCSP.2016.7501015.
  * [9] J. Satyakrishna and R. K. Sagar, “Analysis of smart city transportation using IoT,” 2018 2nd International Conference on Inventive Systems and Control (ICISC), 2018, pp. 268-273, doi: 10.1109/ICISC.2018.8399077.
  * [10] S. A. Shah, D. Z. Seker, M. M. Rathore, S. Hameed, S. Ben Yahia and D. Draheim, “Towards Disaster Resilient Smart Cities: Can Internet of Things and Big Data Analytics Be the Game Changers?,” in IEEE Access, vol. 7, pp. 91885-91903, 2019, doi: 10.1109/ACCESS.2019.2928233.
  * [11] G. Khodabandelou, V. Gauthier, M. El-Yacoubi and M. Fiore, “Population estimation from mobile network traffic metadata,” 2016 IEEE 17th International Symposium on A World of Wireless, Mobile and Multimedia Networks (WoWMoM), 2016, pp. 1-9, doi: 10.1109/WoWMoM.2016.7523554.
  * [12] Douglass, R.W., Meyer, D.A., Ram, M. et al. “High resolution population estimates from telecommunications data”. EPJ Data Sci. 4, 4 (2015). https://doi.org/10.1140/epjds/s13688-015-0040-6
  * [13] K. Zhang, J. Ni, K. Yang, X. Liang, J. Ren and X. S. Shen, “Security and Privacy in Smart City Applications: Challenges and Solutions,” in IEEE Communications Magazine, vol. 55, no. 1, pp. 122-129, January 2017, doi: 10.1109/MCOM.2017.1600267CM.
  * [14] IEEE 802.11-2016 - IEEE standard for information technology–telecommunications and information exchange between systems local and metropolitan area networks–specific requirements - part 11: Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) specifications. IEEE SA - The IEEE Standards Association - Home. (n.d.). Retrieved November 18, 2021, from https://standards.ieee.org/standard/802_11-2016.html.
  * [15] CISCO (2020) Cisco Annual Internet Report (2018-2023). White Paper. https://www.cisco.com/c/en/us/solutions/collateral/executive-perspectives/annual-internet-report/white-paper-c11-741490.pdf
  * [16] Wang, Zhe, Hong, Tianzhen, Piette, Mary Ann, and Pritoni, Marco. “Inferring occupant counts from Wi-Fi data in buildings through machine learning”. United States: N. p., 2019. Web. https://doi.org/10.1016/j.buildenv.2019.05.015.
  * [17] Ouf, Mohamed M., Issa, Mohamed H., Azzouz, Afaf and Sadick, Abdul-Manan 2017, “Effectiveness of using Wi-Fi technologies to detect and predict building occupancy”, Sustainable buildings, vol. 2, pp. 1-10, doi: 10.1051/sbuild/2017005.
  * [18] U. Singh, J. -F. Determe, F. Horlin and P. D. Doncker, “Crowd Forecasting Based on Wi-Fi Sensors and LSTM Neural Networks,” in IEEE Transactions on Instrumentation and Measurement, vol. 69, no. 9, pp. 6121-6131, Sept. 2020, doi: 10.1109/TIM.2020.2969588.
  * [19] Bellini, P. et al. “Wi-Fi based city users’ behaviour analysis for smart city.” J. Vis. Lang. Comput. 42 (2017): 31-45.
  * [20] Kontokosta, C. E., & Johnson, N. (2017). Urban phenology: Toward a real-time census of the city using Wi-Fi data. Computers, Environment and Urban Systems, 64, 144-153. https://doi.org/10.1016/j.compenvurbsys.2017.01.011
  * [21] Buciu I, Pitas I (2004) Application of non-negative and local non negative matrix factorization to facial expression recognition.
  * [22] Shahnaz F, Berry M (2006) Document clustering using nonnegative matrix factorization. Information Processing and Management: An International Journal 42: 373–386
  * [23] Lee, D. D., and Seung, H. S. (1999) “Learning the parts of objects by non-negative matrix factorization” Nature 401(6755), 788-791
  * [24] Gillis, Nicolas. Nonnegative Matrix Factorization. Society for Industrial and Applied Mathematics, 2021.
  * [25] Paatero, P., Tapper, U., 1994. “Positive matrix factorization: a non-negative factor model with optimal utilization of error estimates of data values.” Environmetrics 5, 111-126.
  * [26] Boutsidis C, Gallopoulos E. (2008) “SVD based initialization: A head start for nonnegative matrix factorization” Science Direct, 1350-1362
  * [27] Cichocki A, Anh H.P. (2009) ”Fast Local Algorithms for Large Scale Nonnegative Matrix and Tensor Factorizations” IEICE Transactions on Fundamentals of Electronic Communications and Computer Science
  * [28] Brunton, Steven Lee., and Jose Nathan. Kutz. Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control. Cambridge University Press, 2019.
  * [29] Sawada, Hiroshi & Ono, Nobutaka & Kameoka, Hirokazu & Kitamura, Daichi & Saruwatari, Hiroshi. (2019). A review of blind source separation methods: Two converging routes to ILRMA originating from ICA and NMF. APSIPA Transactions on Signal and Information Processing. 8. 10.1017/ATSIP.2019.5.
  * [30] Golyandina, N., and Zhigljavsky, A. A.,Singular Spectrum Analysis for Time Series, Springer Briefs in Statistics, Springer-Verlag Berlin Heidelberg, 2013.
  * [31] Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011.
  * [32] Lee, Wu-Jung, and Valentina Staneva. “Compact Representation of Temporal Processes in Echosounder Time Series via Matrix Decomposition.” The Journal of the Acoustical Society of America, vol. 148, no. 6, 2020, pp. 3429–3442., doi:10.1121/10.0002670.
  * [33] IEEE 802.11-2016 - IEEE standard for information technology–telecommunications and information exchange between systems local and metropolitan area networks–specific requirements - part 11: Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) specifications. IEEE SA - The IEEE Standards Association - Home. (n.d.). Retrieved September 16, 2021, from https://standards.ieee.org/standard/802_11-2016.html.