11institutetext: Department of Statistics, LMU Munich 22institutetext:
Institute of Statistics, RWTH Aachen 33institutetext: ICB, Helmholtz Zentrum
Munich 44institutetext: CCM, Flatiron Institute 55institutetext: Institute of
Behavioral Science and Technology, University of St.Gallen

# Factorized Structured Regression
for Large-Scale Varying Coefficient Models

David Rügamer 1122 Andreas Bender 11 Simon Wiegrebe 11 Daniel Racek 11
Bernd Bischl 11 Christian L. Müller 113344 Clemens Stachl 55

###### Abstract

Recommender Systems (RS) pervade many aspects of our everyday digital life.
Proposed to work at scale, state-of-the-art RS allow the modeling of thousands
of interactions and facilitate highly individualized recommendations.
Conceptually, many RS can be viewed as instances of statistical regression
models that incorporate complex feature effects and potentially non-Gaussian
outcomes. Such structured regression models, including time-aware varying
coefficients models, are, however, limited in their applicability to
categorical effects and inclusion of a large number of interactions. Here, we
propose Factorized Structured Regression (FaStR) for scalable varying
coefficient models. FaStR overcomes limitations of general regression models
for large-scale data by combining structured additive regression and
factorization approaches in a neural network-based model implementation. This
fusion provides a scalable framework for the estimation of statistical models
in previously infeasible data settings. Empirical results confirm that the
estimation of varying coefficients of our approach is on par with state-of-
the-art regression techniques, while scaling notably better and also being
competitive with other time-aware RS in terms of prediction performance. We
illustrate FaStR’s performance and interpretability on a large-scale
behavioral study with smartphone user data.

##  1 Introduction

From buying products online to selecting a movie to watch, recommender systems
(RS) are part of our everyday life. RS are used to suggest those items that
are most appealing to a given user based on the user’s past preference data or
the similarity of a user to other users. One big advantage of RS is their
scalability, as they allow for modeling thousands of interactions, e.g.,
between users and items, and thereby facilitate individual recommendations in
many dimensions (see, e.g., [46] for a recent implementation framework). Many
RS can be represented as a regression model with the user and item as
covariates. This makes it straightforward to include further features into the
model and extend the method by other structural components.

At the same time, the increasing amount of available data and the possibility
to model increasingly complex data generating processes calls for efficient
methods to fit flexible regression models on large-scale data sets with many
observations and features. In the past, several advanced statistical
regression models have been proposed to incorporate complex feature effects.
One of the most common approaches are generalized additive models (GAMs),
widely considered to be state-of-the-art (SotA) for statistical modeling [39].
These models allow the incorporation of time-varying feature effects and
spatial effects, among others, and can also deal with non-Gaussian outcomes
(see [39] for more details). While well-working adaptions of GAMs for large
data scenarios exist [41], both methodology and software reach their limits
when modelling categorical effects or categorical interactions of several
variables where features comprise hundreds or thousands of categories. An
amalgamation of methods from RS and statistical regression can overcome the
limitations of statistical regression models on large-scale data sets with
many categorical effects and interactions. In this work, we combine smoothing
approaches with factorization terms to overcome the limitations of varying
coefficient models for categorical features with many factor levels. Our idea
arises from the statistical analysis of a large-scale behavioral dataset
analyzed in Section 7. In this dataset, smartphone usage behavior of
participants was tracked for several months. From this dataset, continuous
user activity values aggregated over certain time periods can be obtained
(e.g., average screen time in the morning, afternoon, evening and night) for
every user and every activity. Domain experts are interested in various
structured regression effects, such as the continuous activity levels over
time. While standard regression software allows to fit some of these effects,
several hundred activities and users make it infeasible to fit a model that
learns interaction effects of users and activities or smoothly varying time
effects for one or both of these categorical variables.

#### Our contribution

We propose Factorized Structured Regression (FaStR) for scalable varying
coefficient models. This combined approach has the flexibility that has proven
successful in additive regression models, while also being able to deal with
high-dimensional categorical effects and interactions. More specifically, we
1) derive a general model formulation in (4) to combine GAMs and factorization
approaches, 2) derive a varying factorization interaction in Section 4.1 that
reduces the number of parameters and therefore computations by a factor of
$(I+U)/(IU)$ for given numbers of category levels $I$ and $U$, and 3) propose
an efficient implementation of this fusion approach that a) can reduce the
storage cost by a factor of $1/(I+U)$ and b) circumvents computations
quadratic in $I$ and $U$ by using stochastic optimization, an array
reformulation, as well as dynamic feature encoding. In numerical experiments,
we moreover show that our approach 4) leads to an estimation performance
comparable with a SotA implementation and 5) has the desired computational
complexity. Finally, we 6) demonstrate its interpretability and applicability
to large-scale data sets, where it is also competitive with existing RS
approaches.

##  2 Related Literature

Multiple different RS have been proposed over the last years, many of them
based on matrix factorization (MF) [35, 25] or collaborative filtering (see,
e.g., [37]). While recent methods increasingly rely on neural network-based
factorization or recommendation (see, e.g., [42]), it remains debatable
whether they yield superior results, e.g., with respect to performance and
efficiency [21, 30]. Factorization Machines (FM) represent another line of
research which is closely related to MF. Initially proposed by [29], FM are
based on a linear model problem formulation with pairwise (or order $d\geq 3$)
interactions between all features. In particular, the formulation as
regression model is the basis for extensions to other (prediction) tasks, with
many different FM-type models having been developed in recent years (e.g., [5,
19]). An important influence (context variable) in RS is time. Various methods
for controlling for short- and long-term temporal dynamics, cyclic patterns,
drift or time decay exist [2]. While short-term approaches either divide time
into smaller periods or integrate time features into the factorization of the
neighbourhood, long-term effects are accounted for by some form of distance
calculations between the current and other designated time points. Some
approaches also combine the factorization with the time context, and for
instance assume smoothness in the factorization, e.g., for video pixel
completion [20]. Specific time-aware methods include collaborative filtering
with temporal dynamics [24], dynamic MF [10], temporal regularized MF [43], or
sequence-aware FM [8]. The common ground of these methods is to account for
the time context in the factorization.

#### Statistical approaches and interpretability

Several combinations of statistical approaches and RS have been proposed in
past years. Already in 1999, [11] proposed a (Bayesian) generalized mixed-
effects model as RS Likelihood approximation approach. Multiple probabilistic
versions of matrix factorization have also been proposed (see, e.g., [34]).
[44] proposed the GLMix model that combines the idea of generalized linear
models (GLMs; [27]) with RS for large-scale response prediction. Latest
research puts strong emphasis on understanding and controlling of RS through
improved model explainability [14]. Our work is most similar to RS approaches
which facilitate interpretability by making connections to generalized
additive models (GAMs) [13, 47]. In contrast to our approach, however, past
work does not include smoothing splines directly into the models, nor does it
address varying coefficient models. An exception is the time-varying tensor
decomposition by [45] which is inspired by varying coefficient models. While
similar in motivation, their work does not focus on scaling aspects compared
to classical approaches and only considers approximate varying coefficients
with separately learned basis coefficients. Our approach implements the full
varying coefficient model with exact single-varying coefficients as well as a
doubly-varying coefficient with jointly-learned latent basis coefficients.

##  3 Background

We will first describe the necessary background on factorization approaches in
RS, structured additive regression models, and introduce our notation.

###  3.1 General Notation

In the following, we use $Y\sim\mathcal{F}$ to denote a random outcome value
(e.g., a rating) from distribution $\mathcal{F}$ and its observation
$y\in\mathcal{Y}\subseteq\mathbb{R}$ for which the model generates a
prediction $\hat{y}$. We reserve the indices $i\in\mathcal{I}$ and
$u\in\mathcal{U}$ for two categorical features (exemplarily referred to as
_item_ and _user_) and $t\in\mathcal{T}$ for the context variable _time_ on a
given time interval $\mathcal{T}$. The features associated with $i$ and $u$
are assumed to be binary indicator variables and are only implicitly
referenced using their index. In Section 3.3, we however use an integer
representation to introduce a memory-efficient storage representation. Other
context features are summarized by $\bm{x}\in\mathbb{R}^{p}$. We use $b$, $w$
and $v$ to denote weights in the model that relate to items and users, and
make their dependence explicit by indexing these weights correspondingly with
$i$ and $u$. To distinguish between dependencies of categorical features and
the (continuous) feature $t$, we highlight time-dependency by writing objects
as functions of the time $t$. We assume that we are given a dataset
$\mathcal{D}=\\{(y_{iu}(t),\bm{x}_{iu}(t))\\}_{i\in\mathcal{I},u\in\mathcal{U},t\in\mathcal{T}}$
of total size $N$ and allow observations to be sparse, i.e., for $\mathcal{D}$
to be a true subset of $\mathcal{I}\times\mathcal{U}\times\mathcal{T}$. For
matrix computations in later sections, let
$\bm{A}\odot\bm{B}\in\mathbb{R}^{N\times a\cdot b}$ define the row-wise tensor
product (RWTP) of matrices $\bm{A}\in\mathbb{R}^{N\times a}$ and
$\bm{B}\in\mathbb{R}^{N\times b}$, i.e., a Kronecker product applied to every
pair of rows of both matrices
$[\bm{A}_{[1,]}\otimes\bm{B}_{[1,]}\ldots\bm{A}_{[N,]}\otimes\bm{A}_{[N,]}]$.
Further, for $a=b$, let $\bm{A}\bullet\bm{B}:=(\bm{A}*\bm{B})\bm{1}_{b}$,
where $*$ is the Hadamard product and $\bm{1}_{b}\in\mathbb{R}^{b}$ a vector
of ones. The operation defined by $\bullet$ can be exploited in models with
Kronecker product structures such as array models for fast computation (see
Section 4.2).

