# Deep Graph Learning for Spatially-Varying Indoor Lighting Prediction

Jiayang Bai, Jie Guo$\dagger$, Chenchen Wan, Zhenyu Chen, Zhen He, Shan Yang,
Piaopiao Yu, Yan Zhang and Yanwen Guo$\dagger$ State Key Lab for Novel
Software Technology, Nanjing UniversityNanjingChina

(2022)

###### Abstract.

Lighting prediction from a single image is becoming increasingly important in
many vision and augmented reality (AR) applications in which shading and
shadow consistency between virtual and real objects should be guaranteed.
However, this is a notoriously ill-posed problem, especially for indoor
scenarios, because of the complexity of indoor luminaires and the limited
information involved in 2D images. In this paper, we propose a graph learning-
based framework for indoor lighting estimation. At its core is a new lighting
model (dubbed DSGLight) based on depth-augmented Spherical Gaussians (SG) and
a Graph Convolutional Network (GCN) that infers the new lighting
representation from a single LDR image of limited field-of-view. Our lighting
model builds 128 evenly distributed SGs over the indoor panorama, where each
SG encoding the lighting and the depth around that node. The proposed GCN then
learns the mapping from the input image to DSGLight. Compared with existing
lighting models, our DSGLight encodes both direct lighting and indirect
environmental lighting more faithfully and compactly. It also makes network
training and inference more stable. The estimated depth distribution enables
temporally stable shading and shadows under spatially-varying lighting.
Through thorough experiments, we show that our method obviously outperforms
existing methods both qualitatively and quantitatively.

Lighting, Graph learning, Rendering

$\dagger$ is the corresponding author.

††copyright: acmlicensed††journalyear: 2022††journal: TOG††ccs: Computing
methodologies Mixed / augmented reality††ccs: Computing methodologies Scene
understanding Figure 1. Our method takes a single image as the input and
predicts the full and spatially-varying indoor illumination that can be used
to generate consistent shading and realistic shadows after inserting virtual
objects into the image.

##  1\. Introduction

Predicting indoor illumination from a single image is a challenging but
meaningful task (Einabadi et al., 2021). It is a fundamental step for many
applications, ranging from indoor scene understanding to augmented reality
(AR). However, this problem is highly ill-posed due to the following reasons.
First, indoor illumination stems from a wide range of different luminaires
including spot lights, large area lights and indirect lighting. This is quite
different from outdoor lighting which is dominated by the sun (Liu et al.,
2015; Hold-Geoffroy et al., 2017; Zhang et al., 2019; Hold-Geoffroy et al.,
2019; Yu et al., 2021). Second, a single image captured by an ordinary camera
has limited field-of-view (FoV) that may fail to provide sufficient cues of
indoor light sources. In fact, many indoor light sources have a local impact
and yield spatially-varying lighting effects. This makes them hard to be
identified from a single FoV-limited image. Moreover, strong lighting-
reflectance ambiguities exist in indoor scenes since different combinations of
lighting and surface reflectance may produce the same pixel intensities.

Currently, this problem is addressed in two ways. One line of work assumes the
availability of some prior knowledge about the scene, e.g., scene geometries
(Barron and Malik, 2013; Maier et al., 2017; Weber et al., 2018) obtained by
depth sensors or multi-view inputs. Prior knowledge is effective for
overcoming the ill-posedness of the problem but will limit the applicability
of the methods. A recent trend is to leverage deep learning based solutions
accompanied with large training datasets to infer lighting from a single,
casually captured image (Gardner et al., 2019, 2017b; Garon et al., 2019; Song
and Funkhouser, 2019; Srinivasan et al., 2020; Zhan et al., 2021b, a). The
lighting models used in these solutions can be divided into two categories. In
the first category, the indoor lighting of surrounding environment is stored
in a densely sampled environment map (Gardner et al., 2017b), i.e., image-
based lighting (IBL). Unfortunately, inferring a full environment map directly
is very difficult due to the high dimensionality of the representation.
Alternatively, the second category leverages analytical models with a small
number of parameters to fit real-world illumination. For instance, indoor
illumination can be approximated by a set of Spherical Harmonics (SH) (Cheng
et al., 2018; Garon et al., 2019; Green, 2003) or Spherical Gaussians (SG)
(Gardner et al., 2019; Tsai and Shih, 2006a; Wang et al., 2009; Tsai and Shih,
2006b; Zhan et al., 2021b). Generally, SH works well for low-frequency
lighting, but is plagued with ringing artifacts as the SH order increases.
Gardner et al. (2019) proposed using 2-5 SGs to model high-frequency light
sources in an indoor scene. Although it outperforms SH-based methods,
regressing the positions of floating SGs is unstable for deep learning.
Consequently, this method achieves degraded performance as the number of SGs
increases.

To retain the advantages of SGs in compactly encoding high-frequency lighting
and enable stable training, we introduce _DSGLight_ , a novel SG-based
lighting model augmented with depth values. In this model, $N$ SGs are assumed
to have fixed positions and are evenly distributed over a unit sphere. We also
fix the bandwidth of each SG such that only its color and depth values vary.
In this case, predicting indoor lighting boils down to estimating $N$ RGBD
values, which is more stable and robust than regressing several free SGs on
the sphere. The augmented depth values allow us to generate spatially-varying
shading and shadows with respect to the insert location.

Considering the non-Euclidean nature of DSGLight, we design a graph
convolution network (GCN) for feature extraction and lighting prediction. GCNs
have gained impressive successes in many vision tasks, but are less exploited
in lighting prediction. In our network, we first extract features from an
input image with vanilla convolutional layers. Then, these features are
connected to GCN features with several fully-connected layers. After
information exchanging across neighboring nodes, we finally obtain the color
and depth for each SG. We demonstrate through extensive experiments that our
method outperforms existing methods both qualitatively and quantitatively. The
state-of-the-art performance on indoor lighting prediction allows us to
achieve more realistic shading results when inserting virtual objects into
real indoor scenes (see Fig. 1).

In summary, our main contributions in this paper are:

  * •

A depth-augmented indoor lighting model, i.e., DSGLight, to compactly and
faithfully encode spatially-varying lighting in an indoor scene.

  * •

A dedicated graph convolutional network for predicting DSGLight from a single
image.

  * •

A new dataset with supervision of DSGLight parameters that is constructed to
train our network.

##  2\. Related Work

###  2.1. Indoor lighting prediction

Inferring the illumination from 2D images is a classic problem that has been
extensively studied in the past decades. Early progress on illumination
estimation relied on prior knowledge of scene geometry. For instance, Barron
and Malik (2013) adopted depth sensors to capture geometric information, Wu et
al. (2011) reconstructed geometry with multi-view stereo, and Larsch et al.
(2011) required users to annotate the geometry. Some other works manage to
predict illumination from purely 2D image(s). Khan et al. (2006) presented a
method to flip the given image to approximate the rest of the whole panorama,
which are out of view. Unfortunately, this only works well in some special
cases.

While deep learning is significantly successful in many graphical and vision
tasks, some learning-based methods are proposed to estimate illumination.
Hold-Geoffroy et al. (2017) presented an end-to-end approach that leveraged
CNN to predict outdoor illumination represented by a low-dimensional
analytical model, significantly outperforming previous traditional methods.
Gardner (2017b) also adopted an end-to-end CNN to regress a non-parametric
environment map to represent accurate illumination. Song et al. (2019)
proposed “Neural Illumination” that can predict indoor lighting at a specific
locale. Another recent work (Srinivasan et al., 2020) leveraged a multi-scale
volumetric lighting representation to estimate spatially-coherent illumination
from images.

Close to our work, several methods also leverage multiple SGs to represent
indoor lighting (Gardner et al., 2019; Zhan et al., 2021b). However, these
methods only capture dominated lights and do not account for spatially-varying
environmental lights, resulting in fewer details for rendering specular
objects and inconsistent shadows. Moreover, these methods usually leverage
CNNs to establish the mapping between input images and SGs. However, CNNs,
with a limited receptive field, fail to explore the long-range correlations
between SGs.

###  2.2. Graph convolution networks

To deal with non-Euclidean data like graphs, Graph convolution networks (GCNs)
provide well-suited solutions which can be divided into two types: spatial-
based (Atwood and Towsley, 2015; Micheli, 2009; Niepert et al., 2016) and
spectral-based (Defferrard et al., 2016; Hamilton et al., 2017; Kipf and
Welling, 2017). Spatial-based methods used in many works are mainly based on
the spatial relations of vertices in graphs and they apply spectral
convolutions on graph using ideas from graph signal processing (Shuman et al.,
2012). As transforming the signal to spectral domains can be expensive, many
recent works (Defferrard et al., 2016; Kipf and Welling, 2016) approximate the
spectral convolutions using Chebyshev polynomials. Similarly, our network also
uses Chebyshev polynomials for spectral convolutions in GCN layers. We use
GCNs to extract features from indoor illumination encoded in a non-Euclidean
form and establish the connection between 2D images and a SG-based lighting
model.

##  3\. DSGLight

(a)

(b)

(c)

Figure 2. Illustration of DSGLight. Left: centering directions of 128 evenly
distributed SGs. Middle: fitted DSGLight (RGB channels) representation for an
indoor panorama. Right: fitted DSGLight of depth representation.

(a) Original

(b) DSGLight (RGB channels)

(c)  (Gardner et al., 2019) (3 SGs)

(d)  (Gardner et al., 2019) (5 SGs)

(e)  (Garon et al., 2019) (SH)

Figure 3. Visual comparisons of different lighting representations for an
indoor scene.

###  3.1. Representation

DSGLight is a depth-augmented parametric lighting model that leverages a set
of discrete SG functions to represent the spatially-varying indoor
illumination, as illustrated in Fig. 2. We denote DSGLight by
$\mathcal{M}=\\{m_{i}|i=1,\dots,N\\}$ which contains $N$ vertices. Each vertex
in DSGLight is associated with a spherical Gaussian (Wang et al., 2009)
defined as follows:

(1)
$G(\mathbf{v};\mathbf{a},\lambda,\mu)=\mathbf{a}e^{\lambda(\mathbf{v}\mu-1)}$

where $\mu$ is the lobe axis (the centering direction),
$\lambda\in(0,+\infty)$ is the lobe sharpness (the bandwidth), and
$\mathbf{a}$ is the lobe amplitude ($\mathbf{a}\in\mathbb{R}^{4}$ for RGBD
information). The direction $\mathbf{v}$ denotes the spherical parameter of
the resulting function. Given $m_{i}=\\{\mathbf{a}_{i},\lambda_{i},\mu_{i}\\}$
for $i^{th}$ vertex, the overall indoor illumination encoded in a given
panorama is compressed into an arbitrary spherical function $f(\mathbf{v})$ as
follows:

