# Characterization of Frequent Online Shoppers using Statistical Learning with
Sparsity

Rajiv Sambasivan, Jörg Schad, Arthur Keen, Christopher Woodward, Alexander
Geenen, Sachin Sharma
ArangoDB Inc
San Mateo, CA
{rajiv, joerg, arthur, christopher, alex<EMAIL_ADDRESS>
&Mark Burgess
Chitek-i
<EMAIL_ADDRESS>

###### Abstract

Developing shopping experiences that delight the customer requires businesses
to understand customer taste. This work reports a method to learn the shopping
preferences of frequent shoppers to an online gift store by combining ideas
from retail analytics and statistical learning with sparsity. Shopping
activity is represented as a bipartite graph. This graph is refined by
applying sparsity-based statistical learning methods. These methods are
interpretable and reveal insights about customers’ preferences as well as
products driving revenue to the store.

_K_ eywords Retail Analytics $\cdot$ Matrix Factorization $\cdot$ Statistical
Learning with Sparsity

##  1 Introduction

Understanding customers’ preferences is important in the retail industry [1].
This understanding is used to personalize the shopping experience for the
customer. A recent survey [2] reports that a majority (eighty percent) of the
customers are dissatisfied with the personalization of the shopping
experience. This work illustrates that ideas from retail analytics and
sparsity-based statistical learning methods [3] can be applied to learn these
preferences. Statistical learning methods can be applied to explanatory,
descriptive, or predictive tasks [4]. The perspective in this work is
explanation and description as opposed to prediction. Insights about customer
preferences are inferred by applying interpretable models. A bipartite graph
is commonly used to represent customer shopping behavior, with customers and
store inventory representing the bipartite sets. The incidence matrix of this
graph captures the purchasing activities of the customers shopping at the
store. Matrix factorization applied to this graph yields representations of
the vertex sets of the bipartite graph in Euclidean space. These
representations are called _embeddings_. These representations permit us to
leverage a large body of machine learning and data analysis techniques that
are available for data with a vector space representation to data that has
representation as a graph. Matrix factorization can be viewed through the lens
of dictionary learning. Dictionary learning determines a sparse _basis_ for a
matrix [5] 111A variation, called the _over-complete_ dictionary, where the
opposite is true also exists. However, in this application, a sparser basis is
determined.. In the context of this application, this sparse basis yields
insights about purchase behavior and customer preferences and is the
characterization of frequent shoppers. This characterization can be
represented as a graph and stored in a _graph database_. The embeddings and
the graph characterization can be used to support ad-hoc analytic queries
about frequent shoppers as well as developing applications such as
recommendation systems. The data for this work comes from a real-world online
store. Significant preprocessing and analysis are required to develop the
representation that yields the characterization and the embeddings. The
description of the data for this work is provided in section 2. The
methodology applied to develop the characterization is discussed in section 3.
Research related to this work is discussed in section 4. The relevant
theoretical background is discussed in section 5. The results from applying
the methodology are discussed in section 6. The conclusions from this study
are discussed in 7.

##  2 Data Characterization and Problem Description

The data for this work comes from the UCI machine learning repository [6] and
represents the line items of invoices of sale from an online store in the
United Kingdom selling novelty gift items. A line-item captures the sale
details associated with a single inventory item and the identity of a customer
if the customer has registered with the store. Customers need not register to
purchase at the store. The store serves both retail and wholesale customers.
The data reflects a sample of transactions between December 2010 and December
2011 and reveals $4337$ registered customers. Of these, $447$ customers have
made five or more retail purchases at the store. Most _retail_ customers had
made less than four purchases at the store. Retail customers who make five or
more purchases at the store are determined to be _frequent shoppers_. This
group has a shopping history at the store. The characterization of this group
is developed in this work. The characterization of the other two groups,
_wholesale_ customers, and, customers who have made less than four purchases
at the store, the _infrequent shoppers_ , is not developed in this work.

Figure 1: Bipartite graph representation of frequent shoppers.

The purchase activity of the frequent shopper group can be represented as a
bipartite graph as shown in Figure 1. The _incidence matrix_ of this graph
captures the purchasing activity of the frequent shoppers. The edges of this
graph are _bi-directional_.222The edges from a customer represent his or her
shopping activity at the store. Tracing the edges from an inventory item
provides us the customers purchasing it. _The incidence matrix captures the
relationship between a customer and the store inventory. Each row of the
matrix captures the purchases of a frequent shopper. This matrix is similar to
the document-term incidence matrix encountered in information
retrieval[7][8]_. The columns of the incidence matrix are the store inventory
items purchased by the frequent shopper group. If a customer has purchased an
item, then the corresponding entry in the matrix contains the amount spent by
the customer on the item. If a customer has not purchased an item, that entry
has a $0$ value. The frequent shopper group has purchased $2664$ items. Since
there are $447$ frequent shoppers, the incidence matrix has dimensions of
$(447\times 2664)$. The entries of the incidence matrix are _non-negative_. In
this work, matrix factorization is applied to develop a concise description of
the purchase activity of frequent shoppers. This requires us to determine a
suitable _basis_ for this matrix. A preprocessing step, feature selection, can
enhance the conciseness of the results.

The notion of _Recency-Frequency-Monetary Value_(RFM) is used in retail
analytics to capture the value of a customer [6]. Recency, Frequency, and
Monetary Value are attributes associated with a shopper. The _Recency_
attribute captures the recency of a shopper’s visit to the store. The
_Frequency_ attribute captures the number of times the shopper has made
purchases at the store. The _Monetary Value_ captures the total amount of
money spent by the shopper at the store. In this work, this idea is applied to
perform feature selection. Thus we eliminate purchase activity that is
considered irrelevant. A derived attribute, _RFM score_ , denoted by
$\gamma_{i}$, is defined for the $i^{th}$ customer as in Equation 1.

$\gamma_{i}=w_{recency}Recency_{i}+w_{frequency}Frequency_{i}+w_{monetary\
value}Monetary\ Value_{i}$ (1)

where:

  * •

$w_{recency}$: weight assigned to the Recency attribute

  * •

$w_{frequency}$: weight assigned to the Frequency attribute

  * •

$w_{monetary\ value}$: weight assigned to the Monetary Value attribute

The feature selection step yields a denser and smaller representation of a
customer. _The incidence matrix after the feature selection step represents
the purchase activity of the frequent shoppers that account for customer
value_. Matrix factorization applied to this matrix yields vector
representations of the customers and store inventory items. These vector
representations called _embeddings_ can be used for several novel
applications, for example, examining the similarity of customers, searching
for customers and inventory items in a database, etc. Viewing the
factorization of the matrix as learning a sparse basis provides a useful
perspective [9] [10]. This perspective yields insights into purchase behavior
and customer preferences. This is discussed in section 5.3. The methodology
used to perform feature selection and apply matrix factorization is discussed
in section 3

##  3 Methodology

Figure 2 is a schematic diagram. This process is fairly typical in data
mining. The first step of the process is _exploratory data analysis_. As part
of this step, the characteristics of the dataset are determined. This includes
determining the types of attributes, the different types of populations in
this dataset, and the characteristics of each population. As discussed in
section 2, in this work, we are interested in the group of customers who have
made five or more purchases at the store. This group is called the _frequent
shoppers_ group. The purchasing activities of this group have a representation
as a bipartite graph. The incidence matrix of this representation is
determined. A colab [11] notebook that provides a detailed implementation of
the above steps is available at [12].

Figure 2: Methodology to develop characterization of frequent shoppers.

The incidence matrix is high dimensional. Since a typical customer purchases,
a small set of inventory items, the row capturing the customer’s purchases has
many zeros. The RFM score discussed in section 2 is used to select the
features that capture the customer’s value. A linear model that regresses the
customers’ value (RFM score) versus their purchases is used. In particular, a
LASSO model is used for feature selection. The raw RFM scores are not
distributed normally. Preprocessing the response to be normally distributed is
a common preprocessing step in the development of linear models [13] (see
section 5.1). A Box-Cox transformation [14] is used for this purpose. The
details of developing the LASSO model are provided in section 5.2. The LASSO
model provides us a smaller and denser customer representation. This
representation uses a much smaller set of inventory items to capture the
shopping activity of the frequent shoppers. This representation determines a
bipartite graph that captures most of the customer value of the frequent
shopper segment. The second step of the process applies matrix factorization
to the incidence matrix of this graph. Matrix factorization determines a
Euclidean representation of the nodes of the bipartite graph representing
customer purchase activity (after feature selection). This is achieved by
factoring the matrix into a _dictionary_ matrix and an _affinity_ matrix. The
details of the factorization and the interpretation of these matrices are
provided in section 5.3. This section discusses a method to select the hyper-
parameter associated with the Euclidean representation. Analysis of the
results from matrix factorization is discussed in section 6. The relevant
theoretical background is discussed in section 5.

##  4 Related Work

The data for this work has been analyzed by [15] by treating the _Recency_ ,
_Frequency_ and _Monetary-value_ as the primary attributes that represent the
customer. In contrast, this work represents a customer by his or her
purchases. As results from this work will show, such a representation yields
insights about customer taste and store inventory. In [15], the representation
of the customer is clustered by K-Means [16]. Specific segments that are
deemed as valuable are then profiled using a decision tree [17]. This work
also takes a rigorous data analysis perspective before customer
characterization [12]. Distinct populations of customers are defined. The
scope of this work is limited to one of the identified populations, _frequent
shoppers_. The analysis methodology is based on refining the characterization
in stages using well-established ideas in statistical learning. In [15]
clustering is used to identify groups with similar _customer value_. In this
work, clustering is used to identify groups with similar _shopping
preferences_. This work uses a density-based approach to clustering as opposed
a distance-based clustering.
Real-world dataset present with noise. This refers to points that do not
exhibit clustering tendency from a standpoint of distance. They are not close
to a neighbor. However, because we pick a finite number of clusters, such
points are assigned a cluster. As noted in [18], density-based approaches can
identify datasets with noise. Results observed with this data indicate that
noise was indeed present. In this work, feature selection based on customer
value is performed before the matrix factorization step. Exploring supervised
matrix factorization [19] [20] that can use the notion of customer value in
selecting a sparse representation of a customer is an area of future work.
Comparison of customer preferences identified in this work with other latent-
factor models, such as Latent-Dirichlet Allocation [21], is an area of future
work. Exploring graph reduction techniques, such as graph coarsening [22], to
characterize the frequent shopper group is also an area of future work.

##  5 Theoretical Background

The theoretical background about the ideas discussed in section 3 are provided
in this section. In the analysis, let $n$ be the number of customers in the
frequent shopper group and $m$ be the number of items in the store inventory
before feature selection. Results obtained using the ideas discussed in this
section are reported in section 6. The LASSO and non-negative matrix
factorization discussed in this section determine a _sparse_ basis. The LASSO
determines a _sparse_ basis that is associated with the solution of a linear
model. Matrix factorization determines a _sparse_ basis that can be used to
approximate a matrix. Seeking sparsity is a common theme in the solution of
both these problems.

###  5.1 RFM score transformation

Transforming the RFM score to be normally distributed provides several
advantages. Working with a normally distributed random variable permits us to
analyze customer value distribution using the properties of a well-known
distribution. Maximum likelihood properties based on the Normal distribution
have been extensively studied for linear models. By transforming the RFM score
to be normally distributed, we have a known framework for the evaluation of
the theoretical properties of the model [13]. The Box-Cox transformation is
used to transform the raw RFM score to a normally distributed value. In this
work, the implementation of the Box-Cox transform in scikit-learn [23] is
used. The Box-Cox transformation determines a constant, $\lambda$, that
transforms the raw RFM score of the $i^{th}$ customer, $\gamma_{i}$, to a
value $\gamma^{\prime}_{i}$ using Equation 2.

$\gamma^{\prime}_{i}=\gamma^{\lambda}_{i}=\left\\{\begin{aligned}
&\frac{(\gamma^{\lambda}_{i}-1)}{\lambda},&&\text{for}\ \lambda\neq 0\\\
&\log(\gamma_{i}),&&\text{for}\ \lambda=0\end{aligned}\right.$ (2)

where:

  * •

$\gamma^{\prime}_{i}$: is the Box-Cox transformation of the raw RFM score,
$\gamma_{i}$, of the $i^{th}$ customer

  * •

$\lambda$: is a parameter determined from the Box-Cox transformation

A linear model is used to capture the relationship between the transformed RFM
score of the $i^{th}$ customer and their purchases of store inventory as
described by Equation 3.

$\gamma^{\prime}_{i}=\beta_{0}+\beta_{1}p_{1i}+\beta_{2}p_{2i}+\ldots+\beta_{m}p_{1m}+\epsilon_{i}$
(3)

where:

  * •

$p_{i1},\ldots p_{im}$: are the purchases of store inventory items $1\ldots m$
by the $i^{th}$ customer

  * •

$\epsilon_{i}$: are independent identically distributed errors from a normal
distribution, $\mathbf{N}(0,\sigma^{2})$.

A linear relationship between the transformed RFM score and customer
purchases, as represented in Equation 3, is simple and interpretable. In
matrix form, this relationship can be written as

$\boldsymbol{\gamma^{\prime}}=\mathbf{P}.\boldsymbol{\beta}+\boldsymbol{\epsilon}$
(4)

where:

$\boldsymbol{\gamma^{\prime}}$: is a $n\times 1$ matrix of transformed
customer RFM scores

$\mathbf{P}$: is a $n\times m$ matrix of customer purchases of store inventory

$\boldsymbol{\beta}$: is a $m\times 1$ matrix of linear model coefficients

$\boldsymbol{\epsilon}$: is a $n\times 1$ matrix of errors (from
$\mathbf{N}(0,\sigma^{2})$).

However, this poses a problem. We have $447$ customers ($n=447$) and $2664$
inventory items ($m=2664$). This is an example of an _underdetermined_ system.
We have a linear system with more unknown’s ($\boldsymbol{\beta}$) than
equations (number of customers). _Regularization_ is a technique to find
solutions to an underdetermined system. The _Least Absolute Shrinkage and
Selection Operator_ (LASSO) [24] can be used for this purpose. The details of
the LASSO are provided in section 5.2. When the predictor variables are drawn
from a continuous distribution, the LASSO produces a unique solution (see
[25]).333In this case, the user’s spend on the various items of the store
inventory are the predictor variables. The notion that these are continuous
random variables is a reasonable one.

###  5.2 The Least Absolute Shrinkage and Selection Operator (LASSO)

The LASSO is a regression analysis method that can be used for regularization
and feature selection with linear models. In this work, LASSO is used to
incorporate regularization and feature selection into the linear model
described by Equation 4. The LASSO determines a sparse set of coefficients for
a linear model by minimizing the objective function shown in Equation 5. This
formulation finds a _basis_ that permits us to explain the variation in
customer value with a small number of inventory elements. This approach is
also called _basis pursuit_ [26] in the literature. A detailed exposition of
the LASSO that provides an excellent account of its formulation, solution and
properties is available in [3].

$\min_{\beta\in\mathcal{R}^{m}}\quad\|\boldsymbol{\gamma^{\prime}}-\mathbf{P}.\boldsymbol{\beta}\|+\alpha\|\boldsymbol{\beta}\|$
(5)

where:

  * •

$\boldsymbol{\gamma^{\prime}}$, $\mathbf{P}$ and $\boldsymbol{\beta}$: Are
defined as in Equation 4

  * •

$\alpha$: is a regularization hyper-parameter.

A schematic illustration of the process used to perform feature selection
using the LASSO model is shown in Figure 3. The details of each step of the
process illustrated in Figure 3 are as follows.

Figure 3: Using the LASSO for feature selection.

  1. 1.

Select Regularization Parameter, $\alpha$: Applying the LASSO model for
feature selection requires the selection of the regularization parameter,
$\alpha$. This is selected by considering a range of values of $\alpha$ and
picking the value that produces the lowest cross-validation error [27].

  2. 2.

Select the number of features: Once the $\alpha$ parameter has been selected,
the LASSO model can be computed and used to perform feature selection. The
LASSO solution drives the coefficients in the linear model, $\beta$, to zero.
The absolute value of the coefficient is a direct measure of its ability to
explain customer value. An ordering of these coefficients in descending order
provides a measure of feature importance of an item of inventory in explaining
the customer value of the frequent shopper group. Many of these coefficients
are very small and do not explain customer value. We only want those
coefficients that can account for explaining a meaningful amount of customer
value. To determine the number of coefficients to select from the model, the
following experimental procedure is used. A portion of the dataset is held out
as a test set. The solution from the lasso is initially selected as the set of
features to include in the model. A linear model is developed with these
features and its performance on the test set is noted. The smallest
coefficient is dropped and a model is developed with the trimmed feature set.
The performance of the model on the test set is noted. This process of
dropping the feature associated with the smallest coefficient is repeated and
model performance is noted. We then plot the number of features included in
the model versus test set performance. This plot reveals the optimal number of
features to include in the model. The experimental procedure is described in
Algorithm 1. _This procedure yields a set of $m^{\prime}$ coefficients that
can explain most of the value in the frequent shopper group. Results showed
that $m^{\prime}\ll m$. In other words, a small set of store inventory items
account for most of the customer value_. The details are discussed in section
6, the relevant fact is that such an experiment reveals that there diminishing
returns to including a feature in terms of its ability to explain customer
value. The above procedure gives us an experimental estimate of the number of
store inventory items that contribute towards explaining customer value.

Input: The LASSO solution, $\beta_{\text{LASSO}}$, The hold-out test set

Data: The frequent customer training data

plot data = [ ] /* Compute hold-out test error with LASSO solution */

num_features = number of features in the LASSO solution

fit_linear_model(num_features)

add_holout_error_to_plot_data()

/* Observe the hold-out error as we drop features from the feature set for the
linear model */

while  _num_features $>0$_ do

fit_linear_model(num_features)

add_holout_error_to_plot_data()

/* Drop one feature */

drop_feature_with_smallest_coefficient()

plot_experiment_curve(plot data)

Algorithm 1 An experiment to select number of features from the LASSO solution

  3. 3.

Validate that the assumption of a linear model: We still need to assess if the
proposed linear model for estimating customer value is a reasonable choice.
This assessment is done by:

    1. (a)

Plotting the model estimates of customer value versus the actual customer
value for the test set.

    2. (b)

Examining a probability-probability plot of the errors from the test set
against a normal distribution to evaluate if normality of errors is a
reasonable assumption.

A colab [11] notebook that provides a detailed implementation of feature
selection using the LASSO as illustrated in Figure 3 is available [28].

###  5.3 Non-Negative Matrix Factorization

Feature selection performed using the LASSO provides a small number of
features, $m^{\prime}$ that can account for most of the customer value. The
matrix representing customer purchases after feature selection,
$\mathbf{P}^{\prime}$, has fewer columns ($m^{\prime}$) than the original
matrix of purchases, $\mathbf{P}$, which has $m$ columns. _The matrix
$\mathbf{P}^{\prime}$ is the incidence matrix of a bipartite graph that
captures the purchasing activity of the frequent shopper group_. It should be
noted that the entries $p_{ij}$ of $\mathbf{P}^{\prime}$ represent the
purchase of item $j$ by customer $i$. This is a strictly non-negative
quantity. Non-negative matrix factorization can be used to factorize the
matrix $\mathbf{P}^{\prime}$ into factors $\mathbf{W}$ and $\mathbf{H}$ such
that $\mathbf{P}^{\prime}\approx\mathbf{W}.\mathbf{H}$, as shown in Figure 4.

Figure 4: Matrix factorization factors, number of customers ($n$), number of
inventory items ($m$), latent dimension $k$.

The matrix $\mathbf{W}$ of dimension $n\times k$ is tall and narrow, the
matrix $\mathbf{H}$ of dimension $k\times m$, is shallow and wide. Matrix
factorization approximates the purchase of the $i^{th}$ customer as the
product of two factors, a dictionary element $h_{km}$ and an affinity element
$w_{ik}$, such that the purchases of the customer can be expressed as shown in
Equation 6.

$p_{i}\approx w_{ik}.h_{km}$ (6)

The dimension $k$ is a hyper-parameter to the method and is determined through
an imputation experiment. A third of the entries of the matrix
$\mathbf{P}^{\prime}$ are held out as a test set. A range of values, from $2$
to $20$, are used as candidate values for $k$. For each value of $k$, the mean
square error (MSE) in predicting the entries of the test set is computed. The
value of $k$ that produces the lowest mean square error on the test set is
selected as the optimal value of the latent dimension. We now present two
perspectives of matrix factorization to motivate the essential aspects of the
developing the solution.

  1. 1.

A dictionary learning perspective: The matrix $\mathbf{H}$ represents a
dictionary of purchases. It represents a _basis_ for a latent space of
dimension $k$. The purchases of any customer $i$ can be expressed as a linear
combination of the elements of the dictionary, as expressed by Equation 6.
Consequently, matrix factorization yields a dictionary of purchases that we
can group the store inventory into and a matrix of customer affinities towards
these dictionary elements. The dictionary elements and the affinities have
vector space representations. The rows of the matrix $\mathbf{W}$, the
affinity matrix, represent a $k$ dimensional representation of the customer.
The columns of the dictionary matrix $\mathbf{H}$ represent the $k$
dimensional representation of the store inventory items. These representations
can be used by applications such as similarity search[29]. By examining the
rows of the matrix $\mathbf{H}$ the inventory items that dominate each of the
$k$ dictionary items can be determined. An entry $h_{ij}$ of a dictionary
element can be viewed as the contribution of the $j^{th}$ inventory item to
the $i^{th}$ dictionary element. A high value implies that the $j{th}$
inventory item contributes strongly towards the $i^{th}$ dictionary item. By
examining the entries of $\mathbf{H}$ with a row-wise scan, we can identify
the inventory items that dominate each dictionary element. Similary, an entry
$w_{ij}$ can be interpretted as the affinity that customer $i$ has towards
dictionary element $j$. By examining the $i^{th}$ row of $\mathbf{W}$, we can
identify the preferences of the $i^{th}$ customer towards the dictionary
elements. 444There elements of the dictionary are also called _atoms_[30] and
_parts_ [31] in the literature. Similarly, the affinity component is called
_activations_[32] in the literature.

The problem of determining a non-negative matrix factorization is _ill-posed_.
There can multiple factorizations for a given data matrix
$\mathbf{P}^{\prime}$. There have been several approaches to mitigate the non-
uniqueness of the factorization by imposing different kinds of constraints,
for example, using convexity [33] or sparsity [34]. The perspective in this
work is explanation and description. Occam’s razor [35][36] is applied in this
work. We seek simple or sparse explanations. The sparsity constraint is
imposed in the optimization model used to determine the factorization. This is
discussed next.

  2. 2.

An optimization perspective: Non-negative matrix factorization applied to the
matrix $\mathbf{P}^{\prime}$ requires us to determine matrices $\mathbf{H}$
and $\mathbf{W}$ such that $\mathbf{P}^{\prime}\approx\mathbf{W}.\mathbf{H}$.
These matrices are determined by framing matrix factorization as an
optimization problem. In this work the implementation of matrix factorization
in scikit-learn[23] is used. This implementation uses a flexible cost
function, called _$\beta$ divergence_[34], to compute the matrix
factorization. This cost function is given by Equation 7. Regularization using
a mix of both $\ell^{1}and\ell^{2}$ penalties is possible with this
implementation by supplying the algorithm with a ratio parameter,
$\ell^{1}_{ratio}$. The sparsity constraint is imposed by the choice of the
parameter $\alpha_{m}$. In this work, the $\ell^{1}_{ratio}$ parameter and the
$\alpha_{m}$ parameter are selected by a grid search. A range of values for
the $\ell^{1}_{ratio}$ and $\alpha_{m}$ parameters are considered as candidate
choices. The set of values that minimized the imputation error in the
experiment described above was used for developing the matrix factorization
model.

$\displaystyle\operatorname*{arg\,min}_{\mathbf{W}\in\mathcal{R}^{n\times k},\
\mathbf{H}\in\mathcal{R}^{k\times m}}$
$\displaystyle\Bigg{\\{}\frac{1}{2}\|\mathbf{P}^{\prime}-\mathbf{W}.\mathbf{H}\|+\alpha_{m}.\ell^{1}_{ratio}.\|\mathbf{W}\|_{\ell^{1}}+$
(7)
$\displaystyle\alpha_{m}.\ell^{1}_{ratio}.\|\mathbf{H}\|_{\ell^{1}}+\frac{1}{2}.\alpha_{m}.(1-\ell^{1}_{ratio})\|\mathbf{W}\|^{2}_{Frob}+$
$\displaystyle\frac{1}{2}.\alpha_{m}.(1-\ell^{1}_{ratio})\|\mathbf{H}\|^{2}_{Frob}\Bigg{\\}}$

where:

     * •

$\mathbf{P}^{\prime}$: is a matrix of purchases after feature selection, of
size $n\times m^{\prime}$

     * •

$\mathbf{H}$: is a dictionary matrix of purchases, of size $m^{\prime}\times
k$

     * •

$\mathbf{W}$: is a matrix of customer affinities, of size $n\times k$

     * •

$\alpha_{m}$: is a constant that multiplies the regularization terms.

     * •

$\ell^{1}_{ratio}$: a mixing parameter that specifies the ratio of
$\ell^{1}and\ell^{2}$ penalties to use with the optimization.

     * •

$\|\mathbf{W}\|_{\ell^{1}},\|\mathbf{H}\|_{\ell^{1}}$: are element-wise
$\ell^{1}$ norms, $\sum_{i,j}\|w_{ij}\|$ and $\sum_{i,j}\|h_{ij}\|$ of
$\mathbf{W}$ and $\mathbf{H}$.

     * •

$\|\mathbf{W}\|^{2}_{Frob},\|\mathbf{H}\|^{2}_{Frob}$: are Frobenius norms,
$\sum_{i,j}w_{ij}^{2}$ and $\sum_{i,j}h_{ij}^{2}$ of $\mathbf{W}$ and
$\mathbf{H}$.

A flow chart illustrating the process of applying matrix factorization to
develop the embeddings of the frequent shopper graph and develop a profile of
customer preferences is shown in Figure 5.

Figure 5: Developing embeddings and profiling customer preferences with matrix
factorization.

##  6 Results

This section presents the key results from applying the methodology discussed
in section 3. The _exploratory data analysis_ step identifies the frequent
shopper group. The details of _exploratory data analysis_ are available at
[12]. The purchasing activities of the frequent shopper group have a
representation as a bipartite graph. The incidence matrix of this graph is
used to develop the characterization of the frequent shopper group. The
incidence matrix developed from the raw data set provides a high-dimensional
representation of the customer. A typical _retail_ customer purchases a small
number of inventory items. Therefore, most entries associated with the
customer in the incidence matrix are zeros. Feature selection using the LASSO
[24] is used to develop a denser representation of the customer as discussed
in section 5.2. The results from this step are discussed in section 6.1.
Matrix factorization provides a vector representation of the nodes of the
frequent shopper graph, i.e, of customers and the store inventory. The results
from the matrix factorization step are provided in section 6.2.

###  6.1 Feature Selection with the LASSO

A flow chart illustrating the steps involved in using the LASSO to perform
feature selection is shown in Figure 3. A detailed implementation of this
flowchart as a colab notebook is available at [28]. This section captures the
salient results from executing the process depicted in Figure 3. The
theoretical background is discussed in section 5.2. The parameters used for
the process depicted in Figure 3 and results from the feature selection
process are summarized in Table 1.

Table 1: Parameters and results from feature selection. FS Parameter/Result | Value
---|---
RFM weights | $w_{r}=0.15,$
| $w_{f}=0.15,$
| $w_{m}=0.7$
Regularization parameter, $\alpha$ | 0.22
$R^{2}$ (explained variance) | 0.855
optimal number of features | 75
training MSE | 0.137
test MSE | 0.326

The LASSO drives the coefficients of the coefficients, $\beta$, in Equation 5
to zero. Many of these coefficients are very small and do not contribute in a
meaningful manner to explaining variance in customer value
($\gamma^{\prime}$). The experiment described by Algorithm 1 is used to select
those features that contribute to explaining the variation in customer value.
The results from this experiment are shown in Figure 6.

Figure 6: Optimal number of features to select from LASSO solution.

Since a linear model is used, the magnitude of the coefficients, $\beta$, is a
direct indication of an inventory item in determining customer value. An
inspection of Figure 6 shows that the first $50$ features account for
explaining most customer value. There are diminishing returns in including the
inventory items associated with the coefficients after the $50^{th}$
coefficient in the representation of the customer. About $75$ inventory items
can explain about $86\%$ of the variation in customer value. Therefore,
purchases of these $75$ store items are used to characterize customer behavior
in subsequent analysis. The feature importance of the first $50$ features is
shown in Figure 7.

Figure 7: Feature importance of the first 50 features.

A description of the top ten items in terms of feature importance are shown in
Table 2. _It should be noted that the features in Table 2 are features that
are determined important based on value (RFM score). These are the features
that account for customer value in contrast to customer taste. The dictionary
learned from matrix factorization, discussed in the next section, provides
insights about customer taste._

Table 2: Description of the top 10 items in feature importance. ITEM ID | Description
---|---
20617 | FIRST CLASS PASSPORT COVER
10002 | INFLATABLE POLITICAL GLOBE
10133 | COLOURING PENCILS BROWN TUBE
10135 | COLOURING PENCILS BROWN TUBE
17012D | ORIGAMI ROSE INCENSE/CANDLE SET
10120 | DOGGY RUBBER
15039 | SANDALWOOD FAN
17084N | FAIRY DREAMS INCENSE
15056P | EDWARDIAN PARASOL PINK
10080 | GROOVY CACTUS INFLATABLE

The details of validating that the linear model assumption for customer value
is a reasonable one are provided in section 5.2. The results of executing
these steps are available in [28]. Results showed that a linear model captures
the relationship between the customer value and their purchases accurately. A
comparison of the error distribution from the model to a normal distribution
showed that it is reasonable to assume that the errors from the model
estimates follow a normal distribution.

###  6.2 Matrix Factorization

The flow chart shown in Figure 5 illustrates the steps involved in applying
matrix factorization to the purchase data after feature selection. The first
step in the process is to determine an _optimal_ value of the dimensionality,
$k$, the hyperparameters, $\ell^{1}_{\text{ratio}}$ and $\alpha_{m}$
associated with computing the matrix factorization of the matrix
$\mathbf{P}{\prime}$. A range of values for these parameters were considered
(see Figure 8 for the range of values) and the values that minimized the
imputation error on the held out test set was used to develop the matrix
factorization model. The results from this experiment are shown in Figure 8.
The best set of parameters from the experiment are shown in Table 3. The
implementation of the experiment is available as a colab[11] notebook[37].

Figure 8: Matrix factorization, choosing the best k. Table 3: Matrix Factorization Parameters NMF parameter | Value
---|---
$\ell^{1}_{\text{ratio}}$ | $0.1$
$\alpha_{m}$ | $1.0$
$k$ | $5$

###  6.3 Profiling the Matrix Factorization Solution

Matrix factorization is applied to the incidence matrix with the latent
dimension size, $k$, set to $5$ based on the results from the experiment
discussed in the previous section. Matrix factorization provides matrices
$\mathbf{W}$ and $\mathbf{H}$. The interpretation of these matrices is
discussed in section 5. The matrix $\mathbf{H}$ is called the dictionary. The
dictionary has $5$ elements. The purchases of the $i^{th}$ customer, $p_{i}$,
can be expressed in terms of his/her affinities towards the dictionary
elements, as indicated by Equation 6. This is illustrated in Figure 9 where
each customer’s affinities towards the dictionary elements can be represented
by a pyramid.

Figure 9: Interpretation of Matrix Factorization.

To determine the relative contribution of each inventory item to a dictionary
element, we can normalize the rows of the matrix $\mathbf{H}$ to be of unit
length. This provides a quick view of the relative importance of the inventory
items in a dictionary element. The results of this computation are shown in
Figure 10

Figure 10: Profile of the frequent shopper dictionary.

A review of Figure 10 shows that each dictionary element is dominated by a
small number of inventory items. Table 4 provides the description of the top 3
inventory for each of the five dictionary elements. The dictionary elements
can overlap555Formulations of non-negative matrix factorization that yield
orthogonal basis are also possible. Note for example that the _christmas
bunting_ figures both the first and second dictionary elements in Table 4. The
reader may be tempted to compare Table 2 with Table 4. It should be noted that
dictionary learning is based on the objective function in Equation 7 and is
_unsupervised_ learning. In contrast, the feature importance shown in Table 2
is based on _supervised_ learning with the objective function described in
Equation 5. These facts should be kept in mind when interpreting the results.
The _customer-value_ is not explicitly used in learning the dictionary
representation. Exploring this is an area of future work, as discussed in
section 4. It is possible to select a specific basis (rather than learn it)
when developing the matrix factorization solution. For example, if we like the
set of dictionary elements learned in this work, we can apply it in future
analysis to develop a solution with a new sample of data, see [38] for
details. In other words, incorporating prior knowledge about the dictionary is
an option with non-negative matrix factorization. To the best of our
knowledge, we did not find other work analyzing this data.

Table 4: Profile of the dictionary of purchases. Item_group | Item_ID | Description
---|---|---
Item Group_0 | 11 | FELTCRAFT PRINCESS CHARLOTTE DOLL
Item Group_0 | 26 | 6 RIBBONS RUSTIC CHARM
Item Group_0 | 27 | JUMBO BAG PINK POLKADOT
Item Group_1 | 30 | VINTAGE CHRISTMAS BUNTING
Item Group_1 | 48 | 3 DRAWER ANTIQUE WHITE WOOD CABINET
Item Group_1 | 44 | JUMBO BAG VINTAGE LEAF
Item Group_2 | 37 | BLUE DINER WALL CLOCK
Item Group_2 | 30 | VINTAGE CHRISTMAS BUNTING
Item Group_2 | 41 | LUNCH BAG VINTAGE LEAF DESIGN
Item Group_3 | 74 | SCANDINAVIAN REDS RIBBONS
Item Group_3 | 57 | DOORMAT 3 SMILEY CATS
Item Group_3 | 54 | BOX OF VINTAGE JIGSAW BLOCKS
Item Group_4 | 0 | SPACEBOY GIFT WRAP
Item Group_4 | 39 | JUMBO BAG 50’S CHRISTMAS
Item Group_4 | 23 | POPCORN HOLDER

The rows of the matrix, $\mathbf{W}$, represent the affinities of the
customers towards the dictionary elements. By clustering the matrix
$\mathbf{W}$, we can see how preferences of the store customers towards the
dictionary elements group together. The hdbscan[39][40] algorithm is used for
this purpose. hdbscan can discriminate between _noise_ and clustering
tendency. The algorithm requires us to specify the minimum number of points to
consider aggregating into a cluster. Given the size of the frequent shopper
group ($447$ customers), this is set to $5$. hdbscan produces $5$ clusters in
the data with a large noise cluster ($338$). One of these clusters, Cluster
$3$, is large, with $77$ customers. There is a cluster with $14$ elements -
Cluster $2$. There are three small clusters, Cluster $0$, Cluster $4$, and
Cluster $1$, each with $6$ customers. A size distribution of the clusters is
shown in Figure 11.

Figure 11: Cluster sizes from hdbscan.

A review of the results from hdbscan reveals there is a large group of
customers ($338$) that we cannot assign a cluster to. This is the noise
cluster. The purchase activity of this group does not show clustering
tendencies. There is a smaller group of $109$ 666The size of the smaller group
with clustering tendencies is simply the difference between the size of the
group and the size of the noise cluster. customers whose purchasing activities
show clear clustering tendencies. The cluster centroids of each of these
clusters can be computed. Each centroid can then be normalized to unit length.
An inspection of the coordinates of the normalized centroids reveals the
affinities each cluster has to the dictionary elements. The results from such
a computation are shown in Figure 12. Cluster $0$ has a clear preference for
the dictionary element tagged Item Group_1. Clusters $1$ and Cluster $2$ have
a strong affinity to the dictionary element tagged Item Group_0. These
clusters are over-lapping clusters that are characterized by a strong
preference towards dictionary element Item Group_0 and slightly different and
small affinities towards the other dictionary elements. Cluster $3$, the large
cluster, has affinities for dictionary elements Item Group_1, Item Group_3 and
Item Group_4. There is a strong preference for Item Group_4 in Cluster $3$.
Cluster $4$ is characterized by affinity for dictionary element Item Group_4.
Cluster $4$ does not show strong affinities for Item Group_1 and Item Group_3.
Overall, it is clear that the dictionary element, Item Group_4, is favored by
a majority of customers with clear shopping preferences. The elements of this
dictionary element, see Table 4, capture an important facet of the taste of
the frequent shoppers. Ensuring the supply and availability of the inventory
items that dominate dictionary element, Item Group_4, is important to keep
store customers satisfied.

Figure 12: Cluster centroid affinities.

A graph is an excellent data model for this problem. Both the business problem
and the results from its analysis have very elegant interpretations as graphs.
The customer’s purchase activity could be modeled and interpreted using a
bipartite graph. Similarly, the results from this analysis also have a natural
interpretation as a bipartite graph. The bipartite graph capturing the
purchase activity related to the small set of inventory items capturing
customer value is one of the results from this analysis. The results from
matrix factorization also have an interpretation as a bipartite graph with the
customers and the dictionary elements forming the vertex sets. The matrix
$\mathbf{W}$ represents the edge weights between the customers and the
dictionary vertex sets. These graphs can be stored in a _Graph Database_. For
the storage, we used an ArangoDB graph database. The embeddings of the vertex
sets can also be stored. This graph representation can be used to capture and
query the characteristics of frequent shoppers to the store. For example, an
analyst can query this graph to see if a customer is similar to another
customer. The data for this work covers a particular period. Similar analyses
performed for other periods produce similar results. An analyst can then query
these graphs to see understand how customer preferences change over time. For
example, an analyst could examine how dictionary profiles and customer cluster
profiles change over time (see [37] for examples of such queries).

###  6.4  Discussion of Results

The feature selection models illustrate that a small set of $75$ store
inventory items (out of the $2664$ items in the inventory) account for about
$86\%$ of the variation in customer value. Since a linear model was used, the
elements of this set can be ordered based on their associated coefficients in
the model. The first $50$ of these elements account for most of the customer
value explained by the linear model. Managing the inventory and restocking of
these elements is critical from the standpoint of store operations. Matrix
factorization groups these $75$ inventory items into dictionary elements.
Analysis of the grouping showed that each dictionary element is characterized
by a small set of inventory items. The dictionary elements are not orthogonal.
This implies that there can overlap between the items that characterize the
dictionary elements. Formulations of non-negative matrix factorization that
yield orthogonal bases are possible, for example, [41]. A customer’s purchase
can be viewed in terms of affinity towards these dictionary elements (Equation
6). Clustering of customer affinities revealed that there is a rather large
group whose tastes do not align with other customers. There is a small group
of $109$ customers whose taste aligns with other customers. A profile of the
cluster centroids provided the preferences of each of these groups to the
dictionary elements. The characterization of frequent shoppers developed in
this work can be stored and analyzed as a graph in _Graph database_ such as
ArangoDB. The _embeddings_ of the vertices of the graph are also stored as
properties of graph nodes. These embeddings can be used to develop
applications such as recommendation systems [42]. In summary, the methodology
applied in this work provides embeddings, insights, and interpretable models!

##  7 Conclusion

This work uses sparsity-seeking models to develop the characterization of the
frequent shoppers to a real-world online store. The models used were
interpretable and yielded insights about frequent customers to the store. The
theoretical properties of each of the methods applied are well established.
Methods to apply these techniques to large datasets are also established. The
perspective in this work is analysis as opposed to forecasting or prediction.
To the best of our knowledge, there is no other work analyzing the group of
frequent shoppers associated with this data from a shopping preferences
perspective. The use of interpretable models permits us to identify a set of
facts that form a baseline for future work on this data. Exploring supervised
matrix factorization is an area of future work. A comparison of the results
from this work to results from graph reduction approaches, such as coarsening,
would be interesting. A comparison of the analytical techniques from the
perspective of information provided, cost and computational characteristics
etc. is an area for future work. Characterization of the wholesale shoppers
and infrequent shoppers is also an area of future work.

## References

  * [1] Davenport, Thomas and Dallamule Leonardo and Lucker John.  Know What Your Customers Want Before They Do, 2011.
  * [2] Phibbs, Bob.  What do shoppers really want? Do retailers have a clue?, 2021.
  * [3] Trevor Hastie, Robert Tibshirani, and Martin Wainwright.  Statistical learning with sparsity: the lasso and generalizations.  Chapman and Hall/CRC, 2019.
  * [4] Galit Shmueli.  To explain or to predict?  Statistical science, 25(3):289–310, 2010.
  * [5] Ankur Moitra.  Algorithmic aspects of machine learning.  Cambridge University Press, 2018.
  * [6] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.
  * [7] Manning, Christopher and Raghavan, Prabhakar and Schütze, Hinrich.  An example information retrieval problem, 2008.
  * [8] Christopher D. Manning, Prabhakar Raghavan, and Hinrich Schütze.  Introduction to information retrieval.  Cambridge University Press, 2008.
  * [9] Julian Eggert and Edgar Korner.  Sparse coding and nmf.  In 2004 IEEE International Joint Conference on Neural Networks (IEEE Cat. No. 04CH37541), volume 4, pages 2529–2533. IEEE, 2004.
  * [10] Nicolai Meinshausen.  Multivariate Statistics – non-negative matrix factorisation and sparse dictionary learning, 2021.
  * [11] Google.  Google Colab, 2021.
  * [12] ArangoDB.  Exploratory Data Analysis of a Retail Dataset, 2021.
  * [13] John Fox.  Regression diagnostics: An introduction.  Sage publications, 2019.
  * [14] George EP Box and David R Cox.  An analysis of transformations.  Journal of the Royal Statistical Society: Series B (Methodological), 26(2):211–243, 1964.
  * [15] Daqing Chen, Sai Laing Sain, and Kun Guo.  Data mining for the online retail industry: A case study of rfm model-based customer segmentation using data mining.  Journal of Database Marketing & Customer Strategy Management, 19(3):197–208, 2012.
  * [16] Stuart Lloyd.  Least squares quantization in pcm.  IEEE transactions on information theory, 28(2):129–137, 1982.
  * [17] Leo Breiman, Jerome H Friedman, Richard A Olshen, and Charles J Stone.  Classification and regression trees.  Routledge, 2017.
  * [18] Leland McInnes and John Healy.  Accelerated hierarchical density based clustering.  In 2017 IEEE International Conference on Data Mining Workshops (ICDMW), pages 33–42. IEEE, 2017.
  * [19] Woody Austin, Dylan Anderson, and Joydeep Ghosh.  Fully supervised non-negative matrix factorization for feature extraction.  In IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium, pages 5772–5775. IEEE, 2018.
  * [20] Jaegul Choo, Changhyun Lee, Chandan K Reddy, and Haesun Park.  Weakly supervised nonnegative matrix factorization for user-driven clustering.  Data mining and knowledge discovery, 29(6):1598–1621, 2015.
  * [21] David M Blei, Andrew Y Ng, and Michael I Jordan.  Latent dirichlet allocation.  the Journal of machine Learning research, 3:993–1022, 2003.
  * [22] Andreas Loukas.  Graph reduction with spectral and cut guarantees.  Journal of Machine Learning Research, 20(116):1–42, 2019.
  * [23] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [24] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 1996.
  * [25] Ryan J Tibshirani.  The lasso problem and uniqueness.  Electronic Journal of statistics, 7:1456–1490, 2013.
  * [26] Scott Shaobing Chen and David L Donoho.  Application of basis pursuit in spectrum estimation.  In Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP’98 (Cat. No. 98CH36181), volume 3, pages 1865–1868. IEEE, 1998.
  * [27] Tibshirani, Ryan.  K-Fold Cross Validation, 2021.
  * [28] ArangoDB.  Feature Selection Using the LASSO on a Retail Dataset, 2021.
  * [29] Jeff Johnson, Matthijs Douze, and Hervé Jégou.  Billion-scale similarity search with gpus.  IEEE Transactions on Big Data, 2019.
  * [30] Ronan Hamon, Valentin Emiya, and Cédric Févotte.  Convex nonnegative matrix factorization with missing data.  In 2016 IEEE 26th International Workshop on Machine Learning for Signal Processing (MLSP), pages 1–6. IEEE, 2016.
  * [31] Daniel D Lee and H Sebastian Seung.  Learning the parts of objects by non-negative matrix factorization.  Nature, 401(6755):788–791, 1999.
  * [32] Slim Essid and Alexey Ozerov.  A tutorial on nonnegative matrix factorisation with applications to audiovisual content analysis.  In Tutorial at ICME 2014, 2014.
  * [33] Chris H. Q. Ding, Tao Li, and Michael I. Jordan.  Convex and semi-nonnegative matrix factorizations.  IEEE Trans. Pattern Anal. Mach. Intell., 32(1):45–55, 2010.
  * [34] Cédric Févotte and Jérôme Idier.  Algorithms for nonnegative matrix factorization with the $\beta$-divergence.  Neural computation, 23(9):2421–2456, 2011.
  * [35] Pedro M. Domingos.  The role of occam’s razor in knowledge discovery.  Data Min. Knowl. Discov., 3(4):409–425, 1999.
  * [36] Pedro M. Domingos.  A few useful things to know about machine learning.  Commun. ACM, 55(10):78–87, 2012.
  * [37] ArangoDB.  Matrix Factorization of Frequent Shopper Purchases, 2021.
  * [38] Krithika Manohar, Bingni W Brunton, J Nathan Kutz, and Steven L Brunton.  Data-driven sparse sensor placement for reconstruction: Demonstrating the benefits of exploiting known patterns.  IEEE Control Systems Magazine, 38(3):63–86, 2018.
  * [39] Leland McInnes, John Healy, and Steve Astels.  hdbscan: Hierarchical density based clustering.  Journal of Open Source Software, 2(11):205, 2017.
  * [40] Ricardo JGB Campello, Davoud Moulavi, and Jörg Sander.  Density-based clustering based on hierarchical density estimates.  In Pacific-Asia conference on knowledge discovery and data mining, pages 160–172. Springer, 2013.
  * [41] Megasthenis Asteris, Dimitris Papailiopoulos, and Alexandros G Dimakis.  Orthogonal nmf through subspace exploration.  Advances in Neural Information Processing Systems, 28:343–351, 2015\.
  * [42] Asnat Greenstein-Messica, Lior Rokach, and Michael Friedman.  Session-based recommendations using item embedding.  In Proceedings of the 22nd International Conference on Intelligent User Interfaces, IUI ’17, page 629–633, New York, NY, USA, 2017\. Association for Computing Machinery.


* Kewley et al. (2019) Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511. doi:10.1146/annurev-astro-081817-051832
  * Kewley et al. (2013) Kewley, L. J., Maier, C., Yabe, K., et al. 2013, ApJ, 774, L10. doi:10.1088/2041-8205/774/1/L10
  * Kodra et al. (2023) Kodra, D., Andrews, B. H., Newman, J. A., et al. 2023, ApJ, 942, 36. doi:10.3847/1538-4357/ac9f12
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36. doi:10.1088/0067-0049/197/2/36
  * Kravchenko et al. (2022) Kravchenko, K., Dallilar, Y., Absil, O., et al. 2022, Proc. SPIE, 12184, 121845M. doi:10.1117/12.2629258
  * LaMassa et al. (2019) LaMassa, S. M., Georgakakis, A., Vivek, M., et al. 2019, ApJ, 876, 50. doi:10.3847/1538-4357/ab108b
  * Lamperti et al. (2017) Lamperti, I., Koss, M., Trakhtenbrot, B., et al. 2017, MNRAS, 467, 540. doi:10.1093/mnras/stx055
  * Le Fèvre et al. (2003) Le Fèvre, O., Saisse, M., Mancini, D., et al. 2003, Proc. SPIE, 4841, 1670. doi:10.1117/12.460959
  * Lester et al. (1990) Lester, D. F., Carr, J. S., Joy, M., et al. 1990, ApJ, 352, 544. doi:10.1086/168557
  * Lilly et al. (2007) Lilly, S. J., Le Fèvre, O., Renzini, A., et al. 2007, ApJS, 172, 70. doi:10.1086/516589
  * Maiolino et al. (2020) Maiolino, R., Cirasuolo, M., Afonso, J., et al. 2020, The Messenger, 180, 24. doi:10.18727/0722-6691/5197
  * Markwardt (2009) Markwardt, C. B. 2009, Astronomical Data Analysis Software and Systems XVIII, 411, 251. doi:10.48550/arXiv.0902.2850
  * Masters et al. (2016) Masters, D., Faisst, A., & Capak, P. 2016, ApJ, 828, 18. doi:10.3847/0004-637X/828/1/18
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425. doi:10.1086/155591
  * Matsuoka et al. (2009) Matsuoka, K., Nagao, T., Maiolino, R., et al. 2009, A&A, 503, 721. doi:10.1051/0004-6361/200811478
  * Mignoli et al. (2013) Mignoli, M., Vignali, C., Gilli, R., et al. 2013, A&A, 556, A29. doi:10.1051/0004-6361/201220846
  * Mignoli et al. (2019) Mignoli, M., Feltre, A., Bongiorno, A., et al. 2019, A&A, 626, A9. doi:10.1051/0004-6361/201935062
  * Morisset et al. (2015) Morisset, C., Delgado-Inglada, G., & Flores-Fajardo, N. 2015, Rev. Mexicana Astron. Astrofis., 51, 103. doi:10.48550/arXiv.1412.5349
  * Mouri et al. (1990) Mouri, H., Nishida, M., Taniguchi, Y., et al. 1990, ApJ, 360, 55. doi:10.1086/169095
  * Nagao et al. (2006) Nagao, T., Maiolino, R., & Marconi, A. 2006, A&A, 447, 863. doi:10.1051/0004-6361:20054127
  * Nandra et al. (2015) Nandra, K., Laird, E. S., Aird, J. A., et al. 2015, ApJS, 220, 10. doi:10.1088/0067-0049/220/1/10
  * Noeske et al. (2007) Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43. doi:10.1086/517926
  * Oesch et al. (2023) Oesch, P. A., Brammer, G., Naidu, R. P., et al. 2023, arXiv:2304.02026. doi:10.48550/arXiv.2304.02026
  * Oliva et al. (2001) Oliva, E., Marconi, A., Maiolino, R., et al. 2001, A&A, 369, L5. doi:10.1051/0004-6361:20010214
  * Onori et al. (2017) Onori, F., La Franca, F., Ricci, F., et al. 2017, MNRAS, 464, 1783. doi:10.1093/mnras/stw2368
  * Osterbrock (1989) Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei, by Donald E. Osterbrock. Published by University Science Books, ISBN 0-935702-22-9, 408pp, 1989.
  * Pannella et al. (2015) Pannella, M., Elbaz, D., Daddi, E., et al. 2015, ApJ, 807, 141. doi:10.1088/0004-637X/807/2/141
  * Rayner et al. (2003) Rayner, J. T., Toomey, D. W., Onaka, P. M., et al. 2003, PASP, 115, 362. doi:10.1086/367745
  * Ricci et al. (2022) Ricci, C., Ananna, T. T., Temple, M. J., et al. 2022, ApJ, 938, 67. doi:10.3847/1538-4357/ac8e67
  * Richardson et al. (2014) Richardson, C. T., Allen, J. T., Baldwin, J. A., et al. 2014, MNRAS, 437, 2376. doi:10.1093/mnras/stt2056
  * Riffel et al. (2006) Riffel, R., Rodríguez-Ardila, A., & Pastoriza, M. G. 2006, A&A, 457, 61. doi:10.1051/0004-6361:20065291
  * Riffel et al. (2013) Riffel, R., Rodríguez-Ardila, A., Aleman, I., et al. 2013, MNRAS, 430, 2002. doi:10.1093/mnras/stt026
  * Riffel et al. (2019) Riffel, R., Rodríguez-Ardila, A., Brotherton, M. S., et al. 2019, MNRAS, 486, 3228. doi:10.1093/mnras/stz1077
  * Risaliti et al. (2000) Risaliti, G., Maiolino, R., & Bassani, L. 2000, A&A, 356, 33. doi:10.48550/arXiv.astro-ph/0002169
  * Rodríguez-Ardila et al. (2004) Rodríguez-Ardila, A., Pastoriza, M. G., Viegas, S., et al. 2004, A&A, 425, 457. doi:10.1051/0004-6361:20034285
  * Rodríguez-Ardila et al. (2005) Rodríguez-Ardila, A., Riffel, R., & Pastoriza, M. G. 2005, MNRAS, 364, 1041. doi:10.1111/j.1365-2966.2005.09638.x
  * Sanders et al. (2020) Sanders, R. L., Shapley, A. E., Reddy, N. A., et al. 2020, MNRAS, 491, 1427. doi:10.1093/mnras/stz3032
  * Sargent et al. (2014) Sargent, M. T., Daddi, E., Béthermin, M., et al. 2014, ApJ, 793, 19. doi:10.1088/0004-637X/793/1/19
  * Savage & Sembach (1996) Savage, B. D. & Sembach, K. R. 1996, ApJ, 470, 893. doi:10.1086/177919
  * Schreiber et al. (2015) Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74. doi:10.1051/0004-6361/201425017
  * Seifert et al. (2003) Seifert, W., Appenzeller, I., Baumeister, H., et al. 2003, Proc. SPIE, 4841, 962. doi:10.1117/12.459494
  * Shapley et al. (2015) Shapley, A. E., Reddy, N. A., Kriek, M., et al. 2015, ApJ, 801, 88. doi:10.1088/0004-637X/801/2/88
  * Shields et al. (2010) Shields, G. A., Ludwig, R. R., & Salviander, S. 2010, ApJ, 721, 1835. doi:10.1088/0004-637X/721/2/1835
  * Simcoe et al. (2013) Simcoe, R. A., Burgasser, A. J., Schechter, P. L., et al. 2013, PASP, 125, 270. doi:10.1086/670241
  * Stefanon et al. (2017) Stefanon, M., Yan, H., Mobasher, B., et al. 2017, ApJS, 229, 32. doi:10.3847/1538-4365/aa66cb
  * Steidel et al. (2016) Steidel, C. C., Strom, A. L., Pettini, M., et al. 2016, ApJ, 826, 159. doi:10.3847/0004-637X/826/2/159
  * Storchi-Bergmann et al. (2009) Storchi-Bergmann, T., McGregor, P. J., Riffel, R. A., et al. 2009, MNRAS, 394, 1148. doi:10.1111/j.1365-2966.2009.14388.x
  * Storey & Hummer (1988) Storey, P. J. & Hummer, D. G. 1988, MNRAS, 231, 1139. doi:10.1093/mnras/231.4.1139
  * Sutherland et al. (2018) Sutherland, R., Dopita, M., Binette, L., et al. 2018, Astrophysics Source Code Library. ascl:1807.005
  * Takada et al. (2014) Takada, M., Ellis, R. S., Chiba, M., et al. 2014, PASJ, 66, R1. doi:10.1093/pasj/pst019
  * Thomas et al. (2016) Thomas, A. D., Groves, B. A., Sutherland, R. S., et al. 2016, ApJ, 833, 266. doi:10.3847/1538-4357/833/2/266
  * Thomas et al. (2018) Thomas, A. D., Dopita, M. A., Kewley, L. J., et al. 2018, ApJ, 856, 89. doi:10.3847/1538-4357/aab3db
  * Topping et al. (2020) Topping, M. W., Shapley, A. E., Reddy, N. A., et al. 2020, MNRAS, 495, 4430. doi:10.1093/mnras/staa1410
  * Trump et al. (2011) Trump, J. R., Weiner, B. J., Scarlata, C., et al. 2011, ApJ, 743, 144. doi:10.1088/0004-637X/743/2/144
  * Veilleux (2002) Veilleux, S. 2002, IAU Colloq. 184: AGN Surveys, 284, 111. doi:10.48550/arXiv.astro-ph/0201118
  * Veilleux & Osterbrock (1987) Veilleux, S. & Osterbrock, D. E. 1987, ApJS, 63, 295. doi:10.1086/191166
  * Vito et al. (2018) Vito, F., Brandt, W. N., Yang, G., et al. 2018, MNRAS, 473, 2378. doi:10.1093/mnras/stx2486
  * Vladilo et al. (2011) Vladilo, G., Abate, C., Yin, J., et al. 2011, A&A, 530, A33. doi:10.1051/0004-6361/201016330
  * Wang et al. (2017) Wang, W., Faber, S. M., Liu, F. S., et al. 2017, MNRAS, 469, 4063. doi:10.1093/mnras/stx1148
  * Whitaker et al. (2014) Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104. doi:10.1088/0004-637X/795/2/104
  * Wuyts et al. (2016) Wuyts, E., Wisnioski, E., Fossati, M., et al. 2016, ApJ, 827, 74. doi:10.3847/0004-637X/827/1/74
  * Yung et al. (2021) Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., et al. 2021, MNRAS, 508, 2706. doi:10.1093/mnras/stab2761

##  Appendix A Flux calibration of NIRSpec spectra

###  A.1 Comparison of Balmer and Paschen line ratios

Figure 15: Diagram comparing the H$\beta$/H$\alpha$ emission line ratio to the
Pa$\beta$/H$\alpha$ ratio (top panel) and the Pa$\gamma$/H$\alpha$ ratio
(bottom panel), for CEERS galaxies in the redshift range $1<z<3$. The
predictions of different dust attenuation models are shown with colored lines.
The attenuation AV along the line is specified for the Calzetti law.

In this appendix, we show a representative test performed to assess the
quality of the relative flux calibration of our fully reduced spectra across
the entire wavelength range from $1\mu m$ to $5\mu m$. In this test, we take
the recombination lines in our NIRSpec data, and that are available for the
highest number of sources, namely H$\alpha$, H$\beta$, Pa$\beta$, and
Pa$\gamma$. In Fig. 15, we compare the optical line ratio H$\beta$/H$\alpha$
and the Balmer-Paschen ratio Pa$\beta$/H$\alpha$ (or Pa$\gamma$/H$\alpha$) to
the predictions of different dust attenuation laws, including the Calzetti et
al. (2000) law, an SMC law from Bouchet et al. (1985), and two mixed model
predictions from Calabrò et al. (2018).

Overall, we can see that the observed points occupy a region of the diagrams
that is physically allowed by one or more of the above attenuation models,
with most galaxies showing relatively low dust attenuations and located in the
bottom right part. Other galaxies have slightly different properties and
higher attenuations AV, being consistent with either the Calzetti law, an SMC
law, or the mixed model. Nearly all outliers in Fig. 15 have at least one line
that is undetected (usually H$\beta$ or Pa$\gamma$), and thus their line
ratios are upper limits and cannot be used for constraining the dust
attenuation. As a result, the absence of a significant number of unphysical
line ratios, or systematic deviations from model predictions, indicates that
the relative flux calibration of the spectra across $2\mu m$ in wavelength can
be trusted on average. However, we still suggest that each object should be
checked individually before computing ratios of distant lines or using their
line luminosities for scientific purposes.

##  Appendix B Comparison with other AGN model predictions

Figure 16: BPT diagram, analog of Fig. 7, where the line ratio predictions for
AGNs are obtained using the recent models of Thomas et al. (2016, 2018). For
AGN models with the same ionization parameter (colored as indicated in the
legend), the circles are the predictions obtained with E${}_{\text{peak}}=20$
eV, while the square symbols are derived assuming E${}_{\text{peak}}=100$ eV.
The marker size varies as a function of gas density (from $10^{2}$ to $10^{4}$
cm-3 from the smaller to the larger). The three points with the same marker,
size, and color, are the predictions of three different values of pNT: $0.1$,
$0.25$, and $0.4$. The lines, markers, and colors for the star-forming models
and observations are the same as in Fig. 7.

Figure 17: Near-infrared diagnostic diagrams including Fe2S3-$\beta$,
Fe2S3-$\alpha$, P2S3, and C1S3. This figure is the analog of Fig. 8, but the
emission line predictions for the NLR of AGNs are obtained using the AGN
models described in Thomas et al. (2016, 2018). The markers, sizes, and colors
of AGN models are the same as described in Fig. 16. The star-forming models
are the same as in Fig. 8. The maximum starburst and AGN lines are defined in
Section 4.1 and in Table 3.

We analyze in this appendix all the diagnostic diagrams presented in this
paper, using emission line predictions obtained with alternative AGN SEDs
compared to our default spectral shape (Risaliti et al., 2000). In particular,
we analyze the more recent AGN spectral shapes introduced by Thomas et al.
(2016), which depend on three parameters: E${}_{\text{peak}}$, pNT, and
$\Gamma$, as discussed in Section 3.2. We first analyze the predictions for
the [O iii]/H$\beta$ and [N ii]/H$\alpha$ line ratios (BPT diagram), which are
shown in Fig. 16. We can see that the model predictions, at fixed ionization
parameter, gas density, and metallicity, occupy a similar region compared to
those obtained with our default AGN SED (Fig. 7). Within the same family where
those three parameters ($\log U$, $Z$, and $n_{e}$) are fixed, we can notice
smaller variations as due to the specific properties of the ionizing shape.
The stronger variation is produced by E${}_{\text{peak}}$, for which we show
for clarity only the predictions obtained with the extreme values considered
in this work (E${}_{\text{peak}}=$ $20$ eV and $100$ eV). In detail, a higher
E${}_{\text{peak}}$ yields on average a higher [O iii]/H$\beta$ and [N
ii]/H$\alpha$ ratios, with the increase rate depending on $Z$ and $\log U$.
The enhancement is almost negligible at small (i.e., subsolar) metallicities
and small ionization parameters ($\log U\leq 3.5$), but it can be up to $+0.2$
dex at supersolar metallicity and higher $\log U$. This effect was already
discussed in Thomas et al. (2016), and is due to more energetic and ionizing
photons in SEDs with higher E${}_{\text{peak}}$, which mimics the enhancement
of the ionization parameter and gas temperature.

On the other hand, pNT produces much smaller variations in the predicted line
ratios. A higher pNT increases the fraction of energetic photons in the X-ray
part of the spectrum (Fig. 6), including the soft X-rays that are more easily
absorbed by the nebula. This causes a more extended partially ionized zone and
more emission coming from this region, especially from low-ionization species.
However, pNT does not significantly alter the photoionization model
predictions in the BPT diagram, with variations of the line ratios typically
smaller than $0.1$ dex, even smaller than the effects due to varying gas
density in the nebula. In Fig. 16 and all subsequent figures, model
predictions with different E${}_{\text{peak}}$ are represented with different
symbols (circles and squares), while predictions corresponding to different
values of pNT are drawn with the same symbol and are almost overlapping.

We also analyze the predictions of the OXAF models for all the near-infrared
diagrams introduced in Section 4.1. The results are presented in Fig. 17,
which is the analog of Fig. 8. Similarly to the BPT diagram, also in the near-
IR diagnostics the models of Thomas et al. (2016) produce line ratios that are
comparable to those obtained with our default AGN SED, occupying the same
parameter space with the same limits on [S iii]/Pa$\gamma$, [S iii]/Pa$\beta$,
[Fe ii]/Pa$\beta$, [Fe ii]/Pa$\alpha$, [P ii]/Pa$\beta$, and [C i]/Pa$\beta$.
If we exclude the influence of metallicity and ionization parameter,
E${}_{\text{peak}}$ has the third major effect on the emission line
intensities. In particular, choosing models with higher E${}_{\text{peak}}$
increases all the line ratios listed above, moving all the points toward the
upper right part of the plot, by factors that depend on the metallicity,
ionization parameter, and on the specific diagram. In detail, for models with
$\log U>-2$, the line ratios on both axis increase by $0.2$-$0.5$ dex in the
range $20<E_{\text{peak}}/eV<100$, because of more energetic ionizing photons
available at higher E${}_{\text{peak}}$. In the last bin of $\log U$ shown in
Fig. 17 ($\simeq-2$), the [S iii]/Pa$\beta$ and [S iii]/Pa$\gamma$ ratios
saturate with E${}_{\text{peak}}$, as higher ionization states of sulfur are
populated. On the other hand, the ratios on the x-axis still increase, as
those low-ionization species are very sensitive to the higher contribution of
soft X-ray photons from models with enhanced E${}_{\text{peak}}$ (see Fig. 6),
as also discussed in Section 5.1. As shown by Thomas et al. (2016), a similar
mechanism is responsible for the higher [O i] and [S ii] emission in the
optical.

The influence of pNT on the same near-IR line ratios is in general smaller
than that produced by E${}_{\text{peak}}$, with variations $<0.1$ dex for
models with low and intermediate ionization ($\log U<-2$) in all the diagrams
in Fig. 17, except [C i]/Pa$\beta$, which can vary by up to $\sim 0.3$ dex
from p${}_{NT}=0.1$ to $0.4$. Also for this parameter, at the highest $\log U$
considered here, it mostly affects the line ratios on the x-axis by the same
mechanism explained above for E${}_{\text{peak}}$, with increments of up to
$\sim 1$ dex for [C i]/Pa$\beta$ and $\sim 0.5$ dex for the remaining ratios.
We note that [C i]/Pa$\beta$ is more sensitive to variations of
E${}_{\text{peak}}$ and pNT because it is entirely produced in partially
ionized regions where electron collisions stimulated by non-thermal continuum
radiation are the main excitation mechanism.

Finally, as shown by Thomas et al. (2016), the third parameter $\Gamma$ has a
similar behavior to pNT, but produces even smaller variations of the line
ratios in all the diagrams analyzed above (both in the optical and near-IR),
so we have fixed it to $+2.0$ in the simulations. We note that a more
exhaustive discussion of the effects of each parameter can be found in the
introductory paper of the OXAF models. We also remark that the exact line
ratio predictions of AGN models do not affect the maximum starburst separation
line (derived from star-forming models), hence they do not alter the results
of this paper.