###  3.2 Model-based Recommender Systems

The basic MF model generates its predictions as

$\hat{y}_{iu}=\langle\bm{v}_{1,i},\bm{v}_{2,u}\rangle$ (1)

using a dot product
$\langle\bm{v}_{1,i},\bm{v}_{2,u}\rangle=\bm{v}_{1,i}^{\top}\bm{v}_{2,u}$ of
two latent factors $\bm{v}_{1,i},\bm{v}_{2,u}\in\mathbb{R}^{D}$ from a
$D$-dimensional joint latent factor space. After learning the mapping from
each item $i$ and user $u$ to the respective latent factor vector, the dot
product describes the interplay between user and item and is used to estimate
the outcome $y$ (ratings). If the combination $\mathcal{I}\times\mathcal{U}$
is observed completely, common matrix decomposition approaches such as a
singular value decomposition can be applied. If the matrix containing the
ratings for all user-item combinations is sparse, missing values can be
imputed. This, however, can be inaccurate and computationally expensive. The
common alternative is to use (1) to only model
$(i,u)\in\mathcal{K}\subset(\mathcal{I}\times\mathcal{U})$ in the set
$\mathcal{K}$ of observed combinations. The solutions $\hat{y}$ can be found
by least squares estimation where an additional $L_{2}$-penalty for
$\bm{v}_{i}$ and $\bm{v}_{u}$ is typically added to the objective function
[25]. In order to account for systematic user- and item-level trends, biases
are further added to (1), yielding

$\hat{y}_{iu}=\mu+b_{i}+b_{u}+\langle\bm{v}_{1,i},\bm{v}_{2,u}\rangle,$ (2)

where $\mu$ is a global intercept representing the average rating, and
$b_{i},b_{u}$ represent the item and user tendencies. The latter two bias
terms are again penalized using a ridge penalty. Together with this
penalization, $b_{i}$ and $b_{u}$ can also be interpreted as a random effect
for the item and user (see, e.g., [44]).

#### Time-aware recommender systems

Contexts such as the location or time in which data has been observed can make
a crucial difference (see, e.g., [3]). RS therefore often include a context
dependency. One of the most common context-aware RS are time-aware RS [7].
Time-aware model-based approaches assume the following relationship

$\hat{y}_{iu}(t)=\mu+b_{i}(t)+b_{u}(t)+\bm{v}_{2,u}^{\top}(t)\bm{v}_{1,i},$
(3)

where both biases and the dot product are time-dependent. The rationale behind
a time-varying latent user effect is that users change their behavior over
time, whereas influences of items should be time-independent [25, 7]. While
time is often assumed to be continuous, categorical time-aware models are used
if time information is represented as discrete contextual values.

The time-varying latent user effect in (3) has a similar role as varying-
coefficients in structured additive regression discussed in the following
section.

###  3.3 Structured Additive Regression and Varying Coefficient Model

In statistical modeling, structured additive regression is a technique for
estimating the relationships between a dependent variable (outcome value) and
one or more independent variables (predictors, features). While the most
common form of regression follows an additive structure as introduced in (2)
and (3), in particular including linear effects $\bm{x}^{\top}\bm{\beta}$ of
features $\bm{x}$ or pairwise interactions, factorization terms are usually
not present in classical regression models. Instead, to adapt models for
complex relationships between features and outcome value, smooth, non-linear,
additive terms $f(\cdot)$ of one or more features are included into the model.
These terms are represented by a linear combination of $L$ appropriate basis
functions. A univariate non-linear effect of feature $z$ is, e.g.,
approximated by $f(z)\approx\sum_{l=1}^{L}{B}_{l}(z){w}_{l}$, where
${B}_{l}(z)$ is the $l$-th basis function (such as regression splines,
polynomial bases or B-splines) evaluated at $z$ and $w_{l}$ is the
corresponding basis coefficient. Similarly, tensor product basis
representations allow for two- or moderate-dimensional non-linear
interactions. It is also possible to represent discrete spatial information or
cluster-specific effects in this way (see, e.g., [41]).

One important part of additive regression is the so-called _varying
coefficient_ model [17]. The rationale for these models is the same as for
time-varying RS: effects of features in the model naturally vary over time.
Therefore these models include effects $xf(t)$, such that the effect
(coefficient) of $x$ is given by $f$ evaluated at time $t$, and $f$ is
estimated from the data. A special case is a varying coefficient
$f_{i}(t),i=1,\ldots,I$ where a separate function $f_{i}$ is estimated for all
$I$ levels of a categorical feature. Existing software to model varying
coefficients with smooth time-effects is, however, not scalable to features
with many categories. The bottleneck is an RWTP of the matrix of evaluated
basis functions $\bm{B}:=(B_{1}(t)\ldots B_{L}(t))$ and a (one hot-)encoded
matrix for a categorical variable (e.g., item with $I$ levels).

#### Computational complexity

Assuming equal number of basis functions $L$ for every smooth term $f_{i}$ in
a varying coefficient model with $i=1,\ldots,I$ levels and $N$ observations,
the storage required for the model matrix is $\mathcal{O}(NLI)$ and the
computations $\mathcal{O}(N(LI)^{2})$ (cf. [40]). Similar, for a model with an
interaction effect of, e.g., item and user (with $U$ levels), the storage is
$\mathcal{O}(NLIU)$ and $\mathcal{O}(N(LIU)^{2})$.

##  4 Factorized Structured Regression

In order to address the computational limitations of statistical regression
techniques, we will first introduce the general idea to obtain predictions
from a FaStR model and then go into more specific details and merits. We use
the RS notation to define the model by means of a classical recommendation
task with items, users, time context and an outcome $y$ such as a rating. As
in a typical regression setting, further features $\bm{x}$ might exist that
the modeler is interested in. We assume that the outcome $Y_{iu}(t)$ and all
features $\bm{x}_{iu}(t)$ are observed on a grid for item $i=1,\ldots,I$, user
$u=1,\ldots,U$, time $t=1,\ldots,T$. While our approach also works for
sparsely and irregularly observed data, we assume a grid of observations to
simplify the notation. Conditional on the item, user, time and further
features, $Y_{iu}(t)$ is assumed to follow a parametric distribution
$\mathcal{F}$. We model the expectation of $Y_{iu}(t)$ as

$\begin{split}\mathbb{E}(Y_{iu}(t)|\bm{x}_{iu}(t))&=h(\eta_{iu}(t)),\\\
\eta_{iu}(t)&=\mu+b_{i}+b_{u}+b_{iu}+f^{[0]}(t)+f^{[1]}_{i}(t)+f^{[2]}_{u}(t)+f^{[3]}_{iu}(t)\\\
&\quad+\textstyle\sum_{o=1}^{O}g_{o}(\bm{x}_{iu}(t)).\end{split}$ (4)

Here, $h$ is an activation or response function mapping the additive predictor
$\eta_{iu}(t)$ onto the correct domain (e.g., $\mathbb{R}^{+}$ for a positive
outcome variable). All terms indicated with $b$ are (regularized) bias terms.
$\mu$ is a global bias term, $b_{i}$ an item-specific bias, $b_{u}$ a user-
specific bias and $b_{iu}$ an item-user-specific one. Terms denoted by $f$ are
smooth non-linear functions of time $t$ represented by (penalized) basis
functions. These include a global trend $f^{[0]}$, a subject and activity
trend, $f^{[1]}$ and $f^{[2]}$, respectively, and a joint trend $f^{[3]}$.
Additional covariates $\bm{x}$ can be modeled using other (smooth) functions
$g_{o}$. In Supplemental Material 0.A we provide further details on smoothness
penalties and optimization of the model.

The model in (4) can be seen as an alternative notation for a varying
coefficient model, or also as a time-aware RS with additional exogenous terms.
As $\mathcal{F}$ is not required to be Gaussian, it has, however, a more
general applicability (e.g., binary, count or interval data). What further
distinguishes (4) from existing approaches is the smoothness assumption of
terms denoted with $f$, combined with the efficient implementation of terms
$f^{[1]}$, $f^{[2]}$ and a factorization assumption for $f^{[3]}$. These
aspects are explained in more detail in the following.

###  4.1 Varying Factorized Interactions

For high-dimensional data, such as the mobile phone data in our example,
estimating the 2- or 3-way interaction terms is computationally not feasible.
We thus propose to define $\eta_{iu}(t)$ in (4) using latent factorization
representations. We therefore decompose the discrete interaction term(s) into
an inner product