(2) $f(\mathbf{v})=\sum_{i}^{N}G(\mathbf{v};m_{i}).$

The corresponding depth panorama is treated in a similar way. Currently, we
set $N$ to 128.

To further reduce the number of parameters, we fix the position $\mu$ and the
bandwidth $\lambda$ of each SG. Specifically, 128 nodes are placed evenly over
a unit sphere, and they share the same bandwidth. This allows $\mu$ and $v$ to
be preset for our model. The distribution of these nodes is shown in Fig. 2.
Besides being a more compact representation, such a design stabilizes network
training and inferences. As pointed out in (Gardner et al., 2019), regressing
several floating SGs with floating positions and unconstrained bandwidth is
unstable, especially for many SGs. As a result, performance decreases from 3
to 5 SGs. With fixed position and bandwidth, inferring DSGLight boils down to
predicting $N$ RGBD values, which is easier and more stable than inferring
floating SGs.

To make $N$ SGs cover the whole panorama as much as possible and avoid
overlapping, we set the bandwidth parameter $\lambda$ as follows:

(3) $\lambda=\frac{\ln{0.6}}{\cos(\arctan(2/\sqrt{N}))-1}.$

Details are provided in the supplemental material.

###  3.2. Comparison and analysis

Previous works have proposed many parametric methods to reduce the number of
parameters in representing panoramic lighting, such as SH, SG and IBL. Here,
we compare these representations with our DSGLight. One example comparing
these representations is provided in Fig. 3.

Among these methods, IBL is the most accurate. However, regressing the high-
dimensional IBL directly may lead to local minimal and it is difficult to
converge during training. With much fewer parameters, SH is a low-dimensional
lighting representation and can be easily predicted with a compact decoder.
Unfortunately, SH is too sensitive to its parameters’ variation. Moreover, the
network has a tendency to overfit to the ringing artifacts when employing
high-order SH functions to preserve high-frequency information.

SG is able to represent high-frequency lighting and avoid ringing artifacts.
However, Gardner et al. (2019) show that it is hard to optimize the positions
of SGs, especially when the number of SGs is large. In comparison, DSGLight
represents an indoor panorama with a certain number of SG functions which have
fixed positions and bandwidth. It is more stable and effective than predicting
arbitrary and floating light sources. In addition, as the vertices are spread
evenly on the unit sphere and cover the entire sphere, DSGLight can capture
both dominated light sources and detailed environmental lighting. This allows
high-quality AR rendering with consistent shading and shadows, as well as fine
details on specular objects.

###  3.3. Dataset preparation

To generate the dataset for training, we have to prepare a large number of
images and their corresponding DSGLight parameters. Given a HDR panorama,
generating ground truth DSGLight can be solved as an optimization problem,
i.e., minimizing the difference between illumination/depth constructed by
DSGLight and the ground truth panorama:

(4)
$\displaystyle\boldsymbol{A}^{*}=\underset{\boldsymbol{A}}{\arg\min}\sum_{p\in\mathcal{P}}\left(p-\sum_{i=1}^{N}G\left(F(p);\mathbf{a}_{i},\lambda_{i},\mu_{i}\right)\right)^{2}$

where
$\boldsymbol{A}=\\{\mathbf{a}_{1},\mathbf{a}_{2},\dots,\mathbf{a}_{N}\\}$ is
the set of lobe amplitudes of SGs, $p$ is a pixel in the ground truth
panoramas $\mathcal{P}$, and $F$ maps $p$ to a direction $\mathbf{v}$ in the
sphere. When fixing $\mathcal{P}$, $\mu$ and $\lambda$, the problem boils down
to a linear optimization problem, which can be solved easily by nonnegative
least square algorithms.

We currently perform the optimization on two datasets:

  * •

Laval dataset (Gardner et al., 2017a) contains 2,200 high resolution indoor
panoramas which are fully HDR without any saturation. This dataset also
contains depth values for each panorama.

  * •

SUN360 Dataset (Xiao et al., 2012) contains 12,000 indoor LDR panoramas
captured in a wide range of indoor scenes. For this dataset, we first convert
each LDR panorama into a HDR panorama with the trained network of (Eilertsen
et al., 2017). This dataset does not contain depth values. Therefore, we can
only handle DSGLight with RGB channels when using examples from this dataset.

We further employ Mitsuba renderer (Nimier-David et al., 2019) to sample
several Fov-limited images from each panorama at random elevation between
-20°and 20°, azimuth between -180°and 180°and FoV between 60°and 80°.
Meanwhile, the corresponding panorama is captured by aligning its center with
the view direction.

Figure 4. The overall pipeline of our method. Given an input image, we first
leverage a backbone CNN as a feature extraction module (a) to learn the
feature of nodes in DSGLight. Then, we construct a graph with feature vectors
initialized by the feature extraction module. A GCN (b) is employed to
exchange feature vectors across adjacent nodes and finally predicts the values
(color or depth) of the SG nodes.

##  4\. Deep Graph Learning for DSGLight

Our goal is to recover DSGLight of the environment in which the photograph has
been taken. It does not simply assign light intensity and depth value to each
visible light source, but must extrapolate large portions of the invisible
illumination. Thus, it requires wide context information and must exploit very
specific relationships between light sources to guide estimation. To this end,
we exploit a graph convolutional network (GCN) to use prior knowledge on the
non-Euclidean nature of graph structure of DSGLight and cope with significant
amounts of potential relationships between nodes. Running graph convolutions
updates the indoor-specific perceptual features extracted from the input
image, which is equivalent as driving light sources’ intensities towards the
correct light intensities in the RGB color space.

###  4.1. Network architecture

Our framework consists of two parts: a feature extraction module and a graph
convolutional network (GCN). Fig. 4 shows the architecture of our network.

In the feature extraction module, we can use any CNN based backbone to extract
the features of an input image. Currently, we adopt the VGG-19 (Simonyan and
Zisserman, 2015) as the backbone and generate a $512\times 7\times 11$ feature
matrix given an image with a resolution of $240\times 360$. This feature
matrix is further connected to GCN with two fully-connected layers, which
transform the feature map to the features of each node in DSGLight.

In our framework, a graph convolutional layer takes a graph as input and
outputs a graph with a different feature matrix. It can be written as the
following non-linear function:

(5)
$\mathbf{H}^{l+1}=\sigma\left(\mathbf{E}\mathbf{H}^{l}\mathbf{W}^{l}\right)$

where $\mathbf{H}^{l}$ is the input of the $l^{th}$ layer, $\mathbf{W}^{l}$ is
a trainable weight matrix, $\mathbf{E}$ is the adjacency matrix describing the
graph structure in the matrix form, and $\sigma$ is a non-linear activation
function (ReLU in our experiments). Stacking graph convolutional networks can
combine graph node features and graph topological structural information to
make final predictions.

GCN layer takes the node representations $\mathbf{H}^{l}$ from previous layer
as inputs and outputs new node representations. Thus we construct the graph
$\mathcal{G}=\\{\mathcal{V},\mathcal{F},\mathbf{E}\\}$ as the input of GCN.
GCNs are vertex-based and perform convolution on vertices with edges denoting
connections. $\mathcal{V}$ is the set of nodes corresponding to vertices in
DSGLight. $\mathcal{F}$ is the features of vertices learned from CNN.
$\mathbf{E}$ is the adjacency matrix in Eq. (5). To get the adjacency matrix,
we use k-Nearest Neighbors (k-NN) to construct the directed dynamic edges
between nodes.

Our GCN module comprises 4 graph convolutional layers. Li et al. (2019) stated
that stacking more layers into a GCN would easily lead to the vanishing
gradient problem and over-smoothing. Our experiments show that four GCN layers
are sufficient for multi-scale feature extraction and are stable for training.

###  4.2. Loss functions

In the design of loss functions, we handle color and depth separately. For
depth, we only use L2 loss. For color, we design a hybrid loss that consists
of a reconstruction loss $\mathcal{L}_{R}$, a perceptual loss
$\mathcal{L}_{VGG}$ and a weighted L2 loss $\mathcal{L}_{W}$ for the SG
parameters, _i.e._ ,

(6) $\mathcal{L}=\mathcal{L}_{W}+\alpha\mathcal{L}_{R}+\beta\mathcal{L}_{VGG}$

where $\alpha$ and $\beta$ are the weights to balance different losses. In our
training, we set $\alpha$ and $\beta$ to 0.2 and 0.1, respectively.

Reconstruction loss penalizes the per-pixel intensity distance between the
ground truth and the predicted environment map reconstructed with Eq. (2). The
perceptual loss aims to enhance the perceptual similarity in the predicted
illumination. We adopt a VGG-19 architecture which is pre-trained in the image
classification task. Let $\Phi$ be the output of the last pooling layer in
VGG-19, and $N_{\Phi}$ represents the total number of pixels in this feature
space. The perceptual loss is defined as:

(7)
$\mathcal{L}_{VGG}(p,\hat{p})=\frac{1}{N_{\Phi}}\ell_{2}(\Phi(p),\Phi(\hat{p}))$

where $p$ and $\hat{p}$ denote the ground truth and the predicted environment
map, separately.

L2 loss is used mostly in regression problems but it will lose high-frequency
information in training. Since the DSGLight has a wide range of values, we
reweight SGs in the DSGLight to boost the importance of salient SGs. Inspired
by ACESFilm (Academy Color Encoding System) curve, which distributes the
weights to pixels of HDR images empirically and achieves good tone mapping
result, we use the following weighted function:

(8) $f_{w}(x)=\frac{x(2.51x+0.03)}{x(2.43x+0.59)+0.14}.$

The weighted L2 loss for SG parameters is calculated as :

(9)
$\mathcal{L}_{W}(\boldsymbol{A},\boldsymbol{\hat{A}})=f_{w}(\boldsymbol{A})\ell_{2}(\boldsymbol{A},\boldsymbol{\hat{A}})$

where $\boldsymbol{A}$ and $\boldsymbol{\hat{A}}$ represent the ground truth
and our predicted parameters, respectively.

###  4.3. Discussion on the locality of indoor lighting

Most lights in indoor scenes have a local influence. This means we not only
need to predict the distribution of lighting around current point of view, but
also have to estimate the relative positions of the lights. To handle depth,
we augment DSGLight with depth values. In this case, the lobe amplitude
parameter of each SG in DSGLight contains four channels (three for color and
one for depth).

The training requires a large number of panoramas and their ground truth depth
maps. However, among our two datasets, only Laval dataset has manually
annotated per-pixel depth. From these annotated depth maps, we optimize the
SGs of depth according to Eq. (4) where $\mathcal{P}$ is replaced by the depth
map.

Some previous methods, such as (Li et al., 2020) and (Garon et al., 2019),
also try to handle spatial variations of lighting for indoor scenes. Unlike
ours, they estimate indoor lighting for each inserted point independently,
yielding a local representation of lighting that varies spatially. Although
they produce finer spacial light prediction, their strategy easily incurs
temporal flickering when the inserted virtual object moves in the scene. In
contrast, our global representation of indoor lighting guarantees temporal
consistency for moving objects.

###  4.4. Training details

Our network is implemented on top of Tensorflow (Abadi et al., 2015). We train
it using SGD and the Adam solver with the moment parameters $\beta_{1}=0.9$
and $\beta_{2}=0.999$. The learning rate is initially set to $0.001$ and
halves every 40 epochs. The weights of all layers in our network are
initialized with the Xavier uniform. Training examples are fed into our
network with a mini-batch size of 5. We train the network for 150 epochs which
takes about two days on one NVIDIA V100 GPU.

Input Image |  (Gardner et al., 2017b) |  (Gardner et al., 2019) |  (Srinivasan et al., 2020) |  (Zhan et al., 2021b) | Ours | Ground Truth
---|---|---|---|---|---|---
|  |  |  |  |  |
|  |  |  |  |  |
| 18.98 dB | 19.19 dB | 17.09 dB | 21.57 dB | 21.74 dB |  $\infty$ dB
|  |  |  |  |  |
|  |  |  |  |  |
| 18.34 dB | 18.36 dB | 24.51 dB | 23.19 dB | 25.60 dB |  $\infty$ dB
Figure 5. Rendering the Armadillo scene (purely diffuse) with the indoor illumination predicted by different methods. Quantitative metrics in terms of PSNR for the rendered images are displayed below each method. Input Image |  (Gardner et al., 2017b) |  (Gardner et al., 2019) |  (Srinivasan et al., 2020) |  (Zhan et al., 2021b) | Ours | Ground Truth
---|---|---|---|---|---|---
|  |  |  |  |  |
|  |  |  |  |  |
| 19.50 dB | 17.88 dB | 18.33 dB | 12.30 dB | 23.37 dB |  $\infty$ dB
|  |  |  |  |  |
|  |  |  |  |  |
| 19.59 dB | 20.09 dB | 20.28 dB | 18.18 dB | 28.61 dB |  $\infty$ dB
Figure 6. Rendering the Dragon scene (specular) with the indoor illumination
predicted by different methods. Quantitative metrics in terms of PSNR for the
rendered images are displayed below each method. Depth values are not
concerned in these cases.

##  5\. Evaluation

In this section, we evaluate the performance of our method in indoor
illumination prediction and make comparisons with the state-of-the-arts. More
results and comparisons can be found in the supplemental material.

###  5.1. Dataset and metrics

To generate our dataset, HDR panoramas are subdivided into two parts, one for
training and the other for testing. Therefore, the panoramas in the test set
will not appear during training. As mentioned in Sec. 3.3, we then randomly
sample FOV-limited images and generate the ground truth DSGLights from the
panoramas. The SUN360 dataset contains 50,000 images for training and the test
set contains 2,000 images. The Laval dataset has 10,000 images for training
and 2,000 for testing. Visual comparisons are conducted on both predicted
environment maps and rendered images with virtual object insertion. PSNR is
adopted as the quantitative measure.

###  5.2. Comparisons to previous methods

We predict indoor illumination with different methods and render two scenes in
Fig. 5 and Fig. 6, respectively. For the methods of (Gardner et al., 2019) and
(Zhan et al., 2021b), we implemented and trained their networks with two
datasets separately. Lighthouse (Srinivasan et al., 2020) requires a held-out
perspective view near the input stereo pair and the depth. Since SUN360 lacks
ground truth depth, we only retrain Lighthouse on the Laval dataset. As for
evaluation on SUN360, we use their pre-trained model. More training details of
Lighthouse are provided in the supplemental material. Rendering is performed
with a physically-based Mitsuba renderer (Nimier-David et al., 2019).

Qualitative evaluation on the Laval dataset: In Fig. 5, we render the
Armadillo scene using the predicted illumination of four recent methods, as
well as the ground truth and our method. The renderings of the Armadillo on a
plane highlight the accuracy of estimated light directions. Note that our
generated shadows and ambient lighting are closer to the ground truth, as
compared with other methods.

Gardner et al. (2017b) proposed a CNN-based pipeline to predict a RGB panorama
and a light mask for the indoor illumination. The predicted panorama is
usually very smooth and may have color drift. Therefore, the synthesized
images are quite different from those generated by the ground-truth lighting.
Gardner et al. (2019) use a parametric form based on SGs to encode indoor
light sources, but the centers of the light sources are allowed to vary
wildly, leading to unstable results. For the method of Srinivasan et al.
(2020), it fails to predict the lighting distribution outside the known
regions, leading to a large bias on rendering. Zhan et al. (2021b) have proved
the effectiveness of SGs by generating accurate illumination maps under the
guidance of predicted SG parameters. However, they do not support spatially-
varying lighting prediction. Most importantly, their CNN networks do not
explore the long-range correlation between SGs with a limited receptive field.

Qualitative evaluation on the SUN360 dataset: In Fig. 6, we further compare
these methods by rendering a specular Dragon model lit by the estimated
lighting. Since the model is highly specular, it can clearly reflect the
environment around it. As expected, our method captures details in the
surroundings (e.g., the white wall and the red floor) quite well. The
predicted panorama is plausible and the details provide better visual
coherence. In comparison, Gardner et al. (2017b) capture the hue of the scene,
but the lighting distribution is overly smooth. In contrast, Zhan et al.
(2021b) hallucinate many details different with the ground truth. For the
method of (Gardner et al., 2019), it hardly infers and represents lighting
details.

Quantitative evaluation: We also provide the quantitative metrics in terms of
PSNR for each method in Table 1. In this table, we list the average PSNR
values for the estimated environment map, the rendered Armadillo scene and the
Dragon scene. Obviously, our method achieves higher PSNR values than others
under every testing scenario. This further shows that our method achieves
state-of-the-art performance on indoor lighting prediction.

Table 1. Quantitative comparisons in terms of PSNR (dB). The best scores are
highlighted in bold. -GCN and -$\mathcal{L}_{W}$ represent the network trained
without GCN module and weighted L2 loss, respectively.

Dataset | Method | Map | Armadillo | Dragon
---|---|---|---|---
Laval |  (Gardner et al., 2017b) | 14.45 | 20.60 | 18.35
(Gardner et al., 2019) | 16.49 | 22.21 | 21.14
(Srinivasan et al., 2020) | 17.50 | 21.42 | 23.00
(Zhan et al., 2021b) | 13.80 | 16.63 | 18.76
Ours | 22.54 | 26.16 | 24.00
-GCN | 20.90 | 24.83 | 22.11
-$\mathcal{L}_{W}$ | 21.87 | 20.30 | 22.61
SUN360 |  (Gardner et al., 2017b) | 11.90 | 23.83 | 20.15
(Gardner et al., 2019) | 12.36 | 23.58 | 19.04
(Srinivasan et al., 2020) | 13.63 | 21.10 | 21.04
(Zhan et al., 2021b) | 9.02 | 20.17 | 15.97
Ours | 13.75 | 24.82 | 21.81
-GCN | 13.16 | 23.39 | 21.21
-$\mathcal{L}_{W}$ | 13.70 | 16.85 | 21.54

GT |  |  |
---|---|---|---
Ours |  |  |
-GCN |  |  |
Figure 7. Validating the contribution of GCN to the performance of the
presented model.

###  5.3. Ablation study

To validate the importance of graph convolutional layers, we remove them from
the whole pipeline (-GCN). Now, the features from feature extraction module
are fed to a fully connected layer which predicts the parameters directly.
Fig. 7 reports the performances of these models. We observe that environment
maps generated by the model without GCN module lack fine details for the
scenes. It can not exchange information across neighboring vertices through
GCN layers and vertices work independently. Therefore, it is hard to infer
illumination from the known region of panoramas. Tab. 1 shows that PSNR drops
slightly when our network removes the GCN module. It also shows that weighted
L2 loss is beneficial for improving the performance.

(a)  (Li et al., 2020)

(b) DSGLight

(c) Ground truth

Figure 8.  Demonstration of spatially-varying lighting on real images of
(Garon et al., 2019).

###  5.4. Validation of spatially-varying lighting

As aforementioned, our method supports spatially-varying lighting when
inserting virtual objects into different locations of an indoor scene, due to
the availability of the depth information. To demonstrate this, we insert
virtual objects at different locations in real images of (Garon et al., 2019).
Two examples are provided in Fig. 8. Here, we compare our method with (Li et
al., 2020) that also enables spatially-varying indoor lighting. As seen, our
results are closer to the ground truth. More importantly, our method
guarantees temporal consistency of lighting when the virtual object moves in
the scene, since our DSGLight is unique and global for a given image. However,
Li et al.’s method (Li et al., 2020) estimates local lighting independently
for each inserted position, which may easily cause temporal flickering.

Figure 9. Each scene in the user study is shown as a column, where different
colors indicate that users preferred the corresponding method instead of the
ground truth.

To further show the superiority of our method, we conduct a perceptual user
study on 20 scenes of (Garon et al., 2019). For each scene, we generate a
reference composite by relighting a bunny model with its ground truth light
probe. We compare these results with objects that were relit with estimated
illumination from (Li et al., 2020) and our method. Users are asked to pick
which rendering is more realistic between image pairs rendered with ground
truth and the estimated lighting. We gather responses from 138 unique
participants and the results are shown in Fig. 9. Both of (Li et al., 2020)
and our method beat each other in half of the scenarios, which suggests that
they are of equal visual quality. However, our method outperforms (Li et al.,
2020) slightly on the average confusion, with a 32.00% vs 30.75% confusion.

###  5.5. Limitations

Despite the state-of-the-art performance, our method still has some
limitations. First, DSGLight is an approximation to the detailed environment
map, lacking sharp features. This prevents specular objects from generating
sharp textures reflecting the surrounding. Moreover, the out-of-view depth
generated by our method is not very accurate due to the high ill-posedness.
This would be solved by multiple input images.