$b_{iu}=\langle\bm{v}_{1,i},\bm{v}_{2,u}\rangle=\textstyle\sum_{d=1}^{D}v_{1,i,d}\cdot
v_{2,u,d}$ (5)

with $D\ll\min(I,U)$ latent factors, resulting in the estimation of
$D\cdot(I+U)$ instead of $I\cdot U$ parameters. If $I=U=1000$ and $D=5$, for
example, this reduces the number of parameters by a factor of $100$ from
$10^{6}$ to $10^{4}$. While this is the common approach to model interactions
in factorization approaches, we here propose to proceed in a similar fashion
to model time-dependent interactions and approximate time-varying interactions
by a factorization approach:

$f^{[3]}_{iu}(t)\approx\tilde{f}^{[3]}(t,\bm{V}_{1,i},\bm{V}_{2,u})=\textstyle\sum_{l=1}^{L}B_{l}(t)\textstyle\sum_{d=1}^{D}v_{1,i,l,d}\cdot
v_{2,u,l,d},$ (6)

where $\bm{V}_{\cdot,\cdot}\in\mathbb{R}^{L\times D}$ are matrices with rows
corresponding to the $L$ basis functions for one categorical effect and
columns to the $D$ latent factors. In other words, we approximate the
interaction of the smooth effect $f(t)$ of $t$ and categorical variables $i,u$
by a product of the non-linear basis $\bm{B}$ of dimension $L$ and the two
latent matrices $\bm{B}^{\top}(\bm{V}_{1,i}\bullet\bm{V}_{2,u})$, which can be
computed efficiently for all $N$ rows. The representation via latent factors
requires the estimation of $L\cdot D\cdot(I+U)$ instead of $L\cdot D\cdot
I\cdot U$ parameters (a multiplicative reduction of $(I+U)/(IU)$ in parameters
and computations). This principle is general and can be applied to various
types of additive effects, also of two or higher dimensions such as tensor-
product splines or Markov random field smooths (see [39]).

#### Penalization

In order to enforce smoothness of the varying coefficients in the time-
dimension, a quadratic Kronecker sum penalty $J$ can be added to the loss
function [39]. In a similar manner, we can promote smoothness of the latent
factors $\bm{V}_{iu}=(\bm{V}_{1,i}\bullet\bm{V}_{2,u})\in\mathbb{R}^{L}$ in
our adaption using a symmetric difference penalty matrix
$\bm{P}_{V}\in\mathbb{R}^{L\times L}$. $\bm{P}_{V}$ penalizes the time-
dimension of the factorized varying-coefficients, where the penalized
differences (its entries) depend on the chosen basis $B$. We further allow for
an $L_{2}$-regularization of the latent factors in the $i$\- and
$u$-dimension, yielding

$J=\lambda_{t}\cdot\bm{V}_{iu}^{\top}\bm{P}_{V}\bm{V}_{iu}+\lambda_{iu}\cdot(||\bm{V}_{1,i}||_{F}+||\bm{V}_{2,u}||_{F}),$
(7)

where $\lambda_{t}$ controls the smoothness of the non-linearity in the
direction of the time $t$, $||\cdot||_{F}$ is the Frobenius norm, and
$\lambda_{iu}$ the regularization for items and users.

###  4.2 Efficient Implementation

While the previous section allows to efficiently model (smooth) interactions
of two or more categorical features with many categories, the factorization is
not a solution for coefficients $f^{[1]}$ and $f^{[2]}$ varying with a single
category. Many use cases also require to estimate one effect for each of the
levels of a categorical (interaction) effect. One bottleneck if $I$ and/or $U$
is large, is the computation of their dummy-encoded design matrices. For
example, for item the matrix $\bm{X}_{I}$ of size $N\times I$ contains binary
entries indicating which observation (rows) belong to which item category
(column), and analogous for a user matrix $\bm{X}_{U}$. Second, the execution
of operations involving such large matrices is another computational
bottleneck. Computations become even more challenging if the model includes
interactions, resulting in $\mathcal{O}(NIU)$ storage and
$\mathcal{O}(N(IU)^{2})$ computations (cf. Section 3.3). These interactions
are created by calculating the RWTP between both matrices, i.e.,
$\bm{X}_{I}\odot\bm{X}_{U}\in\mathbb{R}^{N\times I\cdot U}$. To circumvent
creating, storing and processing $\bm{X}_{I}$ and $\bm{X}_{U}$ as a whole, we
propose two simple yet effective implementation tricks explained in the
following.

Stochastic optimization The first bottleneck in computations of varying
coefficient models at scale is the number of observations $N$. We therefore
implement FaStR in a neural network and thereby can make use of stochastic
gradient descent optimization routines with mini-batches of size $M\ll N$.
This reduces the original cost of computations from $\mathcal{O}(N(IU)^{2})$
to $\mathcal{O}(EMIU)$ where $E$ is the number of model updates. It also
allows us to leverage high-performance computing platforms such as TensorFlow
[1] that support GPU computations.

Array computations The second bottleneck is computing the RWTP of two- or
higher-dimensional interaction terms. Our proposal is to use an array
reformulation that does not require to compute the RWTP design matrix in the
first place. More specifically, a two-dimensional interaction effect
$(\bm{X}_{I}\odot\bm{X}_{U})\bm{w}$ with weights $\bm{w}$ can be equally
represented by

$(\bm{X}_{I}\bm{W})\bullet\bm{X}_{U},$ (8)

where $\bm{W}\in\mathbb{R}^{I\times U}$ is a matrix of weights with the $i$th
row and $u$th column being the weight for the interaction of the $i$th level
in $\bm{X}_{I}$ and the $u$th level in $\bm{X}_{U}$. By using (8) instead of a
plain linear effect, we circumvent the construction of the large RWTP and the
storage cost is reduced from $\mathcal{O}(NIU)$ to $\mathcal{O}(N(I+U))$
without increasing the time complexity (as the $\bullet$ operation for $M$
observations is neglectable with $\mathcal{O}(MU)$ compared to the matrix
multiplication with $\mathcal{O}(MIU)$). This array formulation can also be
defined for higher-order interactions [12].

Dynamic feature encoding Although array computations can reduce the storage
problem notably by not constructing the RWTP in the first place, a third
bottleneck is storing the large dummy-encoded matrices $\bm{X}_{I}$ and
$\bm{X}_{U}$ themselves. We circumvent this extra space complexity, by
evaluating categorical features dynamically during network training and only
constructing the one-hot encoding for categorical features on a given mini-
batch $m$. Thereby, only a matrix $\bm{X}^{(m)}_{I}$ of size $M\times I$ needs
to be loaded into memory and the full $N\times I$ matrix is never created
explicitly. This effectively reduces the storage from $\mathcal{O}(NIU)$ to
$\mathcal{O}(N)$ (exactly $2N$ for two integer vectors). While this
potentially results in redundant computations as it will create the encoding
for a specific observation multiple times (if the number of epochs is greater
than 1), deterministic integer encoding is cheap. Hence, the resulting
computation overhead is usually neglectable and both forward- and backward
pass can make use of the sparse representation.

When evaluating the two matrix operators (matrix product and $\bullet$) in (8)
sequentially, encodings can again be created dynamically and the largest
matrix involved only contains $\max(I,U)$ instead of $I\cdot U$ columns.

##  5 Numerical Experiments

Our numerical experiments investigate 1) whether FaStR can estimate
(factorized) varying coefficients as proposed in Section 4.1 with performance
comparable to other SotA methods, and 2) whether the model presented in (4)
can be initialized and fitted with constant memory scaling w.r.t. the number
of factor levels. In addition, in the Supplementary Material 0.B we
investigate whether FaStR can recover GAMs in general using our implementation
techniques. Details on the data generating processes for each simulation
study, method implementation and the computing environment are given in the
Supplementary Material 0.C. In the first experiment, we compare the estimated
and true model coefficients using the mean squared error (MSE) and the
estimated vs. the true functions $f$ using the squared error, integrated and
averaged over the domain of the predictor variables (MISE). We repeat the data
generating process 10 times to account for variability in the data and model
estimation.

###  5.1 Estimation of Factorized Smooth Interactions

We first investigate how well our models can recover smooth factorized terms
from Section 4.1 of two categorical variables $j$ and $u$, with 4 and 5
levels, respectively, i.e., 20 different smooth effects. We use a relatively
small number of levels to be able to fit all possible interactions also in a
classic structured regression approach. We define one of the two true
factorized smooth effects as a varying coefficient term, i.e., $f_{1,iu}(t)$,
and one stemming from an actual factorization, i.e.,
$f_{2}(t,\bm{V}_{1,i},\bm{V}_{2,u})$. We use the Bernoulli and Gaussian
distributions and investigate factorized terms in the distributions’ mean for
different data sizes $N\in\\{500,1000,2000\\}$. While the true model is
generated using three latent factor dimensions, we additionally investigate a
model with six dimensions to see how the misspecification of the latent
dimension size influences estimation performance. FaStR is trained using a
batch size of $250$ with early stopping on 10% of the training data and a
patience of $50$ epochs.