##  6\. Conclusion

In this paper, we have proposed the first graph learning-based strategy for
inferring high-quality and spatially-varying indoor illumination from a single
image. To represent both dominated local lights and detailed environmental
lighting with high fidelity, we encode the panoramic indoor illumination in
DSGLight, a depth-augmented lighting model with 128 SGs evenly distributed
over a sphere. A GCN is adopted to predict non-Euclidean DSGLight from a 2D
image. Through thorough comparisons, we have demonstrated that our method
outperforms previous state-of-the-art methods both qualitatively and
quantitatively.

## References

  * (1)
  * Abadi et al. (2015) Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dan Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. 2015\.  TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems.  http://tensorflow.org/ Software available from tensorflow.org.
  * Atwood and Towsley (2015) James Atwood and Don Towsley. 2015.  Diffusion-Convolutional Neural Networks.  _Arxiv_ (01 2015).
  * Barron and Malik (2013) Jonathan Barron and Jitendra Malik. 2013.  Intrinsic Scene Properties from a Single RGB-D Image.  _Proceedings / CVPR, IEEE Computer Society Conference on Computer Vision and Pattern Recognition. IEEE Computer Society Conference on Computer Vision and Pattern Recognition_ 38 (06 2013).  https://doi.org/10.1109/CVPR.2013.10
  * Cheng et al. (2018) Dachuan Cheng, Jian Shi, Yanyun Chen, Xiaoming Deng, and Xiaopeng. Zhang. 2018.  Learning Scene Illumination by Pairwise Photos from Rear and Front Mobile Cameras.  _Computer Graphics Forum_ 37, 7 (2018), 213–221.
  * Defferrard et al. (2016) Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst. 2016.  Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering. In _Proceedings of the 30th International Conference on Neural Information Processing Systems_ (Barcelona, Spain) _(NIPS’16)_. Curran Associates Inc., Red Hook, NY, USA, 3844–3852.
  * Eilertsen et al. (2017) Gabriel Eilertsen, Joel Kronander, Gyorgy Denes, Rafal K. Mantiuk, and Jonas Unger. 2017\.  HDR image reconstruction from a single exposure using deep CNNs.  _ACM Trans. Graph._ 36, 6 (2017), 178:1–178:15.  https://doi.org/10.1145/3130800.3130816
  * Einabadi et al. (2021) Farshad Einabadi, Jean-Yves Guillemaut, and Adrian Hilton. 2021.  Deep Neural Models for Illumination Estimation and Relighting: A Survey.  _Computer Graphics Forum_ 40, 6 (2021), 315–331.
  * Gardner et al. (2019) Marc-André Gardner, Yannick Hold-Geoffroy, Kalyan Sunkavalli, Christian Gagné, and Jean-François Lalonde. 2019\.  Deep Parametric Indoor Lighting Estimation. In _2019 IEEE/CVF International Conference on Computer Vision, ICCV 2019, Seoul, Korea (South), October 27 - November 2, 2019_. IEEE, 7174–7182.  https://doi.org/10.1109/ICCV.2019.00727
  * Gardner et al. (2017a) Marc-André Gardner, Kalyan Sunkavalli, Ersin Yumer, Xiaohui Shen, Emiliano Gambaretto, Christian Gagné, and Jean-François Lalonde. 2017a.  Learning to predict indoor illumination from a single image.  _ACM Trans. Graph._ 36, 6 (2017), 176:1–176:14.  https://doi.org/10.1145/3130800.3130891
  * Gardner et al. (2017b) Marc-André Gardner, Kalyan Sunkavalli, Ersin Yumer, Xiaohui Shen, Emiliano Gambaretto, Christian Gagné, and Jean-François Lalonde. 2017b.  Learning to Predict Indoor Illumination from a Single Image.  _ACM Trans. Graph._ 36, 6, Article 176 (Nov. 2017), 14 pages.
  * Garon et al. (2019) M. Garon, K. Sunkavalli, S. Hadap, N. Carr, and J. Lalonde. 2019.  Fast Spatially-Varying Indoor Lighting Estimation. In _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_. 6901–6910.  https://doi.org/10.1109/CVPR.2019.00707
  * Green (2003) R. Green. 2003.  Spherical Harmonic Lighting: The Gritty Details.
  * Hamilton et al. (2017) William L. Hamilton, Rex Ying, and Jure Leskovec. 2017\.  Inductive Representation Learning on Large Graphs.  _CoRR_ abs/1706.02216 (2017).  arXiv:1706.02216 http://arxiv.org/abs/1706.02216
  * Hold-Geoffroy et al. (2019) Y. Hold-Geoffroy, A. Athawale, and J. Lalonde. 2019\.  Deep Sky Modeling for Single Image Outdoor Lighting Estimation. In _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_. 6920–6928.  https://doi.org/10.1109/CVPR.2019.00709
  * Hold-Geoffroy et al. (2017) Y. Hold-Geoffroy, K. Sunkavalli, S. Hadap, E. Gambaretto, and J. Lalonde. 2017\.  Deep outdoor illumination estimation. In _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_. 2373–2382.  https://doi.org/10.1109/CVPR.2017.255
  * Hold-Geoffroy et al. (2017) Yannick Hold-Geoffroy, Kalyan Sunkavalli, Sunil Hadap, Emiliano Gambaretto, and Jean-François Lalonde. 2017.  Deep Outdoor Illumination Estimation. 2373–2382.  https://doi.org/10.1109/CVPR.2017.255
  * Karsch et al. (2011) Kevin Karsch, Varsha Hedau, David Forsyth, and Derek Hoiem. 2011.  Rendering Synthetic Objects into Legacy Photographs _(SA ’11)_. Association for Computing Machinery, New York, NY, USA, Article 157, 12 pages.  https://doi.org/10.1145/2024156.2024191
  * Khan et al. (2006) Erum Arif Khan, Erik Reinhard, Roland W. Fleming, and Heinrich H. Bülthoff. 2006.  Image-Based Material Editing _(SIGGRAPH ’06)_. Association for Computing Machinery, New York, NY, USA, 654–663.  https://doi.org/10.1145/1179352.1141937
  * Kipf and Welling (2016) Thomas N. Kipf and Max Welling. 2016.  Semi-Supervised Classification with Graph Convolutional Networks.  _CoRR_ abs/1609.02907 (2016).  arXiv:1609.02907 http://arxiv.org/abs/1609.02907
  * Kipf and Welling (2017) Thomas N. Kipf and Max Welling. 2017.  Semi-Supervised Classification with Graph Convolutional Networks. In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_. OpenReview.net.  https://openreview.net/forum?id=SJU4ayYgl
  * Li et al. (2019) Guohao Li, Matthias Müller, Ali K. Thabet, and Bernard Ghanem. 2019.  DeepGCNs: Can GCNs Go As Deep As CNNs?. In _2019 IEEE/CVF International Conference on Computer Vision, ICCV 2019, Seoul, Korea (South), October 27 - November 2, 2019_. IEEE, 9266–9275.  https://doi.org/10.1109/ICCV.2019.00936
  * Li et al. (2020) Zhengqin Li, Mohammad Shafiei, Ravi Ramamoorthi, Kalyan Sunkavalli, and Manmohan Chandraker. 2020\.  Inverse rendering for complex indoor scenes: Shape, spatially-varying lighting and svbrdf from a single image. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 2475–2484.
  * Liu et al. (2015) Yang Liu, Theo Gevers, and Xueqing Li. 2015.  Estimation of Sunlight Direction Using 3D Object Models.  _IEEE Transactions on Image Processing_ 24, 3 (2015), 932–942.
  * Maier et al. (2017) Robert Maier, Kihwan Kim, Daniel Cremers, Jan Kautz, and Matthias Nießner. 2017.  Intrinsic3D: High-Quality 3D Reconstruction by Joint Appearance and Geometry Optimization with Spatially-Varying Lighting. In _IEEE International Conference on Computer Vision, ICCV 2017, Venice, Italy, October 22-29, 2017_. IEEE Computer Society, 3133–3141.  https://doi.org/10.1109/ICCV.2017.338
  * Micheli (2009) A. Micheli. 2009\.  Neural Network for Graphs: A Contextual Constructive Approach.  _IEEE Transactions on Neural Networks_ 20, 3 (2009), 498–511.  https://doi.org/10.1109/TNN.2008.2010350
  * Niepert et al. (2016) Mathias Niepert, Mohamed Ahmed, and Konstantin Kutzkov. 2016\.  Learning Convolutional Neural Networks for Graphs.  _CoRR_ abs/1605.05273 (2016).  arXiv:1605.05273 http://arxiv.org/abs/1605.05273
  * Nimier-David et al. (2019) Merlin Nimier-David, Delio Vicini, Tizian Zeltner, and Wenzel Jakob. 2019.  Mitsuba 2: A Retargetable Forward and Inverse Renderer.  _Transactions on Graphics (Proceedings of SIGGRAPH Asia)_ 38, 6 (Dec. 2019).  https://doi.org/10.1145/3355089.3356498
  * Shuman et al. (2012) David Shuman, Sunil K. Narang, Pascal Frossard, Antonio Ortega, and Pierre Vandergheynst. 2012\.  The Emerging Field of Signal Processing on Graphs: Extending High-Dimensional Data Analysis to Networks and Other Irregular Domains.  _IEEE Signal Processing Magazine_ 30 (10 2012).  https://doi.org/10.1109/MSP.2012.2235192
  * Simonyan and Zisserman (2015) Karen Simonyan and Andrew Zisserman. 2015.  Very Deep Convolutional Networks for Large-Scale Image Recognition. In _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , Yoshua Bengio and Yann LeCun (Eds.).  http://arxiv.org/abs/1409.1556
  * Song and Funkhouser (2019) Shuran Song and Thomas A. Funkhouser. 2019.  Neural Illumination: Lighting Prediction for Indoor Environments. In _IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2019, Long Beach, CA, USA, June 16-20, 2019_. Computer Vision Foundation / IEEE, 6918–6926.
  * Srinivasan et al. (2020) Pratul P. Srinivasan, Ben Mildenhall, Matthew Tancik, Jonathan T. Barron, Richard Tucker, and Noah Snavely. 2020.  Lighthouse: Predicting Lighting Volumes for Spatially-Coherent Illumination. In _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2020, Seattle, WA, USA, June 13-19, 2020_. IEEE, 8077–8086.  https://doi.org/10.1109/CVPR42600.2020.00810
  * Tsai and Shih (2006a) Yu-Ting Tsai and Zen-Chung Shih. 2006a.  All-frequency precomputed radiance transfer using spherical radial basis functions and clustered tensor approximation.  _ACM Trans. Graph._ 25, 3 (2006), 967–976.  https://doi.org/10.1145/1141911.1141981
  * Tsai and Shih (2006b) Yu Ting Tsai and Zen-Chung Shih. 2006b.  All-frequency precomputed radiance transfer using spherical radial basis functions and clustered tensor approximation.  _ACM Transactions on Graphics_ 25, 3 (1 July 2006), 967–976.  https://doi.org/10.1145/1141911.1141981 null ; Conference date: 30-07-2006 Through 03-08-2006.
  * Wang et al. (2009) Jiaping Wang, Peiran Ren, Minmin Gong, John Snyder, and Baining Guo. 2009.  All-frequency rendering of dynamic, spatially-varying reflectance.  _ACM Trans. Graph._ 28, 5 (2009), 133.  https://doi.org/10.1145/1618452.1618479
  * Weber et al. (2018) Henrique Weber, Donald Prévost, and Jean-François Lalonde. 2018\.  Learning to Estimate Indoor Lighting from 3D Objects. In _International Conference on 3D Vision, Verona, Italy_. 199–207.
  * Wu et al. (2011) C. Wu, B. Wilburn, Y. Matsushita, and C. Theobalt. 2011\.  High-quality shape from multi-view stereo and shading under general illumination. In _CVPR 2011_. 969–976.  https://doi.org/10.1109/CVPR.2011.5995388
  * Xiao et al. (2012) Jianxiong Xiao, Krista A. Ehinger, Aude Oliva, and Antonio Torralba. 2012.  Recognizing scene viewpoint using panoramic place representation. In _2012 IEEE Conference on Computer Vision and Pattern Recognition, Providence, RI, USA, June 16-21, 2012_. IEEE Computer Society, 2695–2702.  https://doi.org/10.1109/CVPR.2012.6247991
  * Yu et al. (2021) Piaopiao Yu, Jie Guo, Fan Huang, Cheng Zhou, Hongwei Che, Xiao Ling, and Yanwen Guo. 2021.  Hierarchical Disentangled Representation Learning for Outdoor Illumination Estimation and Editing. In _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_. 15313–15322.
  * Zhan et al. (2021a) Fangneng Zhan, Changgong Zhang, Wenbo Hu, Shijian Lu, Feiying Ma, Xuansong Xie, and Ling Shao. 2021a.  Sparse Needlets for Lighting Estimation with Spherical Transport Loss. In _Proceedings of the IEEE International Conference on Computer Vision_.
  * Zhan et al. (2021b) Fangneng Zhan, Changgong Zhang, Yingchen Yu, Yuan Chang, Shijian Lu, Feiying Ma, and Xuansong Xie. 2021b.  EMLight: Lighting Estimation via Spherical Distribution Approximation. In _Proceedings of the AAAI Conference on Artificial Intelligence_.
  * Zhang et al. (2019) J. Zhang, K. Sunkavalli, Y. Hold-Geoffroy, S. Hadap, J. Eisenman, and J. Lalonde. 2019.  All-Weather Deep Outdoor Lighting Estimation. In _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_. 10150–10158.  https://doi.org/10.1109/CVPR.2019.01040


# Time Series Clustering for Grouping Products Based on Price and Sales
Patterns

Aysun Bozantaa, Sean Berrya, Mucahit Cevika, Beste Bulutb, Deniz Yigitb,
Fahrettin F. Gonenb, and Ayşe Başara CONTACT M. Cevik Author. Email:
<EMAIL_ADDRESS>aData Science Lab at Ryerson University, Toronto, Canada;
bGetir Perakende Lojistik A.S., Istanbul, Turkey

###### Abstract

Developing technology and changing lifestyles have made online grocery
delivery applications an indispensable part of urban life. Since the beginning
of the COVID-19 pandemic, the demand for such applications has dramatically
increased, creating new competitors that disrupt the market. An increasing
level of competition might prompt companies to frequently restructure their
marketing and product pricing strategies. Therefore, identifying the change
patterns in product prices and sales volumes would provide a competitive
advantage for the companies in the marketplace. In this paper, we investigate
alternative clustering methodologies to group the products based on the price
patterns and sales volumes. We propose a novel distance metric that takes into
account how product prices and sales move together rather than calculating the
distance using numerical values. We compare our approach with traditional
clustering algorithms, which typically rely on generic distance metrics such
as Euclidean distance, and image clustering approaches that aim to group data
by capturing its visual patterns. We evaluate the performances of different
clustering algorithms using our custom evaluation metric as well as Calinski
Harabasz and Davies Bouldin indices, which are commonly used internal validity
metrics. We conduct our numerical study using a propriety price dataset from
an online food and grocery delivery company, and the publicly available
Favorita sales dataset. We find that our proposed clustering approach and
image clustering both perform well for finding the products with similar price
and sales patterns within large datasets.

###### keywords:

time series clustering; image clustering; machine learning; unsupervised
learning; product pricing and sales

##  1 Introduction

E-commerce is one of the fastest-growing channels for global online grocery
delivery, recording a $198.5 billion global market in 2020, and is projected
to reach a market size of $550.7 billion by 2027 (Global Industry Analysts,
2021). The conditions created during the pandemic, the developing technology,
and the changing lifestyles all point to the need for such online platforms.
Thus, the number of companies who have recently participated in the market has
rapidly increased. With the growth of the market and the entry of different
players, dynamic product pricing and sales strategies have become an important
tool for the companies. The strategic and competitive product pricing and
sales might move companies one step ahead of their competitors, and help them
gain a competitive advantage.

Identifying the products with similar price and sales patterns over a given
time period might provide a better understanding of the market, and create
opportunities for better marketing and pricing strategies in the long run. In
this regard, time series clustering can help to obtain the products whose
sales volumes and prices move similarly over a given time period. By examining
the generated clusters, we may also gain insights into the time periods in
which certain products have been discounted or got a price hike.

Existing clustering algorithms commonly use distance metrics that measure the
numerical distance between the data points. These metrics are not capable of
capturing the simultaneous change in the time series patterns. That is, even
if two time series move in opposite directions in a given time period, if they
are numerically close to each other, these metrics will find the distance
between these two time series to be small, and disregard the movement patterns
of the time series. In terms of pricing, this might prevent capturing the
movement patterns that are useful for obtaining the price change points and
directions, and only enable putting time series with close prices in the same
clusters. On the other hand, capturing the price change patterns implies
regardless of the value, if the prices of the products move together in the
same periods (i.e., price drops or price hikes), then they should fall in the
same cluster.

In this study, we investigate clustering methods to group similar grocery
products in terms of price and sales volume change patterns. For this purpose,
we develop a custom distance metric that measures the distance between two
time series based on their movement patterns, which is then used within a
hierarchical clustering algorithm to obtain the product groupings. In
addition, we apply image clustering algorithms that aim to group similar data
instances based on their visual time series patterns. We also experiment with
standard clustering approaches (e.g., $k$-means clustering) by using generic
distance metrics (e.g., Euclidean distance). One of the challenging tasks in
clustering analysis is to measure the performance of the generated clusters.
In particular, when there are no cluster labels (i.e., the ground truth), only
the internal validity metrics can be used to evaluate the quality of the
clusters, which rely on generic distance metrics. However, since the generic
metrics such as Euclidean and Manhattan distances only consider the numerical
closeness of the data points within the clusters, they are not suitable for
our problem, which prompted us to develop a custom evaluation metric.

The contribution of our study can be summarized as follows.

  * •

Our study provides an application of time series clustering for grocery
product grouping based on price and sales volume history, which is an
important practical problem for marketing and sales.

  * •

To the best of our knowledge, image clustering has not been previously used
for the task of clustering products based on price and sales volumes. We find
this approach to be effective in capturing common patterns between different
time series, and hence grouping the grocery products using the price and sales
histories.

  * •

We propose a custom distance metric that measures the distance between two
time series based on their movement patterns. This metric can also be used for
time series clustering tasks other than grocery product grouping. We also
develop a custom evaluation metric to evaluate the performance of the
clustering algorithms, which scores the clusters using our custom distance
metric.

The remainder of the paper is organized as follows. Section 2 provides an
overview of clustering approaches with a particular focus on product
clustering with different techniques. Section 3 provides a summary of the
methodologies used, dataset characteristics, and a detailed discussion of our
preprocessing, algorithms, and evaluation metrics. In Section 4, we present
the results of our numerical study, which includes an illustrative example of
how our custom distance metric works, and a performance comparison of
clustering algorithms with different distance metrics by using various
evaluation criteria. Lastly, Section 5 provides concluding remarks and future
research directions.

##  2 Literature Review

Time series clustering has been used in diverse scientific domains to discover
hidden patterns from complex and massive datasets. Clustering methodology in
time series analysis can be examined with respect to the data representation,
distance metric, and clustering algorithm (Aghabozorgi et al., 2015). The
application of data transformation techniques on the original data may yield
better results in clustering algorithms. Commonly used data representation
techniques in the literature are Discrete Fourier Transformation (DFT),
Singular Value Decomposition (SVD), and Discrete Wavelet Transform (DWT)
(Aghabozorgi et al., 2015). Another component of clustering is the distance
metric that is used within the clustering algorithm. Well-known distance
metrics include Euclidean distance, Dynamic Time Warping (DTW), and Longest
Common Subsequence (LCSS) (Aghabozorgi et al., 2015). The clustering
approaches are mainly classified into six groups, namely, partitioning,
hierarchical, grid-based, model-based, density-based clustering, and multi-
step clustering. On the other hand, for time series analysis, mostly
partitioning and hierarchical clustering techniques have been employed.

Martinez Alvarez et al. (2011) applied time series forecasting on the
electricity demand dataset. Before the time series forecasting, they used
$k$-means clustering to group and label the time series data. They evaluated
the performance of the clustering algorithm by three well-known internal
validity indices, namely, silhouette index, Davies-Bouldin (DB) index, and the
Dunn (DU) index. Bandara et al. (2020) proposed a novel forecasting model
which can be used with different types of RNN models on subgroups of similar
time series. Their clustering technique is feature-based, which, instead of
capturing similarity of point values using a distance metric, employs sets of
global features obtained from a time series to summarize and describe the
salient information of the time series. They noted that feature-based
approaches can be more interpretable and more resilient to missing and noisy
data. Javed et al. (2020) used the time series datasets for the University of
California Riverside (UCR) archive and conducted a benchmark time series
clustering study. They combined the $k$-means, fuzzy $c$-means, agglomerative
clustering, and density-based clustering algorithms with the Euclidean, DTW,
and shape-based distance metrics, and compared alternative clustering
approaches using the rand index.