#### Results

All results show that FaStR can estimate varying coefficients equally well
compared to a classic GAM estimation routine. Figure 1 shows the resulting
M(I)SE values for all data settings. Note that the GAM implicitly assumes as
many latent dimensions as there are factor levels, but can also shrink single
smooth functions to zero. In cases where it is feasible to fit smooth effects
for every factor level (as is the case here), GAM can thus be seen as gold
standard. For the case of a normal distribution, we observe that GAM yields
better results, but also that the performance of FaStR converges to the one
from GAM with increasing number of observations. For the Bernoulli case, our
approach benefits from the optimization in a neural network and even
outperforms the classic GAM estimation routine which requires a multiple of
observations compared to a Gaussian setting for good estimation results (as,
e.g., also found in [32]).

Figure 1: Comparison of M(I)SE values for the estimated partial effects of
different methods (colors) for factorized and non-factorized smooth terms
(x-axis), different distributions (columns) and different data sizes (rows).
Values $>1.5$ are set to the value $1.5$ to improve readability.

###  5.2 Memory Consumption

Finally, we compare the memory consumption of our implementation against the
SotA implementation for big additive models (BAM) in the R [28] package mgcv
[38] for an increasing number of category levels (20, 40, 60, 80) when using a
categorical effect or a varying coefficient based on the representation
proposed in Section 4.2, $N\in\\{1000,2000,4000\\}$ and optimization as in
Section 5.1. While the improvement in memory consumption is expected to be
even larger when using factorized terms instead of interaction terms with
weights for each category combination, we do not use factorization in this
experiment as there is no equivalent available in software for additive
regression. Additionally, we also track the time when running FaStR for $10$
epochs to see if there are notable changes in the time consumption for varying
data generating settings.

#### Results

Figure 2 visualizes the results for all different settings and compares run
times and memory consumption of the two methods for factor variables (single)
and varying coefficient effects (varying). Results show that FaStR has both,
almost constant time and memory consumption while the SotA method requires
exponentially more memory for growing numbers of factor levels (as the whole
encoded matrix must be loaded into memory). These results confirm that our
implementation works as intended to allow for the estimation of varying
coefficient models in large-scale settings.

Figure 2: Memory and time consumption (y-axis; $\log_{10}$-scale) comparison
between the SotA big additive model (BAM) implementation and our method (in
different colours) for an increasing number of categorical levels (x-axis) of
a factor effect (single) and varying coefficient term (varying) for different
data sizes (columns).

##  6 Benchmarks

Although the focus of this work is to provide scalable and interpretable
regression models, prediction performance of our models is also of (secondary)
interest. We aim for a similar performance compared to SotA time-aware RS
techniques, yet without the ambition to outperform these methods. We use the
_MovieLens 10M movie ratings_ benchmark data set [16], which is sparse in
terms of user-item combinations, with _items_ corresponding to rated movies.
In addition, we benchmark models on a subset of the densely observed
_PhoneStudy behavior_ data set [36], analyzed in more detail in Section 7. In
both cases, we use single train-test splits (90%/10% and 70%/30%,
respectively) and evaluate the models predictive performance with the root
mean squared error (RMSE). The different characteristics of both data sets are
given in Table 1.

Table 1: Descriptive statistics of the two benchmark data sets. | Movies | Phone Study
---|---|---
# Observations ($N$) |  $\approx 9$ m |  $\approx 8.7$ m
# Users ($U$) | 69,878 | 342
# Movies/Activities ($I$) | 10,677 | 176
# Unique Time Points ($T$) |  $\approx$ 6.5 m | 348

#### Methods

As comparison we use Bayesian timeSVD and Bayesian timeSVD++ flipped, two
variations of the SVD++ method [23], a latent factor model whose key
innovative feature is the incorporation of implicit user information. Both
Bayesian timeSVD and Bayesian timeSVD++ flipped have been extended to be time-
aware [24] and optimized by Gibbs sampling using a Bayesian reformulation
[31]. Bayesian timeSVD++ flipped integrates both implicit user and item
information and has been reported to be the best-performing model among
multiple SotA methods in a recent benchmark study [31]. The second variation,
Bayesian timeSVD, is still a time-aware latent factor model, yet it does not
incorporate implicit user or item information. As we are mainly interested in
the performance of the proposed time-varying coefficient model, the Bayesian
timeSVD provides a much fairer comparison with FaStR as it does not include
the aforementioned types of implicit information. We use tuning parameter
settings as given in [31] for the two benchmark methods (i.e., we use the
already tuned models).

Table 2: Benchmark results based on the RMSE on the respective test data set for both benchmark data sets (columns) for different methods (rows). Best results are in bold with the best results of the respective other method underlined. | Movies | PhoneStudy
---|---|---
timeSVD | 0.872 | 0.089
timeSVD++ flipped | 0.856 | 0.087
FaStR ($D=10$) | 0.984 | 0.076
FaStR ($D=3$) | 0.975 | 0.080
FaStR ($D=1$) | 0.890 | 0.087
FaStR (w/o $f^{[3]}$) | 1.027 | 0.093

We compare these two models against our method as proposed in (4) and thereby
not only test its predictive performance, but also its capability to scale
well to high-dimensional data sets. We do not tune FaStR extensively, but
perform a small ablation study by testing different latent dimensions ($D\in
1,3,10$) for the factorized varying coefficient term $f^{[3]}$ and by
excluding the whole term. All models use early stopping based on a validation
data split with the same proportion as the train-test split.

#### Results

Table 2 shows the performance of all methods. Interestingly, FaStR is
competitive with the SotA timeSVD approaches, even though the premise of this
paper was merely to develop a scalable variant of the varying coefficient
model, not to propose a method with SotA performance on RS tasks.

##  7 User Behavior Phone Study

Figure 3: One dimension of the latent time-varying coefficients for different
activities (left plot with “locked screen” in blue) and users (right plot with
exemplary user in green) over 87 study days (x-axis). Vertical lines show
Christmas and New Year’s Eve.

We finally turn to the motivating case study. A more detailed description of
the study and data set can be found in the Supplementary Material 0.D. To
analyze the activity levels (a value between $0$ and $1$, indicating the
intensity of the activity in the given aggregation window) of participants in
the study, we model the expected activity levels using user and activity
effects, their interaction based on a factorization, an effect of the day of
the week (Mon - Sun), the time of the day (6-hour windows), a factorized
interaction of users and weekday as well as users and daytime, and a smooth
time-dependent study day effect varying by user and/or activity. All
factorizations use a three-dimensional latent space.

#### Results

Results are generally plausible and in line with prior expectations. While
various model effects are examined in the Supplementary Material, we here
briefly analyze the factorized varying coefficient interaction by analyzing
its latent factors. One of the dimensions is depicted in Figure 3. Most
activities follow a global pattern (darker concentration of lines), while a
few show very specific sequences. For example, the “locked screen” event is
observed less in the days around New Year’s Eve. At the user level, no general
pattern is visible (in part due to the different starting dates of the
participants), but we still observe sudden changes in activity around holidays
such as Christmas.

##  8 Conclusion and Outlook

In this work we presented an amalgamation of structured statistical regression
and RS to allow for large-scale varying coefficient models with many
categorical levels. For this, we leveraged factorization approaches combined
with an efficient neural network-based implementation. Empirical results
confirm the efficacy of our proposal. In order to make the proposed approach
as flexible as commonly used statistical regression software, we used a model-
based point of view and cast the approach as a generalized additive model. An
interesting future avenue and additional advantage of formulating the RS model
as structured additive regression is the possibility to extend this approach
to distributional regression, in order to account, e.g., for aleatoric
uncertainty in the data generating process.

## Acknowledgement

This work has been partially supported by the German Federal Ministry of
Education and Research (BMBF) under Grant No. 01IS18036A.