Unlike other clustering approaches, image clustering determines the visual
patterns of the data and groups similar patterns together. This method
consists of two components: feature learning and clustering. For the feature
learning phase, deep learning techniques have been most commonly used in the
literature (Wu et al., 2019; Ren et al., 2020), whereas, for the clustering
phase, $k$-means clustering is typically employed (Caron et al., 2018). Ahmed
(2015) applied different image clustering algorithms on 14 benchmark image
datasets and observed that the performance of image clustering algorithms
depends mostly on the distribution of the image datasets. The application of
image clustering techniques on the time series data is very limited, and they
have been mostly used on the satellite image time series datasets (Khiali et
al., 2019; Lampert et al., 2019).

Time series modeling is frequently used in sales and marketing, especially for
forecasting future sales and demand (Pavlyshenko, 2019). Accurate forecasting
of future sales and demand allows companies to leverage information from
historical data strategically to increase future revenue (Leonard, 2000).
There is an extensive literature on time series models for sales and demand
forecasting. Abbasimehr et al. (2020) developed multi-layer LSTM networks and
automated hyperparameter tuning approaches for demand forecasting, and showed
their to be effective when compared to some other popular models such as ARIMA
and support vector machines. Pavlyshenko (2019) compared different machine
learning approaches for sales forecasting problem, and found that a model
stacking approach could reliably improve model performance to a significant
degree. For more information on time series models for sales and demand
forecasting, we refer the reader to a recent review by Fildes et al. (2019).
Despite the popularity of time-series forecasting over sales and demand
datasets, there has been very little research done in the area of time series
clustering for these types of data.

Although the time series clustering approaches have been considered in various
domains such as climate (Awad and Qasrawi, 2018), biology (McDowell et al.,
2018), finance (Takahashi et al., 2019), environment and urban planning (Fu et
al., 2019), the applications of time series clustering techniques on product
clustering is very limited. In addition, image clustering approaches have not
been commonly used on time series data, and to the best of our knowledge,
previous studies did not consider clustering approaches for product grouping
based on price and sales volume histories.

##  3 Methodology

In this section, we discuss the datasets, our proposed distance metric for
time series data, the image clustering approach, and the evaluation metrics
for measuring clustering performance. Figure 1 presents the overall
methodology of this study.

Figure 1: Overall methodology for grouping products based on their price
histories, which are represented as distinct time series for each product
(similar steps are followed for product clustering based on sales volume).

###  3.1 Datasets

We consider two datasets in our numerical analysis. The first dataset consists
of a collection of price values over 3,285 products for a period of $n=570$
days. This dataset is obtained from XYZ company111Blinded for peer review an
online food and grocery delivery company that originated in Turkey and
recently expanded its operations to the United Kingdom and the Netherlands.
Daily price values for the products are available as numerical values in the
dataset. The second dataset is the publicly available Favorita sales
data222https://www.kaggle.com/c/favorita-grocery-sales-forecasting/data. The
Favorita dataset was provided by Corporacion Favorita, an Ecuadorian-based
grocery retailer, with the aim of predicting future sales of their products.
Favorita sales data contains 59,038,132 observations over 54 stores, taken
between January 1st 2016 and August 15th 2017. Each observation includes a
date, a store number, item number, the unit sales for the day, as well as
whether the item had a promotional price that day. After aggregating the
distinct products and stores, we end up with 172,130 distinct product store
pairs.

We apply a series of data preprocessing steps for both datasets. We format the
raw data as time series by aggregating all the price/sales values for each
product over time. Real-world time series data is often rife with aberrations.
We remove any product with over 80% of their prices/sales unavailable. For the
remaining products, many of which have missing price values, we fill the
missing prices with the most recently available price for that product. For
Favorita dataset, we fill the missing sales with the mean values.

We next use a min-max scaler to scale the prices/sales of each
product/product-store pair to be between 0.1 and 1. This approach is helpful
since we are interested in relative price/sales fluctuations, and products
with different prices/sales but similar fluctuations are aimed to be grouped
together. After proper cleanup, scaling, and exploratory data analysis, we
note that many of the time series are not suitable to be grouped in any
clusters. We examine the pairwise distances between time series and remove
those that are not likely to be clustered together with any other time series
in the dataset. Sample time series for price and sales data is provided in
Figure LABEL:fig:sample_price and Figure LABEL:fig:sample_sale, respectively.
We note that, while price data shows little fluctuations over time, as
expected, sales values are highly variable.

(a)

(b)

Figure 2: Time series samples.

At the end of the preprocessing phases, 1,160 product price time series
remained in the price dataset, while 134 product-store sales time series
remained in the sales dataset. Those are then discretized for our custom
metric as well as any string-based clustering method that focuses on movement
patterns rather than the actual numerical values. Specifically, we use the
following function that considers equal distances between the normalized
values, which are identified based on preliminary analysis.

$\displaystyle z(x)=\begin{cases}A&x<0.29\\\ B&0.29\leq x<0.47\\\ C&0.47\leq
x<0.65\\\ D&0.65\leq x<0.83\\\ E&0.83\leq x\end{cases}$ (1)

###  3.2 Distance Metrics

We employ three distance metrics, namely, Euclidean Distance, Levenshtein
Distance, and Movement Pattern-based Distance, which we briefly summarize
next.

####  3.2.1 Euclidean Distance

Let $p$ and $q$ be two product time series from our dataset. Then the
Euclidean distance between these products is given by
$D(p,q)=\sqrt{\sum_{i=1}^{n}(q_{i}-p_{i})^{2}}$.

####  3.2.2 Levenshtein Distance

This is a string metric that can be simply defined as the minimum number of
single-character edits required to change one word into the other. In our
case, this corresponds to the number of single-character edits required to
change the discretized prices/sales of one product into the other. Let $p$ and
$q$ be two product time series from our dataset. Then the Levenshtein distance
between these products is given by using the below recursive function:

$\displaystyle D(p,q)=\begin{cases}|p|&\text{if}\ |q|=0\\\\[1.99997pt]
|q|&\text{if}\ |p|=0\\\\[1.99997pt] D(p[1:],q[1:])&\text{if}\
p[0]=q[0]\\\\[1.99997pt] 1+\min\begin{cases}D(p[1:],q[1:])&\\\\[1.99997pt]
D(p,q[1:])&\\\\[1.99997pt] D(p[1:],q)&\\\\[1.99997pt]
\end{cases}&\text{otherwise}\end{cases}$

####  3.2.3 Movement Pattern-based Distance

We propose a custom distance metric, the movement pattern-based distance
(MPBD), which takes into account the direction of change in a time series.
MPBD is developed by iterating through two time series simultaneously and
following the below rules.

  * •

If both time series have constant values at the same time, no distance is
added.

  * •

If both time series move in the same direction, and the same magnitude at the
same time no distance is added.

  * •

If both time series move in the same direction, but at a different magnitude,
absolute value of the difference between the magnitudes is added.

  * •

If both time series move in the opposite directions, a distance value
corresponding to the absolute value of the difference between the magnitudes
multiplied by some weight (e.g., $\omega=2$) is added.

In this manner, we aim to capture products (i.e., time series) that move
similarly even with different price points or relative price points. Let
$d_{p}=p[0]-p[1]$ for product $p$ in the dataset. A naive recursive definition
of the distance function for two products $p$ and $q$ can be obtained as
follows.

$\displaystyle\begin{cases}0&\text{if }|p|=|q|=1\\\ 0+D(p[1:],q[1:])&\text{if
}d_{p}\And d_{q}=0\\\ 0+D(p[1:],q[1:])&\text{if }d_{p}=d_{q}\\\
|d_{p}-d_{q}|+D(p[1:],q[1:])&\text{if }d_{p}\neq d_{q}\And\\\
&sign(d_{p})=sign(d_{q})\\\ \omega\cdot|d_{p}-d_{q}|+D(p[1:],q[1:])&\text{if
}sign(d_{p})\neq sign(d_{q})\\\ \end{cases}$

###  3.3 Clustering Algorithms

We consider two traditional clustering algorithms ($k$-means and hierarchical
clustering) as well as image clustering for grouping the product time series.

####  3.3.1 $k$-means Clustering

$k$-means is a popular clustering algorithm, which aims to minimize the
within-cluster sum of squares (WCSS) measure, that is,

$WCSS=\sum_{i=1}^{k}\sum_{j=1}^{n}(x_{ij}-\mu_{i})^{2}$ (2)

where $x_{ij}$ is item $j$ in cluster $i$ and $\mu_{i}$ is the corresponding
cluster centre of cluster $i$. Items are retroactively assigned to their
nearest cluster center, and then the cluster centers are updated until the
items in the clusters no longer change.

####  3.3.2 Hierarchical Clustering

Hierarchical clustering requires either starting with all items in one cluster
and then breaking that cluster into smaller sub-clusters (divisive
hierarchical), or the exact opposite, i.e., starting with all items in their
own cluster and joining clusters based on some criteria (agglomerative
hierarchical). In our implementations, we use an agglomerative hierarchical
approach with Ward’s linkage criteria.

####  3.3.3 Image Clustering

Image clustering is an increasingly popular research area with many potential
applications. We consider image clustering for time series data obtained over
product price/sale histories. Specifically, for each product in the dataset,
we create an image illustrating the graph of the numerical values over time.
We transform these images into array representations in three-dimensional
space. We then use the VGG16 convolutional neural network (CNN) to extract
feature vectors. Finally, we apply the $k$-means clustering algorithm on the
feature vectors of the time series images.

###  3.4 Evaluation Metrics

Since our dataset does not include ground truth labels (i.e., actual cluster
ids), we rely on internal indices for clustering performance assessment, which
is briefly summarized below.

####  3.4.1 Calinski Harabasz (CH) Index

This index is defined as the ratio of the sum of between-clusters dispersion
and inter-cluster dispersion. Thus a higher score implies that the clusters
are very distinct from each other yet the cluster contents are dense. We
define the between cluster sum of squares as follows.

$\displaystyle BCSS=\sum_{i=1}^{k}(\mu_{i}-\mu)^{2}$ (3)

where $\mu_{i}$ is the cluster centre for cluster $i$ and $\mu$ is the sample
mean for the whole dataset. Following from Eqn. 2 and Eqn. 3, the CH index can
be obtained as $CH=\frac{WCSS}{BCSS}$.