## References

  * [1] Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., Isard, M., et al.: Tensorflow: A system for large-scale machine learning. In: 12th $\\{$USENIX$\\}$ symposium on operating systems design and implementation ($\\{$OSDI$\\}$ 16). pp. 265–283 (2016)
  * [2] Al-Hadi, I., Sharef, N.M., Sulaiman, M.N., Mustapha, N.: Review of the temporal recommendation system with matrix factorization. Int. J. Innov. Comput. Inf. Control 13(5), 1579–1594 (2017)
  * [3] Baltrunas, L., Ricci, F.: Experimental evaluation of context-dependent collaborative filtering using item splitting. User Modeling and User-Adapted Interaction 24(1-2), 7–34 (2014)
  * [4] Baumeister, R., Vohs, K.D., Funder, D.C.: Psychology as the science of self-reports and finger movements: Whatever happened to actual behavior? Perspectives on Psychological Science 2, 396–403 (12 2007)
  * [5] Blondel, M., Fujino, A., Ueda, N., Ishihata, M.: Higher-order factorization machines. In: Lee, D., Sugiyama, M., Luxburg, U., Guyon, I., Garnett, R. (eds.) Advances in Neural Information Processing Systems. vol. 29. Curran Associates, Inc. (2016)
  * [6] Buja, A., Hastie, T., Tibshirani, R.: Linear Smoothers and Additive Models. The Annals of Statistics 17(2), 453 – 510 (1989)
  * [7] Campos, P.G., Díez, F., Cantador, I.: Time-aware recommender systems: a comprehensive survey and analysis of existing evaluation protocols. User Modeling and User-Adapted Interaction 24(1), 67–119 (2014)
  * [8] Chen, T., Yin, H., Nguyen, Q.V.H., Peng, W.C., Li, X., Zhou, X.: Sequence-aware factorization machines for temporal predictive analytics. In: 2020 IEEE 36th International Conference on Data Engineering (ICDE). pp. 1405–1416. IEEE (2020)
  * [9] Chu, E., Keshavarz, A., Boyd, S.: A distributed algorithm for fitting generalized additive models. Optimization and Engineering 14(2), 213–224 (2013)
  * [10] Chua, F.C.T., Oentaryo, R.J., Lim, E.P.: Modeling temporal adoptions using dynamic matrix factorization. In: 2013 IEEE 13th International Conference on Data Mining. pp. 91–100 (2013)
  * [11] Condli, M.K., Lewis, D.D., Madigan, D., Posse, C.: Bayesian mixed-effects models for recommender systems. In: ACM SIGIR. vol. 99 (1999)
  * [12] Currie, I.D., Durban, M., Eilers, P.H.: Generalized linear array models with applications to multidimensional smoothing. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(2), 259–280 (2006)
  * [13] Guo, Y., Su, Y., Yang, Z., Zhang, A.: Explainable recommendation systems by generalized additive models with manifest and latent interactions (2020)
  * [14] Hada, D.V., M., V., Shevade, S.K.: Rexplug: Explainable recommendation using plug-and-play language model. In: Proceedings of the 44th International ACM SIGIR Conference on Research and Development in Information Retrieval. p. 81–91. SIGIR ’21, Association for Computing Machinery, New York, NY, USA (2021)
  * [15] Harari, G.M., Lane, N.D., Wang, R., Crosier, B.S., Campbell, A.T., Gosling, S.D.: Using smartphones to collect behavioral data in psychological science: Opportunities, practical considerations, and challenges. Perspectives on Psychological Science 11, 838–854 (11 2016)
  * [16] Harper, F.M., Konstan, J.A.: The movielens datasets: History and context 5(4) (2015)
  * [17] Hastie, T., Tibshirani, R.: Varying-coefficient models. Journal of the Royal Statistical Society: Series B (Methodological) 55(4), 757–779 (1993)
  * [18] Hester, J.: bench: High Precision Timing of R Expressions (2020), https://CRAN.R-project.org/package=bench, r package version 1.1.1
  * [19] Hong, F., Huang, D., Chen, G.: Interaction-aware factorization machines for recommender systems. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 33, pp. 3804–3811 (2019)
  * [20] Imaizumi, M., Hayashi, K.: Tensor decomposition with smoothness. In: International Conference on Machine Learning. pp. 1597–1606. PMLR (2017)
  * [21] Jin, R., Li, D., Gao, J., Liu, Z., Chen, L., Zhou, Y.: Towards a better understanding of linear models for recommendation. In: Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining. p. 776–785. KDD ’21, Association for Computing Machinery, New York, NY, USA (2021)
  * [22] Kneib, T., Silbersdorff, A., Säfken, B.: Rage against the mean – a review of distributional regression approaches. Econometrics and Statistics (2021)
  * [23] Koren, Y.: Factorization meets the neighborhood: a multifaceted collaborative filtering model. In: Proceedings of the 14th ACM SIGKDD international conference on Knowledge discovery and data mining. pp. 426–434 (2008)
  * [24] Koren, Y.: Collaborative filtering with temporal dynamics. In: Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining. pp. 447–456 (2009)
  * [25] Koren, Y., Bell, R., Volinsky, C.: Matrix factorization techniques for recommender systems. Computer 42(8), 30–37 (2009)
  * [26] Lewin, K.: Defining the ’field at a given time.’. Psychological Review 50, 292–310 (1943)
  * [27] Nelder, J.A., Wedderburn, R.W.: Generalized linear models. Journal of the Royal Statistical Society: Series A (General) 135(3), 370–384 (1972)
  * [28] R Core Team: R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria (2021)
  * [29] Rendle, S.: Factorization machines. In: 2010 IEEE International conference on data mining. pp. 995–1000. IEEE (2010)
  * [30] Rendle, S., Krichene, W., Zhang, L., Anderson, J.: Neural collaborative filtering vs. matrix factorization revisited. In: Fourteenth ACM Conference on Recommender Systems. pp. 240–248 (2020)
  * [31] Rendle, S., Zhang, L., Koren, Y.: On the Difficulty of Evaluating Baselines: A Study on Recommender Systems. arXiv preprint arXiv:1905.01395 (May 2019)
  * [32] Rügamer, D., Kolb, C., Klein, N.: Semi-Structured Deep Distributional Regression: A Combination of Additive Models and Deep Learning. arXiv preprint arXiv:2002.05777 (2020)
  * [33] Ruppert, D., Wand, M., Carroll, R.: Semiparametric Regression. Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press (2003)
  * [34] Shan, H., Banerjee, A.: Generalized probabilistic matrix factorizations for collaborative filtering. In: 2010 IEEE International Conference on Data Mining. pp. 1025–1030 (2010)
  * [35] Srebro, N., Rennie, J.D., Jaakkola, T.S.: Maximum-margin matrix factorization. In: NIPS. vol. 17, pp. 1329–1336. Citeseer (2004)
  * [36] Stachl, C., Au, Q., Schoedel, R., Gosling, S.D., Harari, G.M., Buschek, D., Völkel, S.T., Schuwerk, T., Oldemeier, M., Ullmann, T., Hussmann, H., Bischl, B., Bühner, M.: Predicting personality from patterns of behavior collected with smartphones. Proceedings of the National Academy of Sciences 117, 17680–17687 (7 2020)
  * [37] Thorat, P.B., Goudar, R., Barve, S.: Survey on collaborative filtering, content-based filtering and hybrid recommendation system. International Journal of Computer Applications 110(4), 31–36 (2015)
  * [38] Wood, S.N.: Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1), 3–36 (2011)
  * [39] Wood, S.N.: Generalized additive models: an introduction with R. Chapman and Hall/CRC (2017)
  * [40] Wood, S.N.: Inference and computation with generalized additive models and their extensions. Test 29(2), 307–339 (2020)
  * [41] Wood, S.N., Li, Z., Shaddick, G., Augustin, N.H.: Generalized additive models for gigadata: Modeling the u.k. black smoke network daily data. Journal of the American Statistical Association 112(519), 1199–1210 (2017)
  * [42] Wu, C., Lian, D., Ge, Y., Zhu, Z., Chen, E., Yuan, S.: Fight fire with fire: Towards robust recommender systems via adversarial poisoning training. In: Proceedings of the 44th International ACM SIGIR Conference on Research and Development in Information Retrieval. p. 1074–1083. SIGIR ’21, Association for Computing Machinery, New York, NY, USA (2021)
  * [43] Yu, H.F., Rao, N., Dhillon, I.S.: Temporal regularized matrix factorization for high-dimensional time series prediction. In: Lee, D., Sugiyama, M., Luxburg, U., Guyon, I., Garnett, R. (eds.) Advances in Neural Information Processing Systems. vol. 29. Curran Associates, Inc. (2016)
  * [44] Zhang, X., Zhou, Y., Ma, Y., Chen, B.C., Zhang, L., Agarwal, D.: Glmix: Generalized linear mixed models for large-scale response prediction. p. 363–372. KDD ’16, Association for Computing Machinery, New York, NY, USA (2016)
  * [45] Zhang, Y., Bi, X., Tang, N., Qu, A.: Dynamic tensor recommender systems. Journal of Machine Learning Research 22(65), 1–35 (2021)
  * [46] Zhao, W.X., Mu, S., Hou, Y., Lin, Z., Li, K., Chen, Y., Lu, Y., Wang, H., Tian, C., Pan, X., Min, Y., Feng, Z., Fan, X., Chen, X., Wang, P., Ji, W., Li, Y., Wang, X., Wen, J.R.: Recbole: Towards a unified, comprehensive and efficient framework for recommendation algorithms (2020)
  * [47] Zhuang, H., Wang, X., Bendersky, M., Grushetsky, A., Wu, Y., Mitrichev, P., Sterling, E., Bell, N., Ravina, W., Qian, H.: Interpretable ranking with generalized additive models. In: Proceedings of the 14th ACM International Conference on Web Search and Data Mining. p. 499–507. WSDM ’21, Association for Computing Machinery, New York, NY, USA (2021)

##  Appendix 0.A Penalization and Optimization