####  3.4.2 Davies Bouldin (DB) Index

This index is another popular internal validity metric, which aims to measure
the average similarity of each cluster with a cluster most similar to it. A
lower DB value is an indicator that the clusters are well separated. The
calculation of the DB index is most easily represented in multiple steps as
follows.

  * •

First, we calculate the intra-cluster dispersion for each cluster with the
formula $S_{i}=\sqrt{\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}(x_{ij}-\mu_{i})^{2}}$,
where $n_{i}$ is the number of items in cluster $i$ and $\mu_{i}$ is the
cluster center of cluster $i$.

  * •

Next, we calculate the distance between cluster centers for all clusters with
the formula: $M_{i,j}=\sqrt{(\mu_{i}-\mu_{j})^{2}}$

  * •

The similarity measure between two clusters $i$ and $j$ is thus
$R_{i,j}=\frac{S_{i}+S_{j}}{M_{i,j}}$

  * •

DB index is then given by $DB=\frac{1}{k}\sum_{i=1}^{k}\max_{j}(R_{i,j})$.

####  3.4.3 Movement Pattern-based Index

Let $D_{cust}(x,j)$ be the distance between two products $x$ and $j$ as
defined by our custom distance metric, MPBD. Then a custom evaluation metric,
which we refer to as movement pattern-based index (MPBI) can be obtained as
follows

$\displaystyle
MPBI=\frac{1}{k}\sum_{c=1}^{k}\frac{1}{n}\Big{[}\sum_{i=0}^{n-1}\sum_{j=i+1}^{n}D_{cust}(x_{ci},x_{cj})\Big{]}$

where $x_{{ci}}$ denotes the $i$th product in cluster $c$, $n$ denotes the
number of data points in a given cluster, and $k$ is the number of clusters.
Specifically, this measure is the sum of the clusters’ pairwise distances
divided by the number of products in the same cluster and averaged over all
the clusters. Therefore a lower MPBI value suggests low pairwise distances,
which in turn suggest that the clustered items are close. Additionally, this
index allows us to use the discretized data exclusively without the need for a
numeric dataset.

##  4 Results

In this section, first, the working principle of MPBD is presented in an
illustrative example. Then, the experimental results comparing the
performances of the combination of two clustering algorithms with four
different distance metrics are provided for a clustering task involving
product price and product sales time series. Lastly, we provide cluster
profile results to show which products tend to be grouped together based on
their price and sales histories.

###  4.1 Illustrative Example

We created three scenarios over synthetic price time series to depict the
difference between MPBD, Euclidean distance, and Levenshtein distance. Figure
3 presents the price changes of two products over a ten-day period for three
scenarios. Table 1 shows the normalized distance values between two time
series representing product prices, calculated by using different distance
metrics. In scenario 1, although the price levels of the two products are
entirely different from each other, they show the same price change pattern.
Since MPBD considers the price change similarity instead of the exact
numerical closeness of two prices, its value becomes 0.0 for scenario 1. On
the other hand, because each price value in each time step is different from
each other for scenarios 1 and 2, Levenshtein distance values become 1.0,
which is the maximum value attained for these scenarios. Euclidean distance
for scenario 1 is 0.63 because, although the prices of these two products are
numerically different from each other, their normalized distances are
relatively close (see Eqn. (1)).

(a)

(b)

(c)

Figure 3:  Time series price changes of two products in different scenarios Table 1: Comparison of the distance metrics based on three different scenarios Distance Metric | scenario 1 | scenario 2 | scenario 3
---|---|---|---
MPBD | 0.00 | 0.10 | 0.60
Levenshtein | 1.00 | 1.00 | 0.30
Euclidean | 0.63 | 0.58 | 0.35

Scenarios 1 and 2 are very similar, but in scenario 2, product0 price drops by
one unit (i.e., string distance between B and A) on day 4, while product1
price drops by two units, and, on day 7, the product0 price increases by one
unit and product1 price increases by two units. Although they move in exactly
the same direction during the 10-day period, the magnitudes of the price
movements are different from each other on the 4th and 7th days. According to
the MPBD, they are still very close to each other, however, since MPBD also
considers the magnitudes of the movements, it results in a distance value of
0.1. When compared to scenario 1, since prices get closer numerically on days
4 and 7, the Euclidean distance becomes 0.58 for scenario 2.

In scenario 3, the product prices move in opposite directions on the 3rd, 6th,
and 8th days. Except for these three days, prices move steadily and show the
same movement pattern. Besides, for a 5-day period, the prices continue at the
same level. Accordingly, when measured with the Levenshtein and Euclidean
distance measures, the price difference between these two products seems to be
very small, while the MPBD takes the highest value among these three
scenarios. This is because 3 out of 10 days, the prices of these two products
move in opposite directions.

###  4.2 Performance of the Clustering Algorithms

We evaluate the performances of hierarchical and $k$-means clustering
algorithms with the combination of different distance metrics by using CH, DB,
and the MBPI. For these experiments, we use normalized versions of the
distance values generated by the different distance metrics to avoid any bias
due to the scale of the metrics. We first run the hierarchical clustering
algorithm with the MPBD for a different number of clusters ranging from 0 to
100, and plot the graph of CH and DB scores (see Figure 4). After 15 clusters,
the CH score decreases significantly for the price dataset. Although one of
the lowest DB values seems to be obtained with 13 clusters, the most
appropriate number of clusters is chosen as 15 for the price dataset, since a
very similar value was obtained with 15 clusters. For the sales dataset, after
7 clusters, the CH score decreases significantly. Although DB has lower values
for a higher number of clusters, at 7 clusters, there is a sharp drop in its
value. As a result, we consider $k=15$ clusters for price dataset and $k=7$
clusters for sales dataset in the remaining experiments.

(a)

(b)

(c)

(d)

Figure 4: Results of hierarchical clustering with MPBD

Table 2 presents the performances of $k$-means and hierarchical clustering
algorithms with respect to DB, CH, and MBPI. For the price dataset, we obtain
the highest CH and the lowest DB score by using a $k$-means algorithm with the
Euclidean distance. This is expected because the evaluation criterion for
these two scores is the Euclidean distance, and they consider the numerical
closeness of the product prices. This algorithm is followed by a hierarchical
clustering algorithm using Euclidean distance and Levenshtein distance.
However, when we look at the MBPI, we see that the hierarchical clustering
with the MPBD and image clustering ($k$-means with CNN-based feature learning)
algorithms perform very similarly, and MBPI values are lower compared to other
algorithms. That is, both hierarchical clustering with the MPBD and image
clustering techniques provide better clusters than other algorithms based on
the MBPI, as they both aim to find price change pattern similarities of
products (i.e., time series), and generate the clusters accordingly. For the
sales dataset, the image clustering algorithm, which uses the CNN-based
algorithm for feature extraction and $k$-means for clustering, obtains the
best values for the DB and CH evaluation metrics, while $k$-means with
Euclidean distance obtains the best value for the MBPI metric. Our clustering
approach obtains comparable results with the image clustering approach in
terms of the MBPI metric.

Table 2: Summary results for the clustering algorithms and distance metrics

|  | Price Dataset | Sales Dataset
---|---|---|---
Clustering Alg. | Distance Metric | DB | CH | MBPI | DB | CH | MBPI
$k$-means | Euclidean | 1.7 | 198.4 | 1525.3 | 2.6 | 8.2 | 7046.5
$k$-means | CNN-based | 4.1 | 63.5 | 1494.7 | 2.3 | 13.7 | 7394.8
$k$-means | DTW | 6.4 | 49.2 | 2259.0 | 4.1 | 6.2 | 7481.6
Hierarchical | Euclidean | 2.8 | 178.8 | 1551.5 | 3.9 | 7.2 | 7175.3
Hierarchical | Levenshtein | 3.8 | 149.2 | 1557.1 | 4.1 | 6.7 | 7362.5
Hierarchical | DTW | 51.8 | 4.7 | 2269.0 | 6.1 | 3.6 | 7406.2
Hierarchical | MPBD | 5.8 | 30.4 | 1415.0 | 4.3 | 5.2 | 7125.0

For the price dataset, the number of time series in the clusters that are
generated by both image clustering and the custom clustering approaches are
not as homogeneous as the distribution of time series among the clusters
generated by other algorithms. In this dataset, there are some products whose
prices do not change or rarely change over the given time period. Because the
price values are more stable compared to sales values, number of such products
can be quite large and their prices typically differ from each other. Although
other traditional clustering techniques cannot capture such patterns, both
image clustering and our custom clustering approach captures this pattern, and
place these products in the same cluster. Therefore, in the clusters that are
created by both image clustering and custom clustering, there are two clusters
having more than 200 time series consisting of these products. Since the sales
values fluctuate at a higher rate, it is more difficult to capture above
mentioned pattern. However, we find that the image clustering approach
generates very homogeneous clusters in terms of the number of time series for
the sales dataset.

Figure 5 presents time series in two sample clusters generated by the custom
clustering approach for the price dataset. As well as the custom clustering
approach captures the time series showing almost the same pattern in the same
price level, it also captures the similar patterns that occur in different
price levels. Both in clusters 1 and 5, although there are some overlapping
patterns in the same level (especially in cluster 5), there are similar
patterns in different price levels.

(a)

(b)

Figure 5:  Two sample clusters generated by the hierarchical clustering with
MPBD for the price dataset.

Figure 6 presents two sample clusters generated using the $k$-means algorithm
with Euclidean distance for the sales dataset. We note that this approach
results in time series clusters that are close together spatially, but may not
capture the similarity of the underlying movement patterns.

(a)

(b)

Figure 6:  Two sample clusters generated by $k$-means clustering with
Euclidean distance for the sales dataset.

###  4.3 Cluster Profiling

We create the cluster profiles based on the available features of both
datasets to better understand the cluster contents. Table 3 presents the
cluster-ID, the number of time series, the number of product categories, top
two categories having the highest number of products, average, min, and max
prices of each cluster that are generated by hierarchical clustering with MPBD
metric using the price dataset. The minimum number of categories in clusters
is one and the maximum number of categories is 15. However, the top two
categories cover almost 50% of all products in each cluster. The average
product prices in the clusters vary between 9 and 34, while the minimum prices
of the products are between 1 and 10, and finally, the maximum price values
change between 6 and 120. Some products in the Personal Care category and the
entire Pet category are separated from other products, and form homogeneous
clusters. Other categories that are often in the same cluster are (Snacks,
Water & Beverage), (Snacks, Essential Food) and (Essential Food, Water &
Beverage).