The proposed model can be optimized by (restricted) Maximum Likelihood with
negative log-likelihood as loss function and additional difference penalty
terms $J_{0},\ldots,J_{3}$ for model terms $f^{[0]},\ldots,f^{[3]}$ that
penalize differences in neighboring basis coefficients to enforce smoothness
(see, e.g., [39] for more details). In order to determine the amount of
smoothness (the amount of penalization), every penalty term is controlled with
a respective smoothing parameter $\lambda_{0},\ldots,\lambda_{3}$. Next to
tuning these parameters, a possible option is to opt for equal amount of
smoothness of all involved non-linear functions. This can be realized by
setting the _degrees-of-freedom_ $\text{df}_{j}$ for each smooth term
$j\in\\{0,\ldots,3\\}$ to the same value df. As all involved terms are
penalized linear smoothers [6], we can exploit the one-to-one mapping between
$\text{df}_{j}$ and $\lambda_{j}$, and efficiently calculate the different
$\lambda_{j}$ values corresponding to $\text{df}_{j}\equiv\text{df}$ using the
Demmler-Reinsch Orthogonalization (see Appendix B.1.1 of [33]). Given the
different $\lambda$ values, varying coefficient models such as (4) are convex
optimization problems that can be solved efficiently (see, e.g., [9]).

##  Appendix 0.B Further Numerical Experiments

###  0.B.1 General Estimation Performance

The goal of the experiments in this subsection is to show that our framework
is able to represent effects of GAMs and works well when applying
implementation techniques from Section 4.2 together with neural-based
optimization routines. To this end, we simulate GAM models for different data
sizes $N\in\\{1000,2000,4000\\}$ from different distributions (Bernoulli,
Beta, Gaussian, Poisson) with five non-linear smooth terms, a tensor product
smooth, a varying coefficient effect with a categorical variable $i$ with 10
levels, and an interaction term of two categorical variables $j$ and $u$, each
with 10 levels. FaStR is trained using a batch size of $50$ with early
stopping on 10% of the training data and a patience of $100$ epochs.

#### Results

All results indicate that FaStR can estimate the different effects and
distributions equally well when compared to a SotA estimation technique for
GAMs [38]). Figure 4 shows the resulting mean (integrated) squared error
(M(I)SE) values111See Supplementary Material 0.C for a detailed definition.
for all simulation settings, indicating that FaStR’s estimation of the
splines, factor effects and varying coefficients works well and even seems
less vulnerable to difficult data settings (here the estimation of categorical
variable effects for the Bernoulli distribution).

Figure 4: Comparison of model fits of the SotA GAM implementation (in red) and
our method (in gray) for different model terms (columns), different
distributions (x-axis) and corresponding M(I)SE values visualized by boxplots
over 10 runs and different numbers of observations. Outliers (values larger 3)
are set to 3.

##  Appendix 0.C Details on Numerical Experiments

We here describe the details how data is generated in our numerical
experiments. For details on the computing environment see Section 0.E.

In all simulations, we use the same data generating process to create the
additive predictor $\eta_{iu}(t)$. Depending on the distribution $\mathcal{F}$
to be simulated, we transform $\eta_{iu}(t)$ with the respective $h$-function
to generate the distribution’s mean $\theta_{iu}(t)=h(\eta_{iu}(t))$ and draw
outcome values according to $Y_{iu}(t)\sim\mathcal{F}(\theta_{iu}(t))$. The
additive predictor $\eta$ itself is generated by first drawing numerical
features from a standard uniform distribution and categorical features using a
balanced design with $I$ (or $U$) different categories and $\lfloor
I/N\rfloor$ ($\lfloor U/N\rfloor$) observations per category. Based on the
numeric values, univariate splines are generated using the functions
$f_{uni,1}(x)=-(x/4)^{5}$, $f_{uni,2}(x)=\log(x^{2})/10$,
$f_{uni,3}(x)=\sin(3x)$, $f_{uni,4}(x)=(\sin(3x)\cdot 0.1+1)\cdot
I(x<0)+(-2x+1)/4\cdot I(x\geq 0)$ with indicator function $I(\cdot)$ and
$f_{uni,5}(x)=-x\tanh(3x)\cdot\sin(4x)/4$. Bivariate tensor-product splines
are generated using the function $f_{bi}(x,y)=(\pi^{0.3}\cdot
0.4)\cdot(1.2\cdot\exp(-(x-0.2)^{2}/0.09-(y-0.3)^{2}/0.16)+0.8\cdot\exp(-(x-0.7)^{2}/0.09-(y-0.8)^{2}/0.16))$
and the smooth part of varying coefficient splines using one of the following
functions (depending on the category of the factor variable):
$f_{vc,1}(x)=\cos(3x)$, $f_{vc,2}(x)=\tanh(3x)$, $f_{vc,3}(x)=-(x/4)^{3}$,
$f_{vc,4}(x)=\cos(3x-2)\cdot(-x/3)$, $f_{vc,5}(x)=\exp(x/2)-1$,
$f_{vc,6}(x)=(x/2)^{2}$, $f_{vc,7}(x)=\sin(x)\cos(x)$,
$f_{vc,8}(x)=\sqrt{|x|}$, $f_{vc,9}(x)=-(x/4)^{5}$,
$f_{vc,10}(x)=\log(x^{2})/100$.

We evaluate all non-linear functions by taking the point-wise MSE and
integrate over the whole function domain. More specifically, given a function
$f(t)$ and estimate $\hat{f}(t)$ of feature $t\in\mathcal{T}$ with sampled
values $t_{1},\ldots,t_{n}$, the MISE is calculated by
$\mathbb{E}\int_{\mathcal{T}}(\hat{f}(t)-f(t))^{2}dt\approx
n^{-1}\sum_{i=1}^{n}\Delta_{i}(\hat{f}(t_{i})-f(t_{i}))^{2}$, where
$\Delta_{i}$ are integration weights. We do this on the $n$ simulated training
data points as the goal is to evaluate estimation performance.

##  Appendix 0.D Details Phone Study

###  0.D.1 Motivation

In the behavioral sciences and psychology in particular, the understanding and
prediction of behavior is one of the main goals. For a long time, it has been
assumed that behavior can be seen as a function of a person’s internal factors
(e.g., personality traits, past experiences) and external factors (e.g.,
perceived environment) [26]. However, as trivial as this assumption might
seem, it has proven very difficult for researchers to investigate this
assumption empirically [4]. The currently ongoing digital revolution in
behavioral sciences is driven by the increasing availability of large
quantities of digital behavioral data which can be collected directly with
consumer devices. High-dimensional, fine-grained mobile sensing data has
emerged as one of the most promising sources of information to study human
behavior and experience, unobtrusively, in the wild, and at very large scale
[15]. This development requires new methods to handle large-scale data while
also yielding interpretable models. Here, we demonstrate the usefulness of
FaStR in a first analysis of human behavior inferred from a high-dimensional
mobile sensing data.

###  0.D.2 Data set and model specification

The PhoneStudy dataset [36] consists of behavioral metrics that were obtained
with smartphone sensing. The dataset features various activities (mobility,
calendar events, music listening, calling, power charging, app use and more).
Here we use a subset of the dataset that includes up to 87 days of sensing and
that was collected during a study period from October 28th, 2017 to January
22nd, 2018. After pre-processing, we obtain activities from 342 participants
in Germany in 6-hour windows ([0-6], [6-12], [12-18], [18-24]) throughout each
day of the study period.

###  0.D.3 Results

Based on the model’s estimated effects, WiFi events
(connected/disconnected/disabled), airplane mode, Bluetooth (off) and power
charging are the top categories with highest overall (i.e., user-independent)
activities, while certain apps show the lowest overall activity.

Figure 5: Partial effects of the day of the week (compared to Sunday; upper
row) and the time of the day (compared to hours 18-24) on the mean (left
column) and standard deviation (right column).

Results further show that activity of mobile phone usage notably decreases
towards the end of the week (Figure 5) with activity on Saturday being the
lowest, and users are least active in the [0-6]h window (i.e., during sleep),
closely followed by the window [12-18]h (e.g., during work/school). Moreover,
the results of our factorization approach show that some activities are
specific to certain time points (e.g., phone charging is done more often at
the beginning or end of the day, while airplane mode is more often active
during the night). All of these effects are plausible and demonstrate the
usefulness of FaStR to model contextualized human behavior.

##  Appendix 0.E Computing Environment

For comparability reasons, FaStR is implemented in R and only runs on a CPU
(as it is the case for the reference software for GAMs [38]). The details on
computing environments used for both numerical experiments and benchmark
studies are as follows. For all numerical experiments we used the R version
4.1.2, Python 3.7, TensorFlow 2.5.0rc0 on a personal computer with Intel(R)
Core(TM) i7-8665U CPU @ 1.90GHz, 15.3 GB main memory and Ubuntu 20.04.3 LTS
64-bit as operating system. Timings and memory consumptions were tracked using
bench [18].

While we run our approach on the Phone Study data on the same environment as
described in Section 0.E, we use a cloud server with Intel(R) Xeon(R) Gold
6148 CPU @ 2.40GHz, 92.8 GB main memory and Ubuntu 18.04.3 LTS for the Movies
data. We also use the latter environment to run the timeSVD++ models on both
data sets.