Table 3: Cluster profiling for the price dataset

cluster | size | # of categories | top 2 categories (# of products) | avg. price | min price | max price
---|---|---|---|---|---|---
1 | 9 | 1 | Personal Care: 9 | 34 | 10 | 80
2 | 8 | 1 | Personal Care: 8 | 18 | 8 | 26
3 | 9 | 1 | Pet: 9 | 4 | 3 | 6
4 | 20 | 4 | Home & Life: 7 | 14 | 1 | 30
|  |  | Home Maintenance: 7 |  |  |
5 | 68 | 8 | Personal Care: 18 | 19 | 2 | 63
|  |  | Water & Beverages: 15 |  |  |
6 | 23 | 2 | Home Maintenance: 5 | 19 | 2 | 42
|  |  | Personal Care: 18 |  |  |
7 | 56 | 10 | Personal Care: 12 | 18 | 3 | 120
|  |  | Home & Life: 10 |  |  |
8 | 82 | 8 | Essential Food: 21 | 14 | 1 | 78
|  |  | Personal Care: 20 |  |  |
9 | 22 | 6 | Baby: 6 | 15 | 4 | 40
|  |  | Water & Beverages: 5 |  |  |
10 | 69 | 8 | Essential Food: 27 | 9 | 2 | 37
|  |  | Water & Beverages: 12, Baby: 12 |  |  |
11 | 71 | 11 | Snacks: 20 | 11 | 1 | 40
|  |  | Essential Food: 15 |  |  |
12 | 261 | 12 | Water & Beverages: 57 | 11 | 1 | 110
|  |  | Snacks: 55 |  |  |
13 | 234 | 15 | Essential Food: 46, Milk and Breakfast: 40 | 12 | 1 | 90
|  |  | Snacks: 47 |  |  |
14 | 114 | 13 | Essential Food: 30 | 12 | 1 | 87
|  |  | Water & Beverages: 23 |  |  |
15 | 80 | 10 | Essential Food: 22 | 13 | 1 | 55

When the clusters in Table 3 are examined closely, we find that the products
whose prices do not change or rarely change are ice creams, books and
magazines, frozen and canned goods, nuts, sauces and condiments, and jams.
Another interesting finding is that clusters generated by standard clustering
algorithms and distance measures do not reveal the relationships between
product categories and the price change patterns. However, image clustering
and hierarchical clustering with MPBD tend to generate clusters that are
consisting of a single product category, such as personal care and pet
products. In addition, we observe that the products that are likely to be sold
together (e.g., cereal, avocado, bowl) show similar price change patterns, and
are grouped under the same cluster.

Table 4 presents the similar statistics to those of Table 3 for the sales
dataset, which are obtained by using $k$-means clustering approach (i.e., the
best performing clustering method for this dataset). Note that this dataset
does not include detailed information (e.g., product category names), thereby
limiting this analysis. The number of products in each cluster is fewer than
the number of product-store sales time series. This means that the same
product even if it is sold in different stores shows a similar sales pattern.
Since the number of categories in the data set after preprocessing step is
very small, we cannot make a detailed inference as we did for the price data
set. We observe that cluster four is the cluster with high sales, while
clusters one, three, and seven are low-sales clusters, and clusters two, five
and six are moderate-level sales clusters.

Table 4: Cluster profiling for the sales dataset

cluster | size | # of products | # of stores | # of categories | avg. sales | min sales | max sales
---|---|---|---|---|---|---|---
1 | 49 | 24 | 9 | 1 | 5 | 0 | 53
2 | 22 | 13 | 5 | 1 | 9 | 0 | 80
3 | 22 | 15 | 7 | 2 | 5 | 0 | 58
4 | 19 | 14 | 3 | 1 | 19 | 0 | 109
5 | 16 | 15 | 5 | 2 | 9 | 0 | 67
6 | 5 | 8 | 1 | 1 | 11 | 0 | 39
7 | 1 | 1 | 1 | 1 | 4 | 0 | 21

##  5 Conclusion

In this study, we focus on clustering grocery products based on their price
and sales volume histories, which can be represented as time series. For the
price dataset, this analysis can be used to gather insights about the products
having price hikes and discounts in the same periods. On the other hand, for
the sales dataset, finding the products with similar sales patterns can be
used for various marketing operations and analysis, including which products
to use in promotions and market basket analysis. We develop a custom distance
(MPBD) and also a custom evaluation metric (MBPI) for this purpose. We also
employ the image clustering technique to group the products based on their
price change patterns. We compare these approaches with the traditional
clustering methods and evaluate their performances by using CH, DB, and the
MBPI metrics. The results show that the hierarchical clustering with the MPBD
and the image clustering algorithm achieve comparable results on the time
series product price dataset, and outperform traditional techniques, whereas
high level of fluctuations and low number of time series in the sales dataset
largely prevent having a dominant clustering technique. Overall, we find both
our proposed method and the image clustering as suitable approaches for
clustering over grocery product datasets to identify similar time series based
on their movement patterns. The work presented in this paper can be further
extended by applying both custom and image clustering approaches on different
time series datasets. These clustering approaches could also benefit from
application to similar datasets and tasks to strengthen our findings.
Furthermore, we note that many global forecasting models may be able to
successfully leverage the clustering data generated by the proposed methods in
this study, e.g., for the purposes of increasing predictive power for
forecasting.

## Acknowledgements

The authors would like to thank the XYZ company for supporting this study and
providing data and feedback throughout.

## Data availability statement

The Favorita dataset that supports the findings of this study is available at
the following link:
https://www.kaggle.com/c/favorita-grocery-sales-forecasting/data
Price dataset is propriety and obtained from XYZ company.

## Disclosure statement

No potential conflict of interest was reported by the authors.

## References

  * Abbasimehr et al. (2020) Hossein Abbasimehr, Mostafa Shabani, and Mohsen Yousefi.  An optimized model using lstm network for demand forecasting.  _Computers & Industrial Engineering_, 143:106435, 2020\.  ISSN 0360-8352.  https://doi.org/10.1016/j.cie.2020.106435.
  * Aghabozorgi et al. (2015) Saeed Aghabozorgi, Ali Seyed Shirkhorshidi, and Teh Ying Wah.  Time-series clustering–a decade review.  _Information Systems_ , 53:16–38, 2015.
  * Ahmed (2015) Nasir Ahmed.  Recent review on image clustering.  _IET Image Processing_ , 9(11):1020–1032, 2015\.
  * Awad and Qasrawi (2018) Mohammed Awad and Ibrahim Qasrawi.  Enhanced RBF neural network model for time series prediction of solar cells panel depending on climate conditions (temperature and irradiance).  _Neural Computing and Applications_ , 30(6):1757–1768, 2018.
  * Bandara et al. (2020) Kasun Bandara, Christoph Bergmeir, and Slawek Smyl.  Forecasting across time series databases using recurrent neural networks on groups of similar series: A clustering approach.  _Expert systems with applications_ , 140:112896, 2020.
  * Caron et al. (2018) Mathilde Caron, Piotr Bojanowski, Armand Joulin, and Matthijs Douze.  Deep clustering for unsupervised learning of visual features.  In _Proceedings of the ECCV_ , pages 132–149, 2018.
  * Fildes et al. (2019) Robert Fildes, Shaohui Ma, and Stephan Kolassa.  Retail forecasting: Research and practice.  _International Journal of Forecasting_ , 2019.
  * Fu et al. (2019) Yingchun Fu, Jiufeng Li, Qihao Weng, Qiming Zheng, Le Li, Shu Dai, and Biyun Guo.  Characterizing the spatial pattern of annual urban growth by using time series Landsat imagery.  _Science of The Total Environment_ , 666:274–284, 2019\.
  * Global Industry Analysts (2021) Inc Global Industry Analysts.  Online grocery - global market trajectory & analytics, 2021.
  * Javed et al. (2020) Ali Javed, Byung Suk Lee, and Donna M Rizzo.  A benchmark study on time series clustering.  _Machine Learning with Applications_ , 1:100001, 2020.
  * Khiali et al. (2019) Lynda Khiali, Mamoudou Ndiath, Samuel Alleaume, Dino Ienco, Kenji Ose, and Maguelonne Teisseire.  Detection of spatio-temporal evolutions on multi-annual satellite image time series: A clustering based approach.  _International Journal of Applied Earth Observation and Geoinformation_ , 74:103–119, 2019.
  * Lampert et al. (2019) Thomas Lampert, Baptiste Lafabregue, and Pierre Gançarski.  Constrained distance based k-means clustering for satellite image time-series.  In _IGARSS 2019_ , pages 2419–2422. IEEE, 2019.
  * Leonard (2000) Michael Leonard.  Promotional analysis and forecasting for demand planning: A practical time series approach.  Jan 2000.
  * Martinez Alvarez et al. (2011) Francisco Martinez Alvarez, Alicia Troncoso, Jose C. Riquelme, and Jesus S. Aguilar Ruiz.  Energy time series forecasting based on pattern sequence similarity.  _IEEE TKDE_ , 23(8):1230–1243, 2011.  10.1109/TKDE.2010.227.
  * McDowell et al. (2018) Ian C McDowell, Dinesh Manandhar, Christopher M Vockley, Amy K Schmid, Timothy E Reddy, and Barbara E Engelhardt.  Clustering gene expression time series data using an infinite gaussian process mixture model.  _PLoS Computational Biology_ , 14(1):e1005896, 2018.
  * Pavlyshenko (2019) Bohdan M. Pavlyshenko.  Machine-learning models for sales time series forecasting.  _Data_ , 4(1), 2019.  ISSN 2306-5729.  10.3390/data4010015.  URL https://www.mdpi.com/2306-5729/4/1/15.
  * Ren et al. (2020) Yazhou Ren, Ni Wang, Mingxia Li, and Zenglin Xu.  Deep density-based image clustering.  _Knowledge-Based Systems_ , 197:105841, 2020.
  * Takahashi et al. (2019) Shuntaro Takahashi, Yu Chen, and Kumiko Tanaka-Ishii.  Modeling financial time-series with generative adversarial networks.  _Physica A: Statistical Mechanics and its Applications_ , 527:121261, 2019.
  * Wu et al. (2019) Jianlong Wu, Keyu Long, Fei Wang, Chen Qian, Cheng Li, Zhouchen Lin, and Hongbin Zha.  Deep comprehensive correlation mining for image clustering.  In _Proceedings of the IEEE/CVF ICCV_ , pages 8150–8159, 2019.