##  Appendix 0.F Outlook: Distributional Regression

The presented model in (4) can be generalized as follows:

$Y_{iu}(t)|\bm{x}_{iu}(t)\sim\mathcal{F}({\theta}_{1,ij}(t),\ldots,{\theta}_{K,ij}(t)),$
(9)

where the $K$ parameters $\theta_{k}$ of the distribution $\mathcal{F}$ can be
described by a structured additive (factorized) predictor

$\theta_{k,ij}(t)=h_{k}(\eta_{k,ij}(t)),k=1,\ldots,K.$ (10)

The model described in (9) and (10) is a distributional regression model which
relates some or all distribution parameters to available features in
$\eta_{k},k=1,\ldots,K$ [32, 22]. An interesting property to investigate in
future research is how DR can account for aleatoric uncertainty. This can,
e.g., be useful to check whether the model’s uncertainty varies over different
context variables time.


[a]D. Kostunin

# New insights from old cosmic rays:
A novel analysis of archival KASCADE data

I. Plokhikh  M. Ahlers  V. Tokareva  V. Lenok  P. Bezyazeekov  S. Golovachev
V. Sotnikov  R. Mullyadzhanov  E. Sotnikova

###### Abstract

Cosmic ray data collected by the KASCADE air shower experiment are competitive
in terms of quality and statistics with those of modern observatories. We
present a novel mass composition analysis based on archival data acquired from
1998 to 2013 provided by the KASCADE Cosmic ray Data Center (KCDC). The
analysis is based on modern machine learning techniques trained on simulation
data provided by KCDC. We present spectra for individual groups of primary
nuclei, the results of a search for anisotropies in the event arrival
directions taking mass composition into account, and search for gamma-ray
candidates in the PeV energy domain.

##  1 Introduction

Our knowledge of cosmic rays (CRs) remains sketchy even one century after
their discovery, see e.g. [1, 2]. In particular, the transition between the
dominance of Galactic and extragalactic sources that is expected to occur
somewhere between the CR “knee” and “ankle” is highly uncertain. There are
various approaches to unravel this mystery. Besides the indirect observation
via hadronic $\gamma$-ray and neutrino emission produced in CR interactions,
precision measurements of CR spectra can provide valuable information,
especially when their masses and charges are measured with high resolution.

The charge and mass reconstruction of high-energy CRs above the knee, which
are observable via extended air-showers, is a special problem since the
correlation of primary charge and mass with distribution of secondary
particles is weak. In addition, the reconstruction depends strongly on
hadronic interaction models introducing additional systematic uncertainties.
This is the main reason for skepticism regarding prior results obtained using
simplified methods and interpreted with, by now, outdated models, despite the
huge amount of analyzed data.

In these proceedings, we present the results of a novel analysis of archival
data collected by KASCADE [3] and provided by the KCDC service [4] based on
the latest hadronic interaction models and machine-learning techniques. The
reconstruction of the primary charge allows us to obtain CR spectra of
individual mass groups and to study CR anisotropies in terms of particle
rigidity. We also report the first steps towards the search for photons in
KASCADE data.

We will start in section 2 with a description of our methods based on a random
forest, which is an ensemble machine learning method, using sets of decision
trees on various sub-samples of training data to improve accuracy compared
with single decision tree. We present our results on mass composition in
Section 3 followed by an analysis of rigidity-dependent large-scale
anisotropies in Section 4. We present our results on the photon fraction in
Section 5 before we conclude in Section 6.

##  2 Methods and Data

We used KASCADE preselection data
sets111https://kcdc.iap.kit.edu/datashop/fulldata/, which contain the
following reconstructed air shower properties that we use for the training of
a classifier of primary mass: energy $E$; shower core coordinates $(x,y)$;
arrival direction $(\theta,\phi)$; muon and electron numbers
$\log_{10}N_{\mu}$, $\log_{10}N_{e}$; and shower age $s$. The KDCD service
provides CORSIKA [5] simulations with events generated for five individual
mass groups: H, He, C, Si, Fe. These simulations provide the same properties
as in real data reconstructed using the actual detector response. We have
trained a classifier to return one of the five mass groups based on three
modern hadronic interaction models: QGSJet-II.04 [6], EPOS-LHC [7] and Sibyll
2.3c [8].

Figure 1: Confusion matrix for our classifier trained using three different
hadronic interaction models, QGSJet-II.04 [6], EPOS-LHC [7] and Sibyll 2.3c
[8].

The confusion matrix of the classifier is shown in Fig. 1. One can see that
the matrix has a diagonal structure; however, non-diagonal elements are
heavily contaminated, which indicates systematic uncertainties. It is worth
noting, that these matrices obtained for the simulation data before
application of any quality cuts. In our work we used the quality cuts defined
by KASCADE collaboration [9], suggesting the following: $x^{2}+y^{2}<91$ m,
$\log_{10}N_{\mu}\geq 3.6$, $\log_{10}N_{e}\geq 4.8$, $0.2<s<2.1$,
$\theta<18^{\circ}$. Simulation study does not show any degradation of the
classifier performance, but reconstructed spectra indicate irregularities,
which might point to discrepancies between simulations and data beyond
official quality cuts.

For the study of CR anisotropies it is interesting to go for larger zenith
angles in order to cover a broader declination range. Figure 2 shows the
spectra of primary hydrogen (left panel) and carbon (right panel) mass groups
reconstructed with QGSJet-II.04 for zenith angles beyond $18^{\circ}$. The
spectra suggest that the zenith angle cut can be extended
$\mathcal{O}(30^{\circ})$. For this reason, we present spectra for a
conservative zenith angle cut of $\theta<18^{\circ}$, while we allow for a
somewhat looser cut of $\theta<30^{\circ}$ in our anisotropy study.

##  3 Cosmic-Ray Mass Composition

For the reconstruction of mass-dependent CR spectra we assume full efficiency
of the detector for events that pass the high-quality cuts described above.
The livetime of the detector is obtained from the fit of time differences
between two consecutive events with exponential decay function, i.e. assuming
a Poisson process. The full dataset corresponds to a livetime of
$T_{\mathrm{live}}\simeq 0.42\times 10^{9}$ s. Thus, the total exposure has
the form $\mathcal{E}\simeq\pi\sin^{2}\theta_{\rm max}\times S\times
T_{\mathrm{live}}$, where $\theta_{\rm max}=18^{\circ}$ is the maximum zenith
angle and $S=\pi(91{\rm m})^{2}$ is the surface area of the detector after
quality cuts. Figure 3 shows the CR spectra reconstructed using three
different hadronic interaction models. We show results for individual mass
groups and the total flux. These results are consistent with earlier findings
that reconstructions based on Sibyll show a trend towards heavier mass
compositions.

Figure 4 shows a comparison with recent results by IceCube/IceTop [10] based
on Sibyll 2.1. We have chosen this particular study for comparison, because it
provides results on four mass groups H, He, O, and Fe, which allows us to make
a more accurate and direct comparison of our results on H, He, C, and Fe. The
spectra of light components, H & He, are consistent between experiments within
uncertainties. The differences that we see between C and O and between the
individual spectra for Fe might be related to a “contamination” with Si, that
is not accounted for in the IceCube/IceTop analysis. A detailed study of this
effect is out of scope of this work. Another important point, which should be
kept in mind, is that we used default energy reconstructions from KCDC, which
is not corrected for the mass composition, so our spectra are not completely
unfolded.

Figure 2: Spectra of primary hydrogen (left) and carbon (right) mass groups
reconstructed with QGSJet-II.04 using zenith angles beyond KASCADE quality
cuts. The statistical uncertainties for heavier mass groups are too large to
allow for a study of systematic effects. The zenith bands are selected in
order to obtain equal exposure for each curve. The results indicate that the
zenith angle cut might be accurately pushed to $\mathcal{O}(30^{\circ})$,
thereby increasing the exposure by a factor $\simeq 3$.

Figure 3: Cosmic ray spectra for five individual mass groups and their sum
reconstructed from the full KCDC data (without spectral unfolding) using
different hadronic interaction models. We compare our results to those derived
by KASCADE [9] and IceCube/IceTop [10].

Figure 4: Comparison of spectra for four mass groups as provided by IceCube/IceTop [10] using Sibyll 2.1. It is important to point out that our analysis reconstructs the spectra of five different mass groups instead of only four used in the case IceCube/IceTop. The data corresponding to the spectrum of the Si group (not presented here) should be redistributed between the remaining four groups according to the classifier confusion matrix described above. This lack of events is clearly visible for the C/O and Fe groups, where contamination from Si is more significant than for the lighter components. $\mathcal{R}$ [${\rm V}$] | model | $N_{\rm tot}$ |  $A$ [$10^{-3}$] |  $\alpha$ [∘] |  $p$-value |  $A_{90}$ [$10^{-3}$]
---|---|---|---|---|---|---
$[10^{15.5},10^{16.0}]$ | EPOS LHC | $897,294$ | $10.1^{+5.5}_{-3.5}$ | $251\pm 28$ | $0.10$ | $17.1$
$>10^{16.0}$ | EPOS LHC | $79,140$ | $19.6^{18.0}_{-8.4}$ | $272\pm 48$ | $0.47$ | $44.3$
$[10^{15.5},10^{16.0}]$ | QGSJET-II.04 | $874,416$ | $14.3^{+5.4}_{-3.9}$ | $278\pm 20$ | $0.01$ | $21.1$
$>10^{16.0}$ | QGSJET-II.04 | $74,665$ | $18.7^{+18.5}_{-8.0}$ | $234\pm 51$ | $0.52$ | $44.3$
$[10^{15.5},10^{16.0}]$ | Sibyll 2.3c | $753,824$ | $7.7^{+5.9}_{-3.2}$ | $261\pm 40$ | $0.33$ | $15.6$
$>10^{16.0}$ | Sibyll 2.3c | $65,097$ | $14.3^{+20.5}_{-5.1}$ | $278\pm 67$ | $0.71$ | $42.7$
Table 1: Reconstructed dipole anisotropy using maximum-likelihood techniques
discussed in Refs. [11, 12]. Column 4 and 5 show the best-fit amplitude and
phase of the sidereal dipole anisotropy with 68% uncertainty range. The last
column shows the 90% C.L. upper limit on the amplitude.

Figure 5: Mollweide projections in equatorial coordinates of the reconstructed
anisotropy (top) and pre-trial significance (bottom) for the rigidity bin
$10^{15.5}<\mathcal{R}/{\rm V}<10^{16.0}$ based on QGSJET-II.04. We show the
results for a top-hat smoothing radius of $45^{\circ}$. The grey-shaded area
indicates the unobservable part of the celestial sphere. The dashed line
indicates the projection of the Galactic Plane. The values of pre-trial
significance are shown in units of standard deviations and indicated in red
and blue colors for excesses and deficits, respectively. The location of
maximum pre-trial significance is indicated by the symbol
$\boldsymbol{\times}$. The anisotropy reconstruction is based on a maximum-
likelihood method introduced in Ref. [13].

##  4 Rigidity-Dependent Anisotropy

The reconstruction of CR composition allows us – for the first time – to
analyze the anisotropy of CR arrival direction in terms of rigidity
$\mathcal{R}=pc/Ze$. Table 1 shows the results of the sidereal dipole
anisotropy using data with zenith angle $\theta\leq 30^{\circ}$. We bin the
data into two rigidity bins, $10^{15.5}{\rm V}<\mathcal{R}<10^{16.0}{\rm V}$
and $10^{16}{\rm V}<\mathcal{R}$, based on the average charge of the five
individual mass groups inferred by the composition analysis. The dipole
analysis is based on a maximum-likelihood method following Refs. [11, 12]. We
do not find strong evidence for large-scale anisotropies and place 90% C.L.
upper limits on the dipole amplitude (last column of Table 1).

The most significant excess with a $p$-value of $0.01$ is found for the first
rigidity bin using the composition based QGSJET-II.04. Figure 5 shows the
corresponding relative intensity (top panel) and pre-trial significance
(bottom panel) averaged over a radius of $45^{\circ}$. These all-sky results
are based on a maximum-likelihood method introduced in Ref. [13]. The smoothed
relative intensity is consistent with the best-fit orientation of the dipole
anisotropy. Note that the $45^{\circ}$ smoothing scale reduces the amplitude
of the excess compared to the dipole fit shown in Table 1. These results are
consistent with previous anisotropy measurements; see e.g. Ref. [14].

Figure 6: Performance of the photon classifier developed for KASCADE data.
Left: Distribution of classifier output $\zeta$ for primary hadrons and
photons. Right: Fraction of KASCADE events classified as photons as function
of $\zeta$.

##  5 Towards a Search for PeV Gamma-Rays

The search for high-energy $\gamma$-rays in the PeV domain is of special
interest. The absorption length of PeV $\gamma$-rays via pair production in
the cosmic microwave background is of the order of 10 kpc, comparable to the
distance of the solar system to the Galactic Center. Only Galactic PeVatrons
are visible via this channel, consistent with recent observations by HAWC
[15], Tibet-AS$\gamma$ [16] and LHAASO [17].

Since the expected fraction of PeV $\gamma$-rays are a few orders of magnitude
lower than the hadronic background, the binary classification approach will
not provide significant detection (the classifier provides only 1:10
suppression). In order to increase the efficiency of hadron separation we use
a random forest regressor returning a predicted floating point value
$\zeta\in[0,1]$ which can be treated as class membership probability of
reconstructed event. This allows us to choose a threshold for $\zeta$ that
optimizes the signal-to-background ratio (see left panel of Fig. 6). The
random forest consisting of 1000 trees gives us a suppression power in the
range $10^{2}$–$10^{3}$. When running the classifier on real data we see that
the number of $\gamma$-ray candidates corresponds to the suppression power
(right panel of Fig. 6), but the method requires further optimization.

##  6 Conclusion

We have presented the first results of a novel mass composition analysis based
on archival data of the KASCADE air shower experiment acquired from 1998 to
2013 and provided by the KASCADE Cosmic ray Data Center (KCDC). Using modern
machine learning techniques trained on data features provided by KCDC we have
obtained CR spectra for five mass groups represented by H, He, C, Si, Fe,
using latest hadronic models and machine learning algorithms. This allows us
to perform cross-checks with a state-of-the-art reconstruction recently
published by IceCube/IceTop [10]. For the first time, we performed a
reconstruction of large-scale anisotropy of CRs in the PeV energy domain as
function of rigidity. This intermediate success drives us to move towards
search for ultra-high energy photons in KASCADE data.

The results presented in these proceedings are only the first step in our
reanalysis of archival KASCADE data provided by KCDC. In future work, we plan
to use deep neural networks taking single station responses as input, which
are also provided by KCDC. Moreover, we plan to include KASCADE-Grande data in
order to push towards higher energies.

Last but not least, we were able to outreach our activity and participate in
JetBrains internship program222https://internship.jetbrains.com/projects/994/
and mathematical workshop organized by
NSU333https://bmm.mca.nsu.ru/project/33, thanks to the FAIR-
ness444https://www.go-fair.org/ of KCDC data. We are preparing software and
data release related for this analysis, some tutorial notebooks are already
avalable in Jupyter Hub at IAP KIT555https://jupyter.iap.kit.edu/.

## Acknowledgements

The authors would like to express gratitude to the colleagues from KCDC team.
M.A. acknowledges support from Villum Fonden under project no. 18994. The
development and testing of the classifier was supported by the state contract
with Institute of Thermophysics SB RAS.

## References

  * [1] S. Gabici, C. Evoli, D. Gaggero, P. Lipari, P. Mertsch, E. Orlando, A. Strong, and A. Vittino, Int. J. Mod. Phys. D 28 no. 15, (2019) 1930022, arXiv:1903.11584.
  * [2] R. Alves Batista et al., Front. Astron. Space Sci. 6 (2019) 23, arXiv:1903.06714.
  * [3] T. Antoni et al., (KASCADE Collaboration), Nucl. Instrum. Meth. A 513 (2003) 490–510.
  * [4] A. Haungs et al., Eur. Phys. J. C 78 no. 9, (2018) 741, arXiv:1806.05493.
  * [5] D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz, and T. Thouw, FZKA 6019 (1998) .
  * [6] S. Ostapchenko, Phys. Rev. D 83 (2011) 014018, arXiv:1010.1869.
  * [7] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner, Phys. Rev. C 92 no. 3, (2015) 034906, arXiv:1306.0121.
  * [8] F. Riehn, R. Engel, A. Fedynitch, T. K. Gaisser, and T. Stanev, EPJ Web Conf. 99 (2015) 12001, arXiv:1502.06353.
  * [9] T. Antoni et al., (KASCADE Collaboration), Astropart. Phys. 24 (2005) 1–25, arXiv:astro-ph/0505413.
  * [10] M. G. Aartsen et al., (IceCube Collaboration), Phys. Rev. D 100 no. 8, (2019) 082002, arXiv:1906.04317.
  * [11] M. Ahlers, Astrophys. J. 863 (2018) 146, arXiv:1805.08220.
  * [12] M. Ahlers, Astrophys. J. Lett. 886 no. 1, (2019) L18, arXiv:1909.09222.
  * [13] M. Ahlers, S. Y. BenZvi, P. Desiati, J. C. Díaz-Vélez, D. W. Fiorino, and S. Westerhoff, Astrophys. J. 823 no. 1, (2016) 10, arXiv:1601.07877.
  * [14] M. Ahlers and P. Mertsch, Prog. Part. Nucl. Phys. 94 (2017) 184–216, arXiv:1612.01873.
  * [15] A. U. Abeysekara et al., (HAWC Collaboration), Phys. Rev. Lett. 124 no. 2, (2020) 021102, arXiv:1909.08609.
  * [16] M. Amenomori et al., (Tibet ASgamma Collaboration), Phys. Rev. Lett. 126 no. 14, (2021) 141101, arXiv:2104.05181.
  * [17] Z. Cao et al., (LHAASO Collaboration), Nature 594 (2021) 33–36.