# Dynamic Collaborative Filtering for Matrix- and Tensor-based Recommender
Systems

Albert Saiapin  Delft University of Technology
Delft, The Netherlands
<EMAIL_ADDRESS>Ivan Oseledets  Skolkovo Institute of Science and
Technology
Moscow, Russia
<EMAIL_ADDRESS>Evgeny Frolov  Skolkovo Institute of Science and
Technology
Moscow, Russia
<EMAIL_ADDRESS>
###### Abstract

In production applications of recommender systems, a continuous data flow is
employed to update models in real-time. Many recommender models often require
complete retraining to adapt to new data. In this work, we introduce a novel
collaborative filtering model for sequential problems known as Tucker
Integrator Recommender - TIRecA. TIRecA efficiently updates its parameters
using only the new data segment, allowing incremental addition of new users
and items to the recommender system. To demonstrate the effectiveness of the
proposed model, we conducted experiments on four publicly available datasets:
MovieLens 20M, Amazon Beauty, Amazon Toys and Games, and Steam. Our comparison
with general matrix and tensor-based baselines in terms of prediction quality
and computational time reveals that TIRecA achieves comparable quality to the
baseline methods, while being 10-20 times faster in training time.

###### Index Terms:

Sequence-aware tensor factorization, dynamical modeling, projector-splitting
integrator, collaborative filtering, next-item prediction.

##  I Introduction

Recommender systems have become integral technologies across a multitude of
applications, significantly enhancing user experiences and delivering
substantial business value, as exemplified by industry leaders such as Amazon
and YouTube.

Yet, a fundamental challenge persists in this dynamic landscape: How can we
effectively refine model parameters within the continuous flow of data?
Envision the life cycle of a production-level recommender system—a process
involving two critical stages: using the latest model for predicting user
preferences and the rigorous task of retraining or updating the model to adapt
to the ever-evolving data flow.

At the core of this challenge is the optimization of model retraining. It’s a
quest to adjust model parameters seamlessly, aligning with the data dynamics,
all while preserving the patterns within the structure of model weights. This
is a routine procedure every recommender model performs to be up to date.
However, the latest state-of-the-art models [1, 2, 3] cannot deal with the
problem of adding new users or items to a recommender system because the whole
neural network architecture changes.

In this paper, we choose a different path. We initiate our exploration by
addressing the problem formulation of the PureSVD model [4] and progressively
refine it to enable the effective updating of user and item embeddings in
response to incoming data.

The main idea we start from is to incorporate the notion of time into the
optimization problem formulation and then use the same techniques as the
authors of [5] applied. However, we extend their model by processing new users
and items, and generating embeddings for them in an incremental fashion.

In addition to it, as was mentioned above, sequential models based on Neural
Networks show state of the art quality on the task of the next item prediction
[1, 6, 2]. However, as the authors of [7] showed, properly prepared
factorization models can show the same or superior quality. With this in mind,
we propose a model that can use sequential information about user preferences
based on GA-SATF [7] model. It requires a more complex model formulation
compared PureSVD [4]. Namely, we are using a three-dimensional tensor data
representation as there are three different modes in the data: a user, an
item, and a sequential order.

The main contributions 111https://github.com/AlbMLpy/DynamicCF of this work
can be summarized as follows:

  1. 1.

We develop a matrix-based dynamic collaborative filtering recommender system
using Projector Splitting Integrator (PSI). The model can process new users
and items on the fly.

  2. 2.

We develop a novel tensor-based sequential collaborative filtering recommender
system using Tucker Integrator (TI). The model can process new users and items
dynamically.

  3. 3.

We suggest and experimentally check several methodologies to initialize and
further update embeddings (model factors) for new users and new items.

As a result, the developed models could potentially replace many models in
production, from PureSVD [4] to SASRec [1]. They could be faster than state-
of-the-art neural networks, save money, reduce time consumption, and preserve
the same or comparable accuracy. On the other hand, these models could
emphasize the importance of developing general-purpose mathematical frameworks
and continuing the research in the dynamical modeling and tensor factorization
areas as an alternative to mainstream research directions based on neural
networks.

##  II Related Work

###  II-A Recommender Systems

There are different models in recommender systems field. Some of them rely on
content information [8] while others on collaborative signals [9]. In our
study we will focus on collaborative filtering for implicit data that can be
sequential [10]. Moreover, the authors of [4] suggest to use top-n problem
formulation where the goal of the recommender system is to find a few specific
items which are supposed to be most relevant to the user. As this problem
formulation better reflects commercial applications we would rely on it in
this work.

###  II-B Low-Rank Approximations

This section reviews the research related to matrix and tensor decompositions
which are the basis of this work. The authors of [11] give a comprehensive
overview of general tensor decompositions and mathematical concepts behind
them. Nevertheless, in order to handle big data and compute tensor
decomposition out of it one should use different randomization techniques. One
can find an extensive review of recent advances in randomization for the
computation of a Tucker decomposition and Higher Order SVD in [12]. The
connection between recommender systems field and tensor decompositions can be
found in [13].

However, in order to solve the dynamical problem formulation, which will be
presented later, one have to explore the notion of dynamical low-rank
approximations. Several studies [14, 15] present an increment-based
computational approach for the low-rank approximation of time-dependent data
matrices and solutions to matrix differential equations. It is shown that the
approach gives locally quasi-optimal low-rank approximations. A recent study
[5] concluded that it is possible to integrate dynamical matrix factorization
framework into a recommender system in order to get more stable solution.

The following studies [16, 17] present high-dimensional problem formulation.
The authors try to solve the same problem of the approximation of time-
dependent data tensors and of solutions to tensor differential equations by
tensors of low Tucker rank. This approach works with the increments rather
than the full tensor. In this method, the derivative is projected onto the
tangent space of the manifold of tensors represented in Tucker format.

Another area of low-rank approximations research is connected with incremental
parameters updates. For example, the question of how an SVD may be updated by
adding new rows and/or columns of data, which may be missing values or
contaminated with correlated noise is considered in [18, 19]. The size of the
data matrix need not be known as an SVD is developed as the new data comes in.
The paper [20], describes an online tensor subspace learning algorithm which
models appearance changes by incrementally learning a tensor subspace
representation through adaptively updating the sample mean and an eigenbasis
for each unfolding matrix of the tensor.

###  II-C Sequential Models

Another substantial part of this research is connected with sequential
recommender systems. For example, in [1] the authors proposed a self-attention
based sequential model (SASRec) that allows to capture long-term semantics
(like an RNN), but, using an attention mechanism, makes its predictions depend
on relatively few actions (like a Markov Chain). At each time step, SASRec
seeks to identify which items are ‘relevant’ from a user’s action history, and
use them to predict the next item. Top-n sequential recommendation models
represent each user as a sequence of items he interacted in the past and aims
to predict top-n ranked items that a user will likely interact in a “near
future”. The order of interaction implies that sequential patterns play an
important role where more recent items in a sequence have a larger impact on
the next item.

In order to improve sequential recommender system performance the authors of
[6] suggest to use the internal correlations in data to improve the data
representations via pre-training methods. They consider several self-
supervised objectives to learn the correlations among attribute, item,
subsequence, and sequence by using the mutual information maximization (MIM)
principle which is particularly suitable in their scenario. The experimental
results show better prediction quality results for S3Rec than for all the
baselines including SASRec. One more example of a neural sequential
recommender could be - DuoRec [2]. This model uses a contrastive objective to
address the representation degeneration problem by putting regularization over
sequence representations. This model shows superior performance to SASRec and
S3Rec models.

Previous studies have almost exclusively focused on neural models to model
sequential behaviour of a user with competitive prediction quality. A recent
study [7] develops the model that mimics transformer self-attention with an
alternative and more lightweight approach. It is a tensor factorization-based
model that ingrains the structural knowledge about sequential data within the
learning process. The authors demonstrate how certain properties of a self-
attention network can be reproduced using special Hankel matrix
representation. The resulting model has a shallow linear architecture and
compares competitively to its neural counterpart.

###  II-D Literature Review Conclusions

In this work we concentrate on the top-n problem formulation for sequential
data with implicit feedback. It is more production oriented setting that can
reflect the real performance of the recommender system. The main model we
develop will be based on several distinct components:

  * •

Tensor model for sequential recommendations - GA-SATF - based on [7].

  * •

Computational core of matrix and tensor decompositions - based on [12].

  * •

Dynamical modelling integrators - based on [15, 21]

  * •

Addition of new users and items - based on [18]

##  III Problem Statement

The primary objective of a general recommender system is to leverage prior
information about users, items, and additional contextual knowledge to
accurately predict the most relevant items for a selected user. The relevance
of items is typically measured using a relevance score function denoted as
$f_{R}$, which is estimated from implicit feedback provided by users.

This study is an attempt to address the issue of sequence-aware tensor
factorization for the next item prediction. The data utilized in this research
can be organized in the form of a third order tensor $\mathcal{X}$ such that:

$x_{ijk}=\begin{cases}1,&\text{if j-th item is in position $p^{i}_{k}$ of i-th
user}\\\ 0,&\text{otherwise},\end{cases}$ (1)

where $k=p_{k}^{i}-n_{i}+L$ and $n_{i}$ is the total number of items in the
ordered history of user i. To ensure consistency, any user history with a
length exceeding $L$ is truncated, such that $n_{i}\leq L$ for all users. Note
that, regardless of the length of the original sequence, the most recent item
is always located at position $L$. If a sequence contains fewer that $L$
items, zero padding is employed to fill in the remaining positions up to $L$.

The aforementioned data representation allows to define a tensor model using
relevance score function:

$f_{R}:\text{User}\times\text{Item}\times\text{Position}\rightarrow\text{Relevance}$
(2)

that assigns some relevance score to each triplet of the observed interactions
between users and items considering their respective positions in the user’s
history of interactions. As user consumption behavior may exhibit variability
and the history length could potentially approach the size of the entire item
catalog, computational challenges may arise due to the high dimensionality of
positional encoding. However, it is reasonable to suppose that only a
relatively small number of the most recent items significantly contribute to
explaining current user decisions. Hence, truncating user sequences to a
length of $L$, which is considerably smaller than the catalog size $N$, is a
pragmatic approach.

In order to represent the tensor model (2) we will use the mathematics of Low-
Rank Tensor Decomposition [11]. Apart from that, we suppose that the incoming
tensor data depends on time. That is why, to encode the dynamical nature of
the system we will use the theory from [21, 16]. We would consider the
approximation of time-varying tensors by tensors of low Tucker rank. The
approach to get tensor decomposition factors is a continuous-time updating
technique, which works only with the increments in the tensors rather than the
tensors themselves and does not require computing any decompositions of large
tensors.

Consider a time-varying family of tensors
$\mathcal{X}(t)\in\mathbb{R}^{M\times N\times L}$ for $0\leq t\leq\hat{t}$.
Let $\mathcal{M}_{r}$ denote the manifold of all order-3 tensors of Tucker
rank r = $(r_{1},r_{2},r_{3})$, where $r_{1}\leq M$, $r_{2}\leq N$, $r_{3}\leq
L$. The best approximation to $\mathcal{X}(t)$ in $\mathcal{M}_{r}$ with
respect to Frobenius norm is:

$\|\mathcal{X}(t)-\mathcal{R}(t)\|\rightarrow\min_{\mathcal{R}(t)\in\mathcal{M}_{r}}$
(3)

where
$\mathcal{R}(t)=\mathcal{C}(t)\times_{1}U(t)\times_{2}V(t)\times_{3}W(t)$ \-
Tucker decomposition [11], the factor matrices $U(t)\in\mathbb{R}^{M\times
r_{1}}$, $V(t)\in\mathbb{R}^{N\times r_{2}}$, $W(t)\in\mathbb{R}^{L\times
r_{3}}$ represent user, item, positional embeddings accordingly.

However, this problem formulation does not take into account the connection
between two subsequent steps $t_{1}$ and $t_{2}$. That is why, the authors of
[21] consider instead the dynamical tensor approximation
$\mathcal{Y}(t)\in\mathcal{M}_{r}$ determined from the condition that for
every time step - t the derivative $\dot{\mathcal{Y}}(t)$ which lies in the
tangent space $\mathcal{T}_{\mathcal{Y}(t)}\mathcal{M}_{r}$ can be chosen as:

$\begin{dcases}\|\dot{\mathcal{Y}}(t)-\dot{\mathcal{X}}(t)\|\rightarrow\min_{\dot{\mathcal{Y}}(t)\in\mathcal{T}_{\mathcal{Y}(t)}\mathcal{M}_{r}}&\\\
\mathcal{Y}(0)=\mathcal{R}(0)&\end{dcases}$ (4)

Moreover, in this work we suppose that the number of users $M$ and items $N$
also depends on time - $t$ such that:

$\begin{cases}M=M(t)\text{, }M(t_{1})\leq M(t_{2})\text{, }t_{1}\leq t_{2}&\\\
N=N(t)\text{, }N(t_{1})\leq N(t_{2})\text{, }t_{1}\leq t_{2}&\end{cases}$ (5)

In this regard, we can consider four distinct cases of user and item
interactions:

  1. 1.

Interactions with new users

  2. 2.

Interactions with new items

  3. 3.

Interactions with known users and known items

  4. 4.

Interactions with new users and new items

The overall purpose of the study is to solve the problem (4), (5) that will
allow to develop a tensor-based recommender system that can update user and
items embeddings effectively. The developed model is supposed to work with
data chunks rather than all the data up to the moment and should be able to
add embeddings for new users and items.

##  IV Dynamic Collaborative Filtering

Figure 1: General integrator-based recommender system framework. New entities
mean new users that interacted only with new items. New users represent a set
of new users that interacted with known items. New items contain a set of new
items that were consumed by the known users. New interactions represent
interactions of known entities.

Let us introduce the methodology to solve the problem (4), (5) presented on
Figure 1. Firstly, we have some initial user-item interactions data. Using it
we can find user and item embeddings using SVD or HOOI depending on
dimensionality - which are model parameters. With these parameters the model
can generate recommendations for a user. Secondly, the recommender system is
dynamic because in a day or a week the new data can be generated (e.g. new
users register on a platform, people create more content, people consume new
content). That is why at some point the recommender model should be trained
with the incoming data. In the framework of matrix and tensor decompositions
one would have to either retrain the model from scratch or use time window to
create a new model. In the framework of integrators we suggest to use the four
step procedure depicted on Figure 1.

The next section describes technical details needed for implementation of this
framework in matrix and tensor cases: how to train initial factors; how to add
new entities, users, items; how to add new interactions; how to do the ranking
of items for a particular user.

###  IV-A Matrix Case

First of all, in two-dimensional case one would find the initial user $U$ and
item $V$ embeddings using an SVD of incoming data matrix [4].

Algorithm 1 PSIRec: Update known user/item factors

1:

2:\- $\Delta A$ \- new data collected since the last time-step for known
users/items

3:\- $U_{0}$, $S_{0}$, $V_{0}^{\top}$ \- last state of the model at time $t-1$

4:$U_{1}$, $S_{1}$, $V_{1}^{\top}$ \- new state of the model at time $t$

5:$U_{1}$, $S^{\prime}$ = QR($U_{0}S_{0}+\Delta AV_{0}$)

6:$L=V_{0}(S^{\prime}-U_{1}^{\top}\Delta AV_{0})^{\top}+\Delta A^{\top}U_{1}$

7:$V_{1}$, $S_{1}^{\top}$ = QR($L$)

Secondly, in this research we present several techniques to incrementally add
embeddings for new users and/or new items (step 2-3 on Figure 1). The first
proposition is to use zero embeddings as initialization for further update by
PSI - Algorithm 1. Zero initialization of new factors does conflict with the
orthogonality condition for the factor matrices - $U$, $V$. The second
proposition is to use random initialization for new user/item embeddings. As
the addition of new rows - $\hat{U}$, $\hat{V}$ to matrices - $U$, $V$ would
break the orthogonality condition one have to apply QR decomposition to
recover the structure (example for user factors):

$\begin{dcases}\begin{bmatrix}U\\\ \hat{U}\end{bmatrix}=QR\\\ U\leftarrow Q\\\
S\leftarrow RS\end{dcases}$ (6)

These new embeddings are further updated by Algorithm 1 on the forth step. The
third proposition is to apply the idea of incremental SVD developed by [18].
You can find the implementation details in Algorithm 5. To add new users that
interacted with known items one has to swap $U$ and $V$ in Algorithm 5 and
swap resulting factors accordingly.

Another key thing to remember is that there might be a situation when a new
user that is not represented in a recommender system interacted with a new
item. In this case, user and item embeddings can be updated using an SVD of a
block-diagonal matrix represented in Algorithm 2 which is the first step in
update factors procedure on Figure 1.

Algorithm 2 PSIRec: Add new users and new items

1:

2:\- $\Delta A$ \- data representing new users and new items interactions

3:\- $U_{0}$, $S_{0}$, $V_{0}^{\top}$ \- state of the model

4:$U_{3}$, $S_{3}$, $V_{3}^{\top}$ \- updated state of the model

5:$U_{1}$, $S_{1}$, $V_{1}^{\top}$ = SVD ($\Delta A$)

6:$U_{2},S_{2},V_{2}^{\top}$ = SVD $\left(\begin{bmatrix}S_{0}&0\\\
0&S_{1}\end{bmatrix}\right)$

7:$U_{3}\leftarrow\begin{bmatrix}U_{0}&0\\\ 0&I\end{bmatrix}U_{2}$

8:$S_{3}\leftarrow S_{2}$

9:$V_{3}\leftarrow\begin{bmatrix}V_{0}&0\\\ 0&I\end{bmatrix}V_{2}$

The final step in update factors procedure is to add new interactions. This
step uses Projector-Splitting Integrator (PSI) [15] technique. PSI offers an
incremental scheme that takes only the newly added data to directly update the
model, which eliminates the need to recompute the entire model from scratch.
See Algorithm 1 for implementation details.

Speaking of model inference, it is the same as for a PureSVD model [4]. Hence,
given some user preferences binary vector $p\in\mathbb{R}^{M(t)}$, the list of
top-n recommendations can be generated as:

$\text{toprec}(p,n)=\operatorname*{arg\,max}^{n}VV^{\top}p,$ (7)

where $V$ \- item embeddings matrix, $\text{toprec}(p,n)$ \- returns $n$ item
indices of maximum relevance score. The presented model is called as PSIRec.

###  IV-B Tensor Case

According to GA-SATF model from [7] to mimic the positional attention
mechanism and make the model be able to learn from sequential data we would
use the following data in this section:

$\mathcal{X}(t)\equiv\mathcal{X}(t)\times_{3}A^{\top},$ (8)

where $A\in\mathbb{R}^{L\times L}$ \- attention matrix such that:

$A=\begin{bmatrix}a_{1}&0&\dots&0\\\ a_{2}&a_{1}&\dots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ a_{L}&\dots&a_{2}&a_{1}\end{bmatrix}$ (9)

$a_{L}=L^{-f}$, $f\geq 0$ \- scaling hyper-parameter.

To begin with, to get the initial embeddings for users and items we apply HOOI
[12] to three dimensional sequential data tensor (8).

Secondly, to add new users and/or new items (step 2-3 on Figure 1) to a model
one can apply different strategies which are the same as for a matrix case
IV-A. However, the third proposition should be extended on tensor case which
is presented in Algorithm 6.

Thirdly, in order to add new users that interacted with new items and vise
versa we have to extend Algorithm 2 on tensor case. The incoming data can be
naturally represented as a block tensor and the update procedure is depicted
in Algorithm 3.

Algorithm 3 TIRec: Add new users and new items

1:

2:\- $\Delta\mathcal{A}$ \- data representing new users and new items ordered
interactions

3:\- $U_{1}^{0}$, $U_{2}^{0}$, $U_{3}^{0}$, $\mathcal{C}^{0}$ \- state of the
model

4:$U_{1}^{1},U_{2}^{1},U_{3}^{1},\mathcal{C}^{1}$ \- updated state of the
model

5:$U_{1}$, $V_{1}$, $\mathcal{C}_{1}$ = Tucker2
($\Delta\mathcal{A}\times_{3}(U_{3}^{0})^{\top}$)

6:$\hat{U},\hat{V},\hat{W},\mathcal{C}^{1}$ = HOSVD
$\left(\begin{bmatrix}\mathcal{C}^{0}&0\\\
0&\mathcal{C}_{1}\end{bmatrix}\right)$

7:$U_{1}^{1}\leftarrow\begin{bmatrix}U_{1}^{0}&0\\\
0&U_{1}\end{bmatrix}\hat{U}$

8:$U_{3}^{1}\leftarrow U_{3}^{0}\hat{W}$

9:$U_{2}^{1}\leftarrow\begin{bmatrix}U_{2}^{0}&0\\\
0&V_{1}\end{bmatrix}\hat{V}$

The forth step is to add new interactions in update factors procedure on
Figure 1. This step uses Tucker Integrator (TI) [21] to dynamically update
user, item, positional embeddings. See Algorithm 4 for implementation details.

Speaking of model inference, it is the same as for a GA-SATF model [7]. Hence,
given some user preferences binary matrix $P\in\mathbb{R}^{M(t)\times L}$, the
list of top-n recommendations can be generated as:

$\text{toprec}(P,n)=\operatorname*{arg\,max}^{n}VV^{\top}PSAW\hat{w}_{L}$ (10)

where $V$ \- item embeddings matrix, $W$ \- position (order) embedding matrix,
$A$ \- attention matrix, $\text{toprec}(P,n)$ \- returns $n$ item indices of
maximum relevance score, $\hat{w}=[A^{-\top}W]_{L}$.
$S=[\delta_{k,k^{\prime}+1}]\in\mathbb{R}^{L\times L}$ \- lower shift matrix
that decreases positions of all observed items in $P$ by one. Note that if
length of a user sequence is exactly $L$, the first item in the sequence gets
discarded, which satisfies the length-$L$ requirement for user histories in
tensor construction. The presented model that uses the third proposition (see
IV-A) is called as TIRec. The model that uses the first proposition is called
TIRecA.

Algorithm 4 TIRec: Update known user/item factors

1:

2:\- $\Delta\mathcal{A}$ \- new data collected since the last time-step for
known users/items

3:\- $U_{1}^{0}$, $U_{2}^{0}$, $U_{3}^{0}$, $\mathcal{C}^{0}$ \- state of the
model

4:$U_{1}^{1},U_{2}^{1},U_{3}^{1},\mathcal{C}^{1}$ \- updated state of the
model

5:for $i\in(1,2,3)$ do

6: $Q_{i}^{0}$,
$(S_{i}^{0})^{\top}=\text{QR}\left((\mathcal{C}^{0}_{[i]})^{\top}\right)$

7:
$(V_{i}^{0})^{T}=\left((Q_{i}^{0})^{\top}_{[i]}\bigtimes_{k=1}^{i-1}(U_{k}^{1})^{\top}\bigtimes_{l=i+1}^{3}(U_{l}^{0})^{\top}\right)_{[i]}$

8: $K_{i}^{1}=U_{i}^{0}S_{i}^{0}+\Delta\mathcal{A}_{[i]}V_{i}^{0}$

9: $U_{i}^{1}$, $S_{i}^{1}=\text{QR}\left(K_{i}^{1}\right)$

10: $S_{i}^{0}=S_{i}^{1}$

11: $S_{i}^{1}\leftarrow
S_{i}^{0}-(U_{i}^{1})^{\top}\Delta\mathcal{A}_{[i]}V_{i}^{0}$

12: $\mathcal{C}^{0}=(S_{i}^{1}(Q_{i}^{0})^{\top})_{[i]}$

13:end for

14:$\mathcal{C}^{1}=\mathcal{C}^{0}+\Delta\mathcal{A}\bigtimes_{i=1}^{3}(U_{i}^{1})^{\top}$

##  V Experiments

In this section we describe the general evaluation setup, datasets and
preprocessing steps for performing experiments.

###  V-A Data

In this work, we use the following publicly available datasets in recommender
systems field: Movielens-20M (ML-20M), Amazon Beauty (AMZ-B), Amazon Toys and
Games (AMZ-G), and Steam [1]. We follow the same data initial preparation
procedure as in [7]. Namely, we leave no less than five interactions per each
user and each item for all datasets except Movielens-20M. The explicit values
of ratings are transformed into an implicit binary signal indicating the
presence of a rating. As in the Steam data some users assigned more than one
review to the same items, such duplicates were removed. In Movielens-20M data
we use 20% of the last interactions. The resulting statistics for the datasets
are provided in Table I.

###  V-B Metrics

In order to compare the prediction quality of different methods we would use
the following metrics: HR@n, MRR@n, WJI(u, v)@n defined below:

$\text{HR}@\text{n}=\dfrac{\lvert U_{\text{hit}}^{n}\rvert}{\lvert
U_{\text{all}}\rvert}$ (11)

where $\lvert U_{\text{hit}}^{n}\rvert$ is the number of users for which the
correct answer is included in the top-n recommendation list, $\lvert
U_{\text{all}}\rvert$ is the total number of users in the test set.

$\text{MRR}@\text{n}=\dfrac{1}{\lvert U_{\text{all}\rvert}}\sum_{u=1}^{\lvert
U_{\text{all}}\rvert}\dfrac{\text{relevance}(u)}{\text{rank}(u)}$ (12)

$\text{relevance}(u)$ is 1 if a test item is included in the top-n
recommendation list else 0, $\text{rank}(u)$ is a position of a test item in
the top-n recommendation list, $\lvert U_{\text{all}}\rvert$ is the total
number of users in the test set.

In order to measure the stability of recommendations we would use the method
from [5]. The top-n similarity between two recommendations lists is measured
as the Weighted Jaccard Index (WJI) of their ”bag-of-items” (BOI)
representations $u$, $v$:

$\text{WJI}(u,v)@\text{n}=\dfrac{\sum_{j=1}^{N}\min(w_{j}(u),w_{j}(v))}{\sum_{j=1}^{N}\max(w_{j}(u),w_{j}(v))}$
(13)

where $u=(u_{i})_{i=1}^{N}$ and $u_{i}$ denotes the position of item $i$ in
the list of recommendations; $N$ \- the number of items;
$w_{i}(u)=\dfrac{1}{u_{i}}$ and $w_{i}(u)=0$ if $u_{i}>n$.

###  V-C Evaluation

As in a general Machine Learning field, we preprocess the data and split it
into train, validation and test subsets. Validation is used for hyper-
parameters optimization, whereas test is used for final evaluation. In our
settings test data is divided into chunks - chunk1, chunk2, …, etc
sequentially day by day. You can find data dynamic on Figures 2, 3 for ML-20M
dataset.

Figure 2: ML-20M. How many users/items are processed by a recommender system.
#user - the number of new users; #items - the number of new items;
#new_users_ni - the number of new users that interacted with new items;
#new_items_nu - the number of new items that were consumed only by new users.

The first step is to merge the validation data back into the train data when
the optimal hyper-parameters are found. Then, we can retrain the model using
train data and optimal hyper-parameters. After that, we can evaluate the model
on chunk1. To do that, we have to find a set of users that are present in the
train data and chunk1. The first items these users interacted with in chunk1
are the target we are going to compare against the recommended lists of items
using metrics defined above in V-B. Then we can retrain the model using all
the data up to the moment - train and chunk1 or update the model parameters
using only chunk1 which depends on the model. The same procedure is applied
several times to different chunks sequentially. In the end one can get the
dynamical behavior of the model based on measured metrics.

In order to evaluate a recommender system’s stability we find a set of target
users beforehand. They must be present in train data and in many chunks.
Top-50 users are chosen for that purpose. To measure stability we compare
ranking lists from two subsequent time steps for target users and take an
average by their number.

For Amazon Beauty and Amazon Toys and Games the following settings were used:
70% of initial data goes into train; 0.01% out of train is kept for
validation; 100 chunks (days) were used for dynamical evaluation; maximum
length of a user sequential history - $K=20$. Experimental setup for
Movielens-20M data is: 40% of initial data goes into train; 0.01% out of train
is kept for validation; 200 chunks (days) were used for dynamical evaluation;
maximum length of a user sequential history - $K=20$. For experiments with the
final Steam dataset we used the settings: 80% of initial data goes into train;
0.01% out of train is kept for validation; 50 chunks (days) were used for
dynamical evaluation; maximum length of a user sequential history - $K=8$;

Figure 3: ML-20M. Dynamical nature of data in recommender systems field.

###  V-D Baselines

The main idea of this work is to develop faster alternatives to SVD-based
models and models based on Tucker Decomposition. The first baseline is PureSVD
model (SVD) [4]. It is a strong baseline in recommender systems field. In the
settings of this work, this model is retrained with all the data up to date.
The second baseline is TDRec. This model applies general Tucker Decomposition
to incoming data tensor and finds factors and core out of it. This model is
also retrained from scratch every time step using all the gathered data. The
final baseline is TDRecReinit. This model applies Tucker Decomposition to
incoming data as well. But on the next time step it uses the factors from the
previous iteration that allows for faster and more robust retraining.

###  V-E Hyper-Parameters Grid Search

In this work we compare different models that can be categorized into two
types: matrix-based and tensor-based.

Matrix-based models are PureSVD and PSIRec and they have 1 hyper-parameter:
rank r. For rank values r, their range is (10, …, 300) with step size 10 for
all datasets but Steam. For Steam data, rank values range is (10, …, 400) with
step size 10.

Tensor-based models are TDRec, TDRecReinit, TIRec, TIRecA. These models have 2
main hyper-parameters: rank r = $(r_{1},r_{2},r_{3})$ and scaling parameter
$f$ (see IV-B). For rank values $r_{1}$, $r_{2}$, their range is (32, 64, 100,
128, 256), $r_{3}$ range is (5, 10), $f\in(0,2,4)$ for all datasets but Steam.
For Steam data, rank values $r_{1}\in(64,128,256)$, $r_{2}\in(64,100)$,
$r_{3}\in(2,4)$ and scaling factor $f\in(0,2)$.

In order to find optimal hyper-parameters we fix parameters configuration,
then train the model using training data and measure prediction quality using
validation data. After that we can find the best parameters configuration
based on metric results. The source code to reproduce our work is openly
published online 222https://github.com/AlbMLpy/DynamicCF.

##  VI Results and Discussion

This section presents experimental comparison of various models in terms of
computational time, prediction quality and stability of recommendations. Apart
from that we compare different regimes to add embeddings for new users and/or
items in the proposed frameworks. We present graphical results only for ML-20M
dataset due to space restrictions. The average results for all other data can
be find in the following tables: Table II shows comparison of proposed models
with the baselines V-D; Table IV presents the metrics results for different
options to add new embeddings to a model based on Tucker Integrator; Table III
presents the metrics results for different options to add new embeddings to a
model based on PSI.

###  VI-A Computational Time Results

Figure 4: ML-20M. Computational time dynamics

Let us start from the comparison of the training time per iteration across
different datasets and models. As we can conclude from the Figure 4 the
slowest model in terms of computational time is TIRec. The main problem of
TIRec is the addition of new users or new items Algorithm 6. Being
mathematically correct, it is a bottleneck of the whole method which prevents
scalability. One have to calculate an SVD of a wide dense matrix (see line 5
in Algorithm 6). The second to worst model is TDRec. It is expected as the
model is retrained from scratch every time step with all the data up to the
moment. It is worth mentioning that the fluctuations in the green curve could
be a result of different number of iterations in the inner cycle of TDRec
training procedure (see HOOI [11]). TDRecReinit baseline needs much less
iterations to converge then TDRec because of proper initialization and it
allows to accelerate the model. The results for TIRecA model confirm that
removing Algorithm 6 from TIRec allows to get increase in training speed. This
model is about 25 times faster then TIRec, 9 times faster then TDRecReinit, 20
times faster then TDRec and even 3 times faster then SVD (see Table II). The
best model according to the plot is PSIRec. The results for the other three
datasets, AMZ-B and AMZ-G, Steam are similar to the ML-20M data and can be
found in an accumulated way in Table II.

###  VI-B Prediction Quality Results

Figure 5: ML-20M. Dynamics of HR metric (Cumulative)

Let us further analyse the results of HR metric dynamics for ML-20M dataset on
Figure 5. As one can see the sequential information is important in this
particular dataset as the difference between matrix and tensor models is
significant. In its turn, PSIRec model results are very similar to SVD
baseline. The results for tensor-based models are also close to each other
with best model - TDRecReinit. This graph shows the behavior of integrator-
based models that we wanted to achieve - the prediction quality is close to
the baselines. According to Table II, one can conclude that the sequential
information is not that important for Steam and AMZ-G datasets as PSIRec and
SVD show better HR@5 metric results than tensor counterparts. The best tensor
model on these data is TDRecReinit. For AMZ-B dataset one can see a major
difference in prediction quality for tensor methods: 0.0013 for TIRecA model
against 0.020 in TDRecReinit. This model behavior can be possibly explained by
a significant change in users’ behavior which breaks the assumptions (too
large step size - h) under the Tucker Integrator [21].

###  VI-C Stability Results

Figure 6: ML-20M. Stability of recommendations

This section presents the findings on the stability of different recommender
models in dynamical settings. Consider the example of the ML-20M dataset
illustrated in Figure 6. The WJI@5 value indicates the stability of the model,
where a larger value implies less variation in recommendations between
subsequent time steps. Matrix models exhibit the best stability performance,
with an average value of 0.668 for PSIRec and 0.649 for SVD (refer to Table
II). TIRecA performs better than the other tensor models, surpassing them by
approximately twofold. TDRec is the least stable model, with instances where
the WJI metric reaches 0, indicating a completely different recommendation
list compared to the previous time step. That is not a stable recommendation.
The result on stability for the other datasets can be found in Table II). We
can conclude that the most stable matrix-based model is PSIRec and tensor-
based one is TIRecA most of the time.

###  VI-D Update Methods Results

This section presents the results of different initialization strategies for
new user/item factors in integrator-based framework (see IV-A).

Let’s start with PSIRec model on Figure 7. As we can see random initialization
does not improve the prediction quality of a model. The best one is pure
PSIRec using Algorithms 5, 2 to add new embeddings to a system. The same
behavior holds true for all the datasets (see Table III). It should be noted
that, the fastest models are PSIRec and PSIRec(0;0) as they do not use extra
QR decomposition. That is why, these methods are preferred as fast and
accurate.

Now let’s consider TIRec model results depicted on Figure 8 and Table IV. The
best model in terms of HR metric quality is TIRec on ML-20M data. However, on
Amazon data the models based on random initialization show more accurate
results - TIRecA(0.001;B) and TIRecA(1e-05;0). Speaking of training time, the
best configuration is TIRecA(0;0) as it uses pure Tucker Integrator to update
the user/item factors.

Figure 7: ML-20M. Comparison of different initialization strategies for PSIRec
model. PSIRec($\sigma$;$Y$) depicts different model configurations, where
$\sigma$ taken from $\mathcal{N}(0,\,\sigma^{2})$ initialization of factors;
$Y=0$ \- zero initialization in the first step of Figure 1, $Y=\text{B}$ \-
use Algorithm 2. The results are represented using a sliding window with step
1 and window length 20. Figure 8: ML-20M. Comparison of different
initialization strategies for TIRec model. TIRecA($\sigma$;$Y$) depicts
different model configurations, where $\sigma$ taken from
$\mathcal{N}(0,\,\sigma^{2})$ initialization of factors; $Y=0$ \- zero
initialization in the first step of Figure 1, $Y=\text{B}$ \- use Algorithm 3.
The results are represented using a sliding window with step 1 and window
length 20.

##  VII Conclusion

The main purpose of this work was to design and implement a novel sequential
recommender system in the framework of tensor data representation that would
be able to update its parameters effectively using incoming data and would
save appropriate prediction quality. This model was developed and it is called
- TIRecA. Comprehensive set of experimental results confirms that this model
can be 10-20 times faster than the general models based of Tucker
decomposition. Speaking of metric results, the model showed comparable
performance to all the baselines under different settings. However, it should
be noted that the model could degrade in the long run and in order to restore
the performance one have to retrain the model from scratch periodically.

Aforementioned results provide a good starting point for discussion and
further research in the area of sequential recommender systems. One possible
future research topic could be combination of the proposed TIRecA model and
model based on Neural architecture. The idea is to use the training speed of
tensor model and combine it with prediction quality of contemporary neural
models. Another interesting topic would be an adaptation of TIRecA model to a
four-dimensional tensor data used in [7] in order to accelerate the LA-SATF
model training.

## Acknowledgment

This work was supported by the Russian Science Foundation under Grant
22-21-00911.

## References

  * [1] W. Kang and J. J. McAuley, “Self-attentive sequential recommendation,” _CoRR_ , vol. abs/1808.09781, 2018. [Online]. Available: http://arxiv.org/abs/1808.09781
  * [2] R. Qiu, Z. Huang, H. Yin, and Z. Wang, “Contrastive learning for representation degeneration problem in sequential recommendation,” in _Proceedings of the Fifteenth ACM International Conference on Web Search and Data Mining_ , ser. WSDM ’22. New York, NY, USA: Association for Computing Machinery, 2022, p. 813–823. [Online]. Available: https://doi.org/10.1145/3488560.3498433
  * [3] A. Petrov and C. Macdonald, “Effective and efficient training for sequential recommendation using recency sampling,” in _Proceedings of the 16th ACM Conference on Recommender Systems_ , ser. RecSys ’22. New York, NY, USA: Association for Computing Machinery, 2022, p. 81–91. [Online]. Available: https://doi.org/10.1145/3523227.3546785
  * [4] P. Cremonesi, Y. Koren, and R. Turrin, “Performance of recommender algorithms on top-n recommendation tasks,” in _Proceedings of the Fourth ACM Conference on Recommender Systems_ , ser. RecSys ’10. New York, NY, USA: Association for Computing Machinery, 2010, p. 39–46. [Online]. Available: https://doi.org/10.1145/1864708.1864721
  * [5] O. Olaleke, I. V. Oseledets, and E. Frolov, “Dynamic modeling of user preferences for stable recommendations,” _CoRR_ , vol. abs/2104.05047, 2021. [Online]. Available: https://arxiv.org/abs/2104.05047
  * [6] K. Zhou, H. Wang, W. X. Zhao, Y. Zhu, S. Wang, F. Zhang, Z. Wang, and J.-R. Wen, “S3-rec: Self-supervised learning for sequential recommendation with mutual information maximization,” in _Proceedings of the 29th ACM International Conference on Information & Knowledge Management_, ser. CIKM ’20. New York, NY, USA: Association for Computing Machinery, 2020, p. 1893–1902. [Online]. Available: https://doi.org/10.1145/3340531.3411954
  * [7] E. Frolov and I. Oseledets, “Tensor-based sequential learning via hankel matrix representation for next item recommendations,” 2022.
  * [8] M. J. Pazzani and D. Billsus, _Content-Based Recommendation Systems_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2007, pp. 325–341. [Online]. Available: https://doi.org/10.1007/978-3-540-72079-9_10
  * [9] J. B. Schafer, D. Frankowski, J. Herlocker, and S. Sen, _Collaborative Filtering Recommender Systems_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2007, pp. 291–324. [Online]. Available: https://doi.org/10.1007/978-3-540-72079-9_9
  * [10] T. Donkers, B. Loepp, and J. Ziegler, “Sequential user-based recurrent neural network recommendations,” in _Proceedings of the Eleventh ACM Conference on Recommender Systems_ , ser. RecSys ’17. New York, NY, USA: Association for Computing Machinery, 2017, p. 152–160. [Online]. Available: https://doi.org/10.1145/3109859.3109877
  * [11] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” _SIAM Review_ , vol. 51, no. 3, pp. 455–500, 2009.
  * [12] S. Ahmadi-Asl, S. Abukhovich, M. G. Asante-Mensah, A. Cichocki, A. H. Phan, T. Tanaka, and I. Oseledets, “Randomized algorithms for computation of tucker decomposition and higher order svd (hosvd),” 2021.
  * [13] E. Frolov and I. V. Oseledets, “Tensor methods and recommender systems,” _CoRR_ , vol. abs/1603.06038, 2016. [Online]. Available: http://arxiv.org/abs/1603.06038
  * [14] O. Koch and C. Lubich, “Dynamical low-rank approximation,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 29, no. 2, pp. 434–454, 2007. [Online]. Available: https://doi.org/10.1137/050639703
  * [15] C. Lubich and I. Oseledets, “A projector-splitting integrator for dynamical low-rank approximation,” 2013. [Online]. Available: https://arxiv.org/abs/1301.1058
  * [16] O. Koch and C. Lubich, “Dynamical tensor approximation,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 31, no. 5, pp. 2360–2375, 2010. [Online]. Available: https://doi.org/10.1137/09076578X
  * [17] C. Lubich, “Time integration in the multiconfiguration time-dependent hartree method of molecular quantum dynamics,” _Applied Mathematics Research eXpress_ , vol. 2015, p. abv006, 06 2015.
  * [18] M. Brand, “Incremental singular value decomposition of uncertain data with missing values,” in _Computer Vision — ECCV 2002_ , A. Heyden, G. Sparr, M. Nielsen, and P. Johansen, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2002, pp. 707–720.
  * [19] V. Kalantzis, G. Kollias, S. Ubaru, A. N. Nikolakopoulos, L. Horesh, and K. Clarkson, “Projection techniques to update the truncated svd of evolving matrices with applications,” in _Proceedings of the 38th International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, M. Meila and T. Zhang, Eds., vol. 139. PMLR, 18–24 Jul 2021, pp. 5236–5246. [Online]. Available: https://proceedings.mlr.press/v139/kalantzis21a.html
  * [20] W. Hu, X. Li, X. Zhang, X. Shi, S. Maybank, and Z. Zhang, “Incremental tensor subspace learning and its applications to foreground segmentation and tracking,” _International Journal of Computer Vision_ , vol. 91, pp. 303–327, 02 2011.
  * [21] C. Lubich, B. Vandereycken, and H. Walach, “Time integration of rank-constrained tucker tensors,” 2017. [Online]. Available: https://arxiv.org/abs/1709.02594

TABLE I: Datasets statistics after initial preprocessing Dataset | n_users | n_items | n_actions | density%
---|---|---|---|---
ML-20M | 31670 | 26349 | 4000052 | 0.48
AMZ-B | 22363 | 12101 | 198502 | 0.07
AMZ-G | 19412 | 11924 | 167597 | 0.07
Steam | 281206 | 11961 | 3484694 | 0.10
TABLE II: Results of evaluation. All metrics are averaged by days and computed for top-n recommendations with n = 5 |  | ML-20M | AMZ-B | AMZ-G | Steam
---|---|---|---|---|---
Time | PSIRec | 0.376 | 0.722 | 0.591 | 0.657
| SVD | 7.129 | 4.420 | 1.513 | 5.413
| TDRec | 22.064 | 10.028 | 12.550 | 130.491
| TDRecReinit | 10.223 | 5.288 | 5.303 | 52.141
| TIRec | 27.681 | 39.237 | 31.858 | -
| TIRecA | 1.729 | 1.253 | 1.284 | 17.410
HR | PSIRec | 0.019 | 0.013 | 0.016 | 0.022
| SVD | 0.020 | 0.015 | 0.017 | 0.022
| TDRec | 0.026 | 0.020 | 0.012 | 0.015
| TDRecReinit | 0.027 | 0.020 | 0.013 | 0.020
| TIRec | 0.025 | 0.012 | 0.011 | -
| TIRecA | 0.025 | 0.013 | 0.010 | 0.009
MRR | PSIRec | 0.008 | 0.007 | 0.008 | 0.011
| SVD | 0.009 | 0.008 | 0.009 | 0.011
| TDRec | 0.013 | 0.010 | 0.006 | 0.007
| TDRecReinit | 0.013 | 0.011 | 0.006 | 0.010
| TIRec | 0.012 | 0.006 | 0.006 | -
| TIRecA | 0.013 | 0.007 | 0.005 | 0.004
WJI | PSIRec | 0.668 | 0.673 | 0.663 | 0.683
| SVD | 0.649 | 0.641 | 0.624 | 0.672
| TDRec | 0.092 | 0.030 | 0.029 | 0.111
| TDRecReinit | 0.241 | 0.231 | 0.167 | 0.488
| TIRec | 0.236 | 0.214 | 0.229 | -
| TIRecA | 0.429 | 0.224 | 0.210 | 0.539
TABLE III: Comparison of different initialization strategies for PSIRec model. All metrics are averaged and computed for top-n recommendations with n = 5. PSIRec($\sigma$;$Y$) depicts different model configurations, where $\sigma$ taken from $\mathcal{N}(0,\,\sigma^{2})$ initialization of factors; $Y=0$ \- zero initialization in the first step of Figure 1, $Y=\text{B}$ \- use Algorithm 2. |  | ML-20M | AMZ-B | AMZ-G | Steam
---|---|---|---|---|---
Time | PSIRec | 0.376 | 0.722 | 0.591 | 0.657
| PSIRec(0.001;0) | 0.677 | 1.397 | 1.164 | 0.690
| PSIRec(0.001;B) | 0.676 | 1.423 | 1.206 | 0.707
| PSIRec(0;0) | 0.326 | 0.731 | 0.619 | 0.417‘
| PSIRec(1e-05;0) | 0.666 | 1.389 | 1.185 | 0.690
| PSIRec(1e-05;B) | 0.636 | 1.404 | 1.194 | 0.700
| PSIRec(1e-08;0) | 0.631 | 1.393 | 1.178 | 0.674
| PSIRec(1e-08;B) | 0.690 | 1.409 | 1.181 | 0.686
HR | PSIRec | 0.019 | 0.013 | 0.016 | 0.022
| PSIRec(0.001;0) | 0.018 | 0.012 | 0.015 | 0.022
| PSIRec(0.001;B) | 0.018 | 0.012 | 0.015 | 0.022
| PSIRec(0;0) | 0.018 | 0.013 | 0.016 | 0.022
| PSIRec(1e-05;0) | 0.018 | 0.013 | 0.016 | 0.022
| PSIRec(1e-05;B) | 0.018 | 0.013 | 0.016 | 0.022
| PSIRec(1e-08;0) | 0.018 | 0.013 | 0.016 | 0.022
| PSIRec(1e-08;B) | 0.018 | 0.013 | 0.016 | 0.022
MRR | PSIRec | 0.008 | 0.007 | 0.008 | 0.011
| PSIRec(0.001;0) | 0.009 | 0.007 | 0.008 | 0.010
| PSIRec(0.001;B) | 0.008 | 0.006 | 0.008 | 0.011
| PSIRec(0;0) | 0.008 | 0.007 | 0.008 | 0.010
| PSIRec(1e-05;0) | 0.008 | 0.007 | 0.008 | 0.011
| PSIRec(1e-05;B) | 0.009 | 0.007 | 0.008 | 0.011
| PSIRec(1e-08;0) | 0.008 | 0.007 | 0.008 | 0.010
| PSIRec(1e-08;B) | 0.009 | 0.007 | 0.008 | 0.010
WJI | PSIRec | 0.668 | 0.673 | 0.663 | 0.683
| PSIRec(0.001;0) | 0.667 | 0.668 | 0.670 | 0.679
| PSIRec(0.001;B) | 0.673 | 0.674 | 0.665 | 0.676
| PSIRec(0;0) | 0.666 | 0.672 | 0.665 | 0.678
| PSIRec(1e-05;0) | 0.669 | 0.673 | 0.660 | 0.678
| PSIRec(1e-05;B) | 0.669 | 0.677 | 0.667 | 0.678
| PSIRec(1e-08;0) | 0.664 | 0.674 | 0.666 | 0.678
| PSIRec(1e-08;B) | 0.672 | 0.673 | 0.666 | 0.680
TABLE IV: Comparison of different initialization strategies for TIRec model. All metrics are averaged and computed for top-n recommendations with n = 5. TIRecA($\sigma$;$Y$) depicts different model configurations, where $\sigma$ taken from $\mathcal{N}(0,\,\sigma^{2})$ initialization of factors; $Y=0$ \- zero initialization in the first step of Figure 1, $Y=\text{B}$ \- use Algorithm 3. |  | ML-20M | AMZ-B | AMZ-G | Steam
---|---|---|---|---|---
Time | TIRec | 27.681 | 39.237 | 31.858 | -
| TIRecA | 1.729 | 1.253 | 1.284 | 17.410
| TIRecA(0.001;0) | 1.677 | 1.799 | 1.360 | 27.723
| TIRecA(0.001;B) | 2.079 | 2.171 | 2.246 | 27.847
| TIRecA(0;0) | 1.397 | 1.300 | 1.044 | 16.855
| TIRecA(1e-05;0) | 1.673 | 1.816 | 1.350 | 27.665
| TIRecA(1e-05;B) | 2.091 | 2.188 | 2.261 | 27.837
| TIRecA(1e-08;0) | 1.675 | 1.810 | 1.345 | 27.305
| TIRecA(1e-08;B) | 2.086 | 2.187 | 2.274 | 27.911
HR | TIRec | 0.025 | 0.012 | 0.011 | -
| TIRecA | 0.025 | 0.013 | 0.010 | 0.009
| TIRecA(0.001;0) | 0.024 | 0.014 | 0.010 | 0.009
| TIRecA(0.001;B) | 0.024 | 0.015 | 0.011 | 0.009
| TIRecA(0;0) | 0.024 | 0.014 | 0.011 | 0.009
| TIRecA(1e-05;0) | 0.024 | 0.014 | 0.012 | 0.009
| TIRecA(1e-05;B) | 0.024 | 0.014 | 0.011 | 0.009
| TIRecA(1e-08;0) | 0.024 | 0.014 | 0.011 | 0.009
| TIRecA(1e-08;B) | 0.024 | 0.014 | 0.011 | 0.009
MRR | TIRec | 0.012 | 0.006 | 0.006 | -
| TIRecA | 0.013 | 0.007 | 0.005 | 0.004
| TIRecA(0.001;0) | 0.012 | 0.007 | 0.005 | 0.004
| TIRecA(0.001;B) | 0.012 | 0.008 | 0.005 | 0.004
| TIRecA(0;0) | 0.012 | 0.007 | 0.005 | 0.004
| TIRecA(1e-05;0) | 0.012 | 0.007 | 0.006 | 0.004
| TIRecA(1e-05;B) | 0.012 | 0.008 | 0.006 | 0.004
| TIRecA(1e-08;0) | 0.012 | 0.007 | 0.005 | 0.004
| TIRecA(1e-08;B) | 0.012 | 0.007 | 0.006 | 0.004
WJI | TIRec | 0.236 | 0.214 | 0.229 | -
| TIRecA | 0.429 | 0.224 | 0.210 | 0.539
| TIRecA(0.001;0) | 0.424 | 0.223 | 0.210 | 0.536
| TIRecA(0.001;B) | 0.425 | 0.219 | 0.214 | 0.542
| TIRecA(0;0) | 0.427 | 0.226 | 0.206 | 0.539
| TIRecA(1e-05;0) | 0.427 | 0.222 | 0.208 | 0.547
| TIRecA(1e-05;B) | 0.425 | 0.220 | 0.210 | 0.534
| TIRecA(1e-08;0) | 0.427 | 0.225 | 0.206 | 0.545
| TIRecA(1e-08;B) | 0.429 | 0.219 | 0.209 | 0.542
Algorithm 5 PSIRec: Add new users or new items

1:

2:\- $\Delta A$ \- data representing new users/items as columns

3:\- $U_{0}$, $S_{0}$, $V_{0}^{\top}$ \- state of the model

4:$U_{1}$, $S_{1}$, $V_{1}^{\top}$ \- updated state of the model

5:J, K = QR($\Delta A-U_{0}U_{0}^{\top}\Delta A$)

6:$U^{\prime},S^{\prime},V^{\prime\top}$ =
SVD$\left(\begin{bmatrix}S_{0}&U_{0}^{\top}\Delta A\\\
0&K\end{bmatrix}\right)$

7:$U_{1}\leftarrow[U_{0}J]U^{\prime}$

8:$S_{1}\leftarrow S^{\prime}$

9:$V_{1}\leftarrow\begin{bmatrix}V_{0}&0\\\ 0&I\end{bmatrix}V^{\prime}$

Algorithm 6 TIRec: Add new users or new items

1:

2:\- $\Delta\mathcal{A}$ \- data tensor representing new users or new items
interactions

3:\- $U_{1}^{0}$, $U_{2}^{0}$, $U_{3}^{0}$, $\mathcal{C}^{0}$ \- state of the
model

4:\- $m$ \- target mode: 1 for new users, 2 for new items

5:$U_{1}^{1},U_{2}^{1},U_{3}^{1},\mathcal{C}^{1}$ \- updated state of the
model

6:$\mathcal{C}_{[m]}^{0}\left(U_{3}^{0}\otimes
U_{2}^{0}\right)^{\top}=\tilde{U}\tilde{S}\tilde{V}^{\top}$

7:$\hat{U}=U_{1}^{0}\tilde{U}$

8:$\bar{V},\bar{S},\bar{U}^{\top}=$
update_vectors_2d($\Delta\mathcal{A}_{[m]}^{\top}$, $\tilde{V}$, $\tilde{S}$,
$\hat{U}^{\top}$)

9:$\mathcal{C}_{[m]}^{1}=\bar{S}\bar{V}^{\top}\left(U_{3}^{0}\otimes
U_{2}^{0}\right)$


11institutetext:  Institute of Computing Science, Faculty of Computing and
Telecommunications, Poznan University of Technology, Piotrowo 2, 60-965
Poznan, Poland 22institutetext: MNM Bioscience Inc., 1 Broadway, Cambridge, MA
02142 33institutetext: Institute of Bioorganic Chemistry of the Polish Academy
of Sciences, Zygmunta Noskowskiego 12/14, 61-704 Poznan, Poland
44institutetext: Physics/Biology Department, Faculty of Physics, Adam
Mickiewicz University, Uniwersytetu Poznanskiego 2, 61-614 Poznan, Poland
44email<EMAIL_ADDRESS>
# Random Similarity Forests

Maciej Piernik 1122 0000-0001-7281-284X Dariusz Brzezinski 112233
0000-0001-9723-525X Pawel Zawadzki 2244 0000-0002-9032-2315

###### Abstract

The wealth of data being gathered about humans and their surroundings drives
new machine learning applications in various fields. Consequently, more and
more often, classifiers are trained using not only numerical data but also
complex data objects. For example, multi-omics analyses attempt to combine
numerical descriptions with distributions, time series data, discrete
sequences, and graphs. Such integration of data from different domains
requires either omitting some of the data, creating separate models for
different formats, or simplifying some of the data to adhere to a shared scale
and format, all of which can hinder predictive performance. In this paper, we
propose a classification method capable of handling datasets with features of
arbitrary data types while retaining each feature’s characteristic. The
proposed algorithm, called Random Similarity Forest, uses multiple domain-
specific distance measures to combine the predictive performance of Random
Forests with the flexibility of Similarity Forests. We show that Random
Similarity Forests are on par with Random Forests on numerical data and
outperform them on datasets from complex or mixed data domains. Our results
highlight the applicability of Random Similarity Forests to noisy, multi-
source datasets that are becoming ubiquitous in high-impact life science
projects.

##  1 Introduction

Over the last decades, machine learning has matured to the point in which
researchers are spoiled for choice with classifiers for scalar data. There are
hundreds of classifiers to choose from and their performance has been tested
on thousands of datasets spanning across various domains [6]. Increasingly,
however, researchers are more and more often faced with datasets of complex
composition, where numerical features are not explicitly available. In such
cases, there are two main options to choose from. Either to transform complex
objects into many simple features, which can then be used by any traditional
feature-based classifier (e.g., Random Forests [3]) or a dedicated one (e.g.,
CNNs for images [9]), or to rely on a domain-specific distance measure and use
distance-based classifiers (e.g., Similarity Forests [15]). However, there is
a third option that relies on the fact that many complex data types can be
broken down into both traditional features (e.g., numerical or categorical)
and simpler but still complex substructures for which good distance measures
exist. For example, sequences can be decomposed into subsequences, sets of
elements, distributions of elements, or distributions of lengths of
subsequences. These complex substructures can then be combined together with
simple, numerical or categorical, features and form a mixed type dataset to
better capture the detailed characteristics of each example. Unfortunately,
the currently available classification methods inherently enable only one of
the first two types of analyses, on scalar data or complex objects. This
limitation may prohibit us from using the full potential hidden in the data,
as it forces us to either reduce all complex structures into simple features
(most probably with information loss) or process all features with a single,
complex distance measure (which is prone to the curse of dimensionality).

Consider a scenario of building a tumor genome classifier predicting patients
responding to a given treatment based on whole-genome sequencing (WGS) [1].
Although, theoretically, a distance-based approach can be used in this
scenario, it is unlikely to produce good results because of the curse of
dimensionality and computational infeasibility (as the reference human genome
is ~3.2 billion base pairs long). Moreover, this approach also neglects
decades of research in the field of oncology and genomics, which provide a
wide array of potentially useful biomarkers. A more natural approach in this
scenario is to use a feature-based classifier that can make use of existing
biomarkers as well as explore new ones. Currently however, this approach
necessitates a drastic decomposition of the data into a fixed set of scalar
features. At the same time, the information in the human genome can be
decomposed into many intermediate structures carrying significantly more
information than their scalar counterparts. Examples of such structures
include variant distributions, sequential patterns, gene sets, interaction
graphs. All of these structures have multiple distance measures available and
could be used in conjunction with existing numerical biomarkers as a single
dataset to produce more informed predictions. However, to the best of our
knowledge, there does not exist any classifier capable of handling data with
mixed, arbitrarily defined types of features.

In this paper, we propose a new classifier, called Random Similarity Forests
(RSF), capable of handling datasets with mixed, arbitrarily defined feature
types. Our approach requires a distance measure to be provided for each
feature and works as a blend of Random Forests and Similarity Forests. By
using a mixture of feature types, Random Similarity Forests inherit the
benefits of both feature-based and similarity-based methods. In the paper, we
intuitively describe and formally define the proposed classifier. We also
perform sensitivity tests with respect to the number of features analyzed in
each tree node and the number of trees in the forest. The method is then
experimentally evaluated against Random Forests and Similarity Forests on
datasets with scalar features, complex objects with a distance measure, and
mixed, arbitrarily defined features with distance measures. Our analysis
highlights the characteristic properties of each approach and discusses their
suitability for different dataset types.

##  2 Related work

Random forests [3] are one of the most popular classifiers and have been shown
to offer the best quality of predictions across a wide range of datasets when
compared to other classifiers [6]. However, like many other approaches, their
use is limited to scalar data and, therefore, requires feature extraction when
dealing with complex objects. Due to this fact, Similarity Forests [15] have
been recently proposed as an alternative to Random Forests, especially when
regular features are absent but similarities or distances between examples are
attainable.

When dealing with distance-based classification, kernel methods are a popular
group to consider, especially since there have been studies exploring the
relationship between Random Forests and kernel methods [16]. In general,
distance-based classification is often applied when dealing with complex data,
e.g., for microbiome, image, or time series classification. A different way of
dealing with complex data in classification using distances is through a
combination of clustering and distance-based feature extraction. This approach
relies on unsupervised clustering of the data and later encoding the clusters
as features, either through aggregate inter- and intra-cluster distances [17]
or simply by encoding the distances to each cluster center as separate
features [13].

Recently, a different combination of feature-based and distance-based methods
has been proposed by Liang et al. [10], where the authors first use Random
Forests to select the most important features according to their impurity-
based importance and later use these features in a k-Nearest Neighbors
classifier with dynamically selected distance measures.

Although the described existing approaches create forests or calculate
distances between objects, they either require scalar features or are based on
a single distance calculated for entire examples. Similarly, integrative
approaches from biomedical studies that combine data from different domains
either require a common feature representation, create separate data type
specific models, or perform data re-scaling and simplification [14]. To the
best of our knowledge, Random Similarity Forests are the only single-
classifier model capable of dealing with datasets with a mixture of
arbitrarily defined features without any data transformations.

##  3 Random Similarity Forest

###  3.1 Preliminaries

Given a dataset of training examples
$\mathcal{X}=\\{\boldsymbol{x_{1}},\boldsymbol{x_{2}},\ldots,\boldsymbol{x_{n}}\\}$
and their corresponding class labels
$\boldsymbol{y}=\\{y_{1},y_{2},\ldots,y_{n}\\}$, the task of a classifier is
to predict the class label $\hat{y}$ of each unlabeled example
$\boldsymbol{\hat{x}}$ (i.e., an example for which the class is unknown).
Every example $\boldsymbol{x_{i}}$ is a list
$(x_{i1},\ldots,x_{ij},\ldots,x_{ip})$ of the same length $p$, where each
position $j$ holds a value a given example has for a particular feature
$F_{1},F_{2},\ldots,F_{p}\in\mathcal{F}$. We also denote all values of
$\mathcal{X}$ across a given feature $F_{j}$ by $\boldsymbol{x_{.j}}$. Each
class label $y_{i}$ falls into one of several categorical class values
$\mathcal{C}=\\{c_{1},c_{2},\ldots,c_{m}\\}$. This is the standard definition
of a classification problem. In our case, each feature can be from a
different, arbitrary data domain. Additionally, every feature
$F_{1},F_{2},\ldots,F_{p}\in\mathcal{F}$ has an associated distance measure
$\delta_{1},\delta_{2},\ldots,\delta_{p}\in\Delta$. Therefore, the distance
between any two examples $\boldsymbol{x_{i}}$, $\boldsymbol{x_{j}}$ for a
given feature $F_{l}$ is equal to $\delta_{l}(x_{il},x_{jl})$, or
$\delta_{l}(\boldsymbol{x_{i}},\boldsymbol{x_{j}})$ for short.

###  3.2 The Algorithm

Just like Random Forests and Similarity Forests, Random Similarity Forests
rely on bagging of a user-defined number (max_trees) of single-tree
classifiers. A single Random Similarity Tree is constructed in a top-down
fashion by recursively splitting the examples in each node into two exclusive
subsets. The recursion stops when a given node contains only examples from a
single class or when one of the early-stopping conditions is met (discussed in
Section 3.4). Each split is calculated using a similarity-based 1-dimensional
projection of the examples in a given node. First, a subset of max_features
candidate features is selected, and afterward, max_pairs pairs of examples are
picked at random for each feature. The examples in a pair are selected so that
they come from different classes and have a different value on a given
feature.

For a given pair of examples $\boldsymbol{x_{p}}$, $\boldsymbol{x_{q}}$, the
distance between all other examples
$\boldsymbol{x_{i}}\in\mathcal{X}_{\mathit{node}}\setminus\\{\boldsymbol{x_{p}}$,
$\boldsymbol{x_{q}}\\}$ and the two selected examples is evaluated for each
selected feature $F_{l}$ in order to create a projection
$\mathcal{P}(\boldsymbol{x_{i}})$ of each example $\boldsymbol{x_{i}}$ into a
direction defined by $\boldsymbol{x_{p}}$ and $\boldsymbol{x_{q}}$. As shown
by Sathe and Aggarwal [15], for the purposes of constructing a split, this
projection can be approximated with
$\mathcal{P}(\boldsymbol{x_{i}})\propto\delta_{l}(\boldsymbol{x_{q}},\boldsymbol{x_{i}})-\delta_{l}(\boldsymbol{x_{p}},\boldsymbol{x_{i}})$,
as this approximation preserves the order of the original projection and the
order is the only information necessary to construct a split. Therefore, for
each pair of examples $\boldsymbol{x_{p}}$, $\boldsymbol{x_{q}}$, all
remaining examples $\boldsymbol{x_{i}}$ in a given node are ordered by
$\delta_{l}(\boldsymbol{x_{q}},\boldsymbol{x_{i}})-\delta_{l}(\boldsymbol{x_{p}},\boldsymbol{x_{i}})$
along a given feature $F_{l}$.

This projection approximation can be viewed as a dynamic feature. After the
projection has been computed, it can be used to perform a split just like any
typical numerical feature would be used in a regular decision tree. The split
point is selected as the point which minimizes the weighted average of the
Gini index of the children nodes. For a given node $N_{i}$, the Gini index is
defined as $G(N_{i})=1-\sum_{c\in\mathcal{C}}{p_{c}^{2}}$, where $p_{c}$ is
the fraction of examples of class $c$ in a given node. Therefore, given two
nodes $N_{i}$, $N_{j}$ consisting of $n_{i}$, $n_{j}$ examples, respectively,
the weighted Gini index is calculated as:

$G(N_{i},N_{j})=\frac{n_{i}G(N_{i})+n_{j}G(N_{j})}{n_{i}+n_{j}}\text{.}$ (1)

Here, we opted for the Gini-index-based splitting strategy used also in
Similarity Forests, however, it is important to note that after the projection
is done we are in fact in possession of a dynamically calculated numerical
feature. Therefore, other measures, such as entropy, can also be used.

Given the order defined with the projection and an impurity measure, all
possible splits along a given feature are evaluated and the split that
minimizes the impurity is selected. After performing this procedure for all
selected features in a given node, the feature producing the best split is
selected. Subsequently, two new tree nodes are created and the examples
falling below the split threshold are placed in one node while the remaining
examples are placed in the other node, and the whole procedure is repeated
until all leaf nodes are pure or an early-stopping criterion is met.

The above-described process is illustrated in Fig. 1. The pseudocode of
training a single Random Similarity Tree is presented in Algorithm 1.

Figure 1: An example illustrating the training process of Random Similarity
Forests. Gray areas symbolize the bootstrap samples and subset of features
going into each tree. Feature $F_{3}$ (A) is selected to perform a split in
one of the trees (B). The distance-based projection (C) for that feature is
calculated and the split point yielding the purest split is selected. The
reference points (distributions) $\boldsymbol{x_{p}},\boldsymbol{x_{q}}$ used
to create the projection are underlined in red and blue in the table.
Algorithm 1 Random Similarity Tree: $\mathit{RST}(\cdot)$

0: $\mathcal{X}$: training examples; $\boldsymbol{y}$: examples’ class labels;
$\boldsymbol{\delta}$: distance measures for each feature; max_features:
number of features picked randomly at each node; max_pairs: number of example
pairs picked at random for each feature at each node

0: a trained Random Similarity Tree model

1: if $\mathit{earlyStopping()}$ or $len(unique(\boldsymbol{y}))=1$ then

2: $\mathit{setLeaf()}$

3: return $\mathit{self}$

4: else

5: for $1..\mathit{max\\_features}$ do

6:
$F_{j}\leftarrow\mathit{random}(\\{F_{j}\in\mathcal{F}:\mathit{var}(\boldsymbol{x_{.j}})>0\\})$

7:
$c_{1}\leftarrow\operatorname*{\arg\\!\min}_{c\in\mathit{unique}(\boldsymbol{y})}\\{\mathit{var}(\\{x_{ij}:y_{i}=c\\})\\}$

8:

9: for $1..\mathit{max\\_pairs}$ do

10:
$\boldsymbol{x_{p}}\leftarrow\mathit{random}(\\{\boldsymbol{x_{i}}:y_{i}=c_{1}\\})$

11:
$\boldsymbol{x_{q}}\leftarrow\mathit{random}(\\{\boldsymbol{x_{i}}:y_{i}\neq
c_{1}\wedge x_{ij}\neq x_{pj}\\})$

12: for $\boldsymbol{x_{i}}\in\mathcal{X}$ do

13:
$P(\boldsymbol{x_{i}})\leftarrow\delta_{j}(\boldsymbol{x_{q}},\boldsymbol{x_{i}})-\delta_{j}(\boldsymbol{x_{p}},\boldsymbol{x_{i}})$

14: for $\mathit{thr}\in\mathit{unique}(P)$ do

15:
$\boldsymbol{y_{\mathit{left}}}\leftarrow\\{y_{i}\in\boldsymbol{y}:P(\boldsymbol{x_{i}})\leq\mathit{thr}\\}$

16:
$\boldsymbol{y_{\mathit{right}}}\leftarrow\\{y_{i}\in\boldsymbol{y}:P(\boldsymbol{x_{i}})>\mathit{thr}\\}$

17: $g\leftarrow
G(\boldsymbol{y_{\mathit{left}}},\boldsymbol{y_{\mathit{right}}})$

18: if $(g<best_{g})$ or ($g=best_{g}$ and
$|len(y_{\mathit{left}})-len(y_{\mathit{right}})|<best_{bal}$) then

19:
$update\\_best(g,\mathit{thr},F,\boldsymbol{x_{p}},\boldsymbol{x_{q}},P,bal)$

20:

21:
$\mathcal{X}_{\mathit{left}},\boldsymbol{y_{\mathit{left}}}\leftarrow\\{\boldsymbol{x_{i}}\in\mathcal{X},y_{i}\in\boldsymbol{y}:P(\boldsymbol{x_{i}})\leq\mathit{thr}\\}$

22:
$\mathcal{X}_{\mathit{right}},\boldsymbol{y_{\mathit{right}}}\leftarrow\\{\boldsymbol{x_{i}}\in\mathcal{X},y_{i}\in\boldsymbol{y}:P(\boldsymbol{x_{i}})>\mathit{thr}\\}$

23:
$N_{\mathit{left}}\leftarrow\mathit{RST}(\mathcal{X}_{\mathit{left}},\boldsymbol{y_{\mathit{left}}},\boldsymbol{\delta},\mathit{max\\_features},\mathit{max\\_pairs})$

24:
$N_{\mathit{right}}\leftarrow\mathit{RST}(\mathcal{X}_{\mathit{right}},\boldsymbol{y_{\mathit{right}}},\boldsymbol{\delta},\mathit{max\\_features},\mathit{max\\_pairs})$

25: return $\mathit{self}$

After all max_trees Random Similarity Trees are computed, the trained Random
Similarity Forest is ready to make predictions. For each node in each tree,
given an unlabeled example $\boldsymbol{\hat{x}}$, the example is projected
along the feature and the two examples stored in that node and assigned to one
of the child nodes based on the stored split point. Following this procedure,
once the example reaches a leaf node, it returns a weighted prediction of the
majority class in that node. After predictions from all trees in the forest
are made, the example is assigned to the class with the highest weighted
average from all single-tree predictions.

###  3.3 Computational Complexity

Let us now discuss the computational complexity of Random Similarity Forests
(RSF). Unsurprisingly, it is strongly related to the complexity of Similarity
Forests (SF), which is $\mathcal{O}(n\cdot\log{n})$. However, RSF additionally
adds the cost of checking multiple features at each node, so the complexity is
raised by that factor to $\mathcal{O}(p\cdot n\cdot\log{n})$, which is
equivalent to the complexity of Random Forests (RF). However, it is worth
noting that in practice the value of the max_features parameter for RSF may
need to be higher than for RF as it is inherently more random as it randomly
picks both features and pairs of examples. Moreover, even though the number of
pairs is a parameter, we recommend picking only a single pair for each
feature, as suggested by Sathe and Aggrawal [15].

Another cost hidden in the computational complexity is the cost of calculating
the distances between examples. It is an essential factor because, depending
on the feature type, it can either be negligible or a very costly operation.
For simple numerical features, the distance computation can be usually omitted
as for any metric it would produce a projection equivalent to the original
feature in terms of the example order. On the other hand, more complex data
types may require substantial computation, e.g., edit distance has quadratic
complexity for sequential data and is NP-complete for graphs. Therefore, the
use of such metrics may be prohibitive for longer sequences or larger graphs.

###  3.4 Discussion

Since Random Similarity Forests draw from Random Forests and Similarity
Forests, they combine their advantages and possess unique properties:

  1. 1.

By using feature-oriented similarity metrics to split data points, Random
Similarity Forests can be used to classify datasets characterized by numerical
as well as complex data features. Random Forests cannot classify complex data
objects, whereas Similarity Forests treat entire examples as data objects,
hiding the characteristics of individual features.

  2. 2.

By sampling feature subsets for each node split and by analyzing each feature
separately, Random Similarity Forests are robust to noisy datasets with a lot
of features. Similarity Forests treat entire examples as data objects, which
can hinder their performance on datasets with many irrelevant features.

  3. 3.

Analyzing similarities based on each feature separately allows Random
Similarity Forests to use several similarity measures to analyze the same
objects and let the forest decide which measures are the most useful for
classification.

In addition to combining the advantages from both Random Forests and
Similarity Forests, Random Similarity Forests also inherit some of their
limitations. The first limitation comes from Similarity Forests and it is the
fact that it is designed for binary classification problems. However, one can
easily adapt it to multi-class classification or regression problems
analogously to the method proposed for Similarity Forests [4], or by using the
one-vs-all approach.

Another characteristic worth discussing is the effect of feature variance on
node splits. If one would select the pair of examples $(\boldsymbol{x_{p}}$,
$\boldsymbol{x_{q}})$ with the same value of feature $F_{j}$, then they would
produce a completely random ordering of examples for the split (since if
$x_{pj}=x_{qj}$ then
$\forall{\boldsymbol{x_{i}}}:\delta_{j}(\boldsymbol{x_{q}},\boldsymbol{x_{i}})-\delta_{j}(\boldsymbol{x_{p}},\boldsymbol{x_{i}})=0$).
Similar problems appear in Random Forests when the feature values are the same
for examples of different classes, whereas Similarity Forests do not suffer
from this particular limitation so often as they always consider each example
as a whole. For datasets with many features with low variance, this reduces
the effective number of features tested at each node and promotes high-
variance features. That is why Random Similarity Forests will perform better
on datasets with features with high variance. This limitation can be overcome
by tuning the max_features hyperparameter. As shown by Geurts et al. [7], the
best value for this parameter highly depends on the number of correlated and
irrelevant features: lower values work better with highly correlated features
while higher values with irrelevant features. Indeed, as will be shown in
Section 4.3, the proposed default value of max_features for Random Similarity
Forests is $0.5p$, which is higher than $\sqrt{p}$ commonly used for Random
Forests [7].

Other issues worth mentioning are the early-stopping conditions (pre-pruning)
and pruning (post-pruning). There are many possible options to stop the tree
construction process before it produces perfectly pure leaf nodes. These could
include, e.g., limiting the maximum depth a tree can reach, specifying the
minimum number of elements required to perform a split or the minimum number
of elements to form a leaf node. Similarly, there are multiple ways in which a
tree can be pruned, with arguably the most popular method being minimal cost
complexity pruning [2]. However, pre-and post-pruning methods are not specific
to the proposed algorithm or any other forest algorithm and are out of the
main scope of this paper.

##  4 Experimental evaluation

In this section, we experimentally evaluate the Random Similarity Forests
algorithm and illustrate its usefulness in real-world scenarios. Since the
proposed classifier is related to both Random Forests and Similarity Forests,
the nature of the evaluation is highly comparative. Our goal is to showcase
the strengths and limitations of our proposal and indicate when its use would
be the most beneficial and when it might be preferable to rely on one of the
alternatives. First, we compare the classification performance of the three
classifiers on ten publicly available datasets with simple numerical features.
Next, we compare the three approaches on complex data, i.e., datasets
consisting of arbitrarily defined objects with distance measures. For this
purpose, we use sequences of sets, time series, and graph datasets. Finally,
we illustrate the usefulness of our proposal by training the classifiers on
datasets consisting of mixtures of complex object-like features and simple
numerical features.

###  4.1 Experimental Setup

In all of the experiments, we compare Random Similarity Forests (RSF) against
Random Forests (RF) and Similarity Forests (SF) using the area under the ROC
curve (AUC). We chose AUC since we focus on binary classification, it is skew
invariant, assesses the classifiers’ ranking abilities, and it has been shown
to be statistically consistent and more discriminating than accuracy [8]. To
evaluate each approach on each dataset, we rely on 10 repetitions of
stratified 2-fold cross-validation (10$\times$2CV). Afterwards, the results
undergo a series of statistical tests [5]. We rely on the Friedman statistic
with a post-hoc Nemenyi test to distinguish between the compared approaches
across all datasets. Afterwards, we check for differences on each individual
dataset. First, we verify the normality assumption using the Shapiro-Wilk
test. If the assumption is met, we use the ANOVA statistic with a post-hoc
Tukey HSD test. If the assumption is not met, we once again resort to the non-
parametric Friedman statistic with a post-hoc Nemenyi test [5]. All tests were
performed at significance level $\alpha=0.05$.

The code used to perform the experiments was written in Python with parts of
the implementation of Random Similarity Forests written in Cython to achieve
higher efficiency. For Random Forests, we use the scikit-learn implementation
[12], whereas for Similarity Forests, we rely on the implementation described
in [4]. Unless stated otherwise, we rely on the default hyperparameters for
each classifier, as specified in their implementations.

###  4.2 Datasets

Since our experimental analysis focuses on data with scalar, complex, and
mixed types of features, we use three groups of datasets, each corresponding
to a given type of features. Additionally, we have a separate set of datasets
dedicated to performing sensitivity analyses. The choice of datasets was made
a priori and independently of the results obtained with our methods.

#### Datasets for sensitivity analysis

For sensitivity analysis, we used binarized versions of 12 numerical datasets:
ten from the UCI Machine Learning111https://archive.ics.uci.edu/ml/ repository
and two introduced by Breiman; for a detailed description of these datasets
see [7]. We chose these datasets as they were used to discuss the impact of
different hyperparameter values for Random Forests and Extra Trees [7].
Importantly, these datasets are independent of those used for assessing the
predictive performance of RF, SF, and RSF, to avoid biasing the comparative
analysis.

#### Scalar data

For the analysis of scalar data, we used ten publicly available datasets with
numeric features. The datasets exhibit various conditions in terms of the
number of features, learning sample size, and feature variance, and can be
accessed through the LIBSVM dataset
repository.222https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/ We used
the scaled versions of the datasets, as required by Similarity Forests. Each
dataset was additionally preprocessed by removing features containing missing
values. Table 1 presents the main characteristics of each dataset after
preprocessing.

#### Complex data

For the analysis of complex data (i.e., data where each example is an object
without any explicitly defined numerical features), we used three groups of
datasets consisting of sequences of sets, time series, and graphs,
respectively.

Sequences of sets are a good candidate for complex data because they combine
set information, which can be easily encoded with numerical features, and
sequential information, which is very difficult to encode as scalars but
manifests itself through distance measures on whole objects. Therefore, we
used a synthetic data
generator333https://gingerbread.shinyapps.io/SequencesOfSetsGenerator/ to
create three datasets of sequences of sets: items, lengths, and order. The
sequences in the items dataset have the same mean length and mean set size,
and are differentiated by the distributions of the selected items between the
two classes. Sequences in the lengths dataset are generated with identical set
size and item distributions but are differentiated by the lengths between the
classes. Finally, all the sequences in the order dataset are generated from
the same length, set size, and item distributions, and are only differentiated
by the order of the elements in one of the classes.

Time series datasets were taken from the UCR Time Series Classification
Archive444https://www.cs.ucr.edu/~eamonn/time_series_data_2018/, with
distances between objects calculated with dynamic time warping. The graph
datasets were taken from the TU-Dortmund Graph Kernel
Benchmarks555https://ls11-www.cs.tu-
dortmund.de/staff/morris/graphkerneldatasets with distances between objects
calculated using the Ipsen-Mikhailov distance. In both cases, we limited our
selection to binary classification problems, with varying numbers of examples
and features (length, sets, nodes, edges).

#### Mixed data

For the analysis of mixed data (i.e., data where each example is described by
arbitrarily defined features with a distance measure for each feature), we
used four groups of datasets consisting of time series, graphs, multi-omics
data, and genome sequencing data.

Time series datasets consisted of three distance measures (euclidean, cosine,
dynamic time warping) computed for series from the UCR
Repository.666https://www.cs.ucr.edu/~eamonn/time_series_data_2018/ Similarly,
the graph datasets consisted of four distance measures (portrait divergence,
Jaccard distance, degree divergence, Ipsen-Mikhailov distance) computed for
graphs from the TU-Dortmund repository777https://ls11-www.cs.tu-
dortmund.de/staff/morris/graphkerneldatasets using the netrd package [11]. The
multiomics datasets consisted of Alzheimer’s Disease and BRCA-mutated breast
cancer patients, described by three sets of numeric expression
data.888https://github.com/txWang/MOGONET Finally, the last group of datasets
represented tumor samples of ovarian and breast cancer patients obtained using
whole genome sequencing.999https://github.com/MNMdiagnostics/dbfe In this
group of datasets, each sample is characterized by 50 gene amplifications
represented as simple numerical features and 6 distributions of lengths of
various DNA structural variants.

###  4.3 Sensitivity Analysis

Using the 12 datasets proposed in [7], we have analyzed the effect of
max_features on Random Similarity Forests. The parameter denotes the number of
potential splits screened at each node during the growth of a Random
Similarity Tree. It may be chosen in the interval $[1,\ldots,p]$ and the
smaller its value the stronger the randomization of the trees. Figure 2
compares the AUC of Random Forests (RF) and Random Similarity Forests (RSF)
for increasing values of max_features.

Figure 2: Evolution of mean AUC of Random Forests (RF) and Random Similarity
Forests (RSF) with varying max_features on 12 datasets. Dashed line shows the
proposed RSF default value of max_features=$0.5p$.

One can notice that for all the datasets the performance of RSF first
increases and then declines as max_features approaches the total number of
features in a dataset $p$. As mentioned in Section 3.4, this phenomenon is
related to the number of correlated and irrelevant features in a dataset:
lower values work better with highly correlated features while higher values
of max_features work better with irrelevant features. This is due to the fact
that a higher value of max_features leads to a better chance of filtering out
the irrelevant variables. These findings are in accordance with [7], were
similar trends were noticed for the Extra Trees classifier. It is worth
noting, however, that compared to RF, the increase of predictive performance
with max_features is slower. This stems from the fact that RSF samples pairs
of examples and creates different projections, thus introducing additional
randomness. Therefore, we propose to use a default value of
max_features=$0.5p$, as opposed to $\sqrt{p}$, which is the commonly used
default for Random Forests and Extra Trees [7].

Figure 3: Mean AUC of Random Forests (RF), Random Similarity Forests (RSF) and
Similarity Forests (SF) with varying max_trees on the Waveform dataset.

We have also analyzed the max_trees parameter, responsible for the number of
trees in the ensemble (Fig. 3). The results are in accordance with studies on
different tree ensembles [7]: the higher the value of max_trees, the better
the ensemble’s accuracy. It can be noticed that all three algorithms (RF, SF,
RSF) have the same speed of convergence. In the remaining experiments, for all
three classifiers, we use max_trees=100, which is large enough to ensure
convergence of the ensemble effect on all our datasets.

###  4.4 Evaluation on Scalar Data

The aim of the first set of experiments was to assess the predictive
performance of Random Similarity Forests (RSF) on scalar data. Since scalar
features are the natural data domain of Random Forests (RF), the aim is also
to verify whether the proposed classifier will be able to match the
performance of Random Forests. Since Similarity Forests (SF) are also capable
of dealing with scalar data and were reported to produce results on par with
or even better than Random Forests [15], we include this approach in this
comparison as well. The results of these experiments are presented in Table 1.

[b]

Table 1: Dataset characteristics and mean AUC performance of Random Forests
(RF), Similarity Forests (SF), and Random Similarity Forests (RSF) from 10
$\times$ 2-fold cross-validation experiments. Features for sequences of sets
denoted as [mean set size]/[mean length]. Features for graphs denoted as [mean
number of nodes]/[mean number of edges].

Dataset | Type | Features | #Ex. | #Features | Distance | AUC
---|---|---|---|---|---|---
RF |  RSF |  SF
australian | scalar | numeric | 690 | 8 | euclidean |  0.92∘ |  0.91 |  0.92
breast | 683 | 10 |  0.99 |  0.99 |  1.00∘
diabetes | 768 | 6 |  0.80 |  0.80 |  0.80
german | 1000 | 18 |  0.73 |  0.73 |  0.73
heart | 270 | 9 |  0.88 |  0.87 |  0.89∘
ionosphere | 351 | 2 |  0.84 |  0.84 |  0.84
liver | 145 | 4 |  0.77 |  0.75 |  0.76
sonar | 208 | 59 |  0.89 |  0.88 |  0.91∘
splice | 1000 | 60 |  0.99∘ |  0.99∘ |  0.90
svmguide3 | 1243 | 21 |  0.84∘ |  0.85∘ |  0.77
items | complex |  sequences of sets | 400 | 20/20 |  edit distance |  1.00∙b |  0.91∘ |  0.85
length | 400 | 20/20 |  0.86b |  0.91∘ |  0.92∘
order | 400 | 20/20 |  0.56b |  0.79∙ |  0.73∘
computers |  time series | 500 | 720 |  dynamic time warping |  0.68 |  0.77∘ |  0.78∘
housetwenty | 159 | 2000 |  0.92 |  0.99∘ |  0.99∘
toeseg | 268 | 277 |  0.85 |  0.96∘ |  0.97∘
cox2 |  graph | 467 | 41/43 |  ipsen-mikhailov |  0.65b |  0.70∘ |  0.71∘
mutag | 188 | 18/20 |  0.87b |  0.92∘ |  0.92∘
proteins | 1113 | 39/73 |  0.77b |  0.95∙ |  0.91∘
ptcfm | 349 | 14/14 |  0.58b |  0.63∘ |  0.61∘
beetlefly | mixed | time series | 40 | 512 |  euclidean, cosine, dtw |  0.84 |  0.85 |  0.81m
birdchicken | 40 | 512 |  0.89 |  0.90 |  0.87m
bzr | graphs | 405 | 36/38 |  portrait, degree, jaccard, ipsen |  0.73b |  0.90∙ |  0.87m∘
dhfr | 756 | 42/45 |  0.78b |  0.85∙ |  0.84m∘
rosmap |  multiomics (numeric) | 351 | 600 |  euclidean, cosine |  0.78 |  0.79∘ |  0.76m
brca | 875 | 2503 |  0.89 |  0.91∘ |  0.91m∘
wgs_ovarian |  distributions, numeric | 219 | 56 |  euclidean |  0.71s |  0.76∙ |  0.71s
wgs_her2+tnbc | 286 | 56 |  0.99s∘ |  0.99∘ |  0.82s
wgs_her2+her2- | 742 | 56 |  0.99s∘ |  0.99∘ |  0.83s
wgs_her2-tnbc | 780 | 56 |  0.90s∘ |  0.90∘ |  0.76s

  * $\bullet$

result significantly better than both competitors

  * $\circ$

result significantly better than one of the competitors

  * $b$

RF used a bag of words representation of the data

  * $m$

maximum performance obtained by SF out of all available distance measures

  * $s$

RF and SF used data with distribution features encoded as sums

The results clearly demonstrate that all three approaches are capable of
producing high quality results on scalar data. Nevertheless, there are some
differences between the reported performance of the approaches. For most
datasets the results are very close, and either of the approaches (RF, SF,
RSF) could be used. However, on datasets splice and svmguide3 RSF produces
very similar results to RF, but both approaches significantly outperform SF by
a large margin. This outcome is probably due to the fact that although both
RSF and SF are distance-based, the former still treats each feature
separately, while the latter calculates distances over the whole feature space
which, in this case, clearly obfuscates the decision boundary. The
significantly worse performance of SF on the splice and svmguide3 datasets
showcases that it should be used with more caution than RF and RSF, which
confirms recent findings of Czekalski and Morzy [4].

###  4.5 Evaluation on Complex Data

The second experiment was designed to assess the performance of Random
Similarity Forests on complex data. This type of data should be a natural
environment for Similarity Forests, so analogously to the first experiment,
the aim of this comparison is to verify if RSF can match SF on this kind of
data. Importantly, this is also the first reported test of Similarity Forests
on actual complex data without explicitly defined features, since the
algorithm was only tested on scalar data in the original paper [15] and the
follow-up study [4].

In this group of experiments, both SF and RSF rely on the same set of distance
measures for projections. For sequences of sets the algorithms use edit
distance (with Jaccard distance as the cost function for element relabeling),
for time series they use dynamic time warping, whereas for graphs the Ipsen-
Mikhailov distance is used. We also wanted to include RF in this comparison,
however, RF is incapable of processing varied sized sequences and graphs in
their original form. Therefore, for this classifier, we transform sequences of
sets and graphs to a bag-of-words representation, where each sequence of
sets/graph is represented as counts of items/nodes. Since time series were of
the same length for a given dataset, RF could use individual time points as
features, without any preprocessing. The results of this set of experiments
are presented in Table 1.

We can observe that for the items dataset a simple bag-of-words representation
is able to easily distinguish between the two classes, with RF significantly
outperforming SF and RSF. This result clearly demonstrates that for easy
problems, relying on approaches based on complex distance measures may
actually be detrimental to the quality of predictions. In this case, the
sequential information only obfuscates the true source of class information.
On the length dataset, where examples in the two classes are differentiated
only by the lengths of the sequences, RSF and SF both significantly
outperformed RF. This is expected as the information about the length of the
sequences is very easy to capture through edit distance. Finally, on the order
dataset, where the elements in one of the classes are ordered while in the
other they are not, RF was barely able to find any regularities and was
significantly outperformed by both SF and RSF.

Results for real-world time series and graphs confirm the above observations.
On each time series and graph dataset RSF and SF significantly outperformed
RF. The differences between RSF and SF were much smaller and, with the
exception of the proteins dataset, were not statistically significant.

###  4.6 Evaluation on Mixed Data

The final set of experiments was designed to evaluate Random Similarity
Forests on datasets with a mix of simple numerical features and complex
object-like features. These are the types of datasets that neither Random
Forest nor Similarity Forest can process out of the box. Here, RF used only
the numerical features, whereas SF treated the examples as complex objects. In
cases where multiple distance measures were available for the data (time
series, graphs, multiomics), RSF used each measure as an independent feature
whereas SF used each of them individually. Therefore, for each dataset with
multiple distances, SF was ran multiple times, each time with a different
measure. The best performing distance was later taken into account during
statistical tests.

The results clearly show that RSF is the best choice for mixed data. For 7 out
of 10 datasets RSF achieved the best performance and was never statistically
significantly worse than any other approach. It is worth noting that on the
bzr and wgs_ovarian datasets, RSF outperforms RF and SF by a large margin.
These are two stark examples that on some datasets the combination of numeric
and complex data can offer much more than any of the component representations
alone, and that multiomics data are potentially a good training ground for
algorithms like RSF. Interestingly, there is no clear runner-up in this set of
experiments. On the beetlefly, birdchicken, rosmap, and wgs_ datasets SF is
the significantly underperforming algorithm, whereas on bzr, dhfr, brca RF
takes last spot. This suggests that depending on the application domain,
numeric or complex features play a more important role. Algorithms like RSF
can take advantage of this property and steer towards the more rewarding
feature types for a given dataset.

###  4.7 Summary of the Results

The main takeaways from the experiments are as follows.

  * •

Due to the increased randomness, Random Similarity Forests usually need to
evaluate more features (max_features) during splits than Random Forests.

  * •

The proposed algorithm was able to match the performance of Random Forests and
match or outperform Similarity Forests on datasets consisting of scalar
features.

  * •

Random Similarity Forest was also able to match or outperform Similarity
Forests and outperform Random Forests on complex data.

  * •

Random Similarity Forest was the only algorithm capable of using a mixture of
scalar and complex features. This property helped Random Similarity Forests
outperform Random Forests and Similarity Forests on datasets with both types
of features. The experiments have shown that for real-world omics data
transforming complex features to scalars or using a single example-wide
distance measure may result in information loss.

##  5 Conclusions

In this paper, we have addressed the problem of classifying data with a
mixture of scalar and complex features by proposing a new classifier, called
Random Similarity Forest. Like other decision forest algorithms, such as
Random Forests and Similarity Forests, the proposed algorithm relies on
bagging of decision trees, but in contrast to other approaches splits the tree
nodes according to distance-based projections of single feature values. This
allows for a very flexible classification approach, in which traditional
numerical features can be used alongside complex object-like features with
domain-specific distance measures, taking full advantage of each data type.
The proposed method was experimentally evaluated against Random Forests and
Similarity Forests on datasets with scalar, complex, and mixed data features,
showcasing its capabilities across all possible data domain configurations.

Experimental results on 30 datasets show that Random Similarity Forests can be
considered a safe alternative to the well-established Random Forests on scalar
datasets and are on par or exceed the performance of Similarity Forests on
complex data. Most importantly, however, the experiments show that the
proposed classifier was the only one inherently capable of dealing with
datasets consisting of a mixture of scalar and complex features, and
outperformed both of the alternatives on most datasets with mixed feature
types. In situations were numeric and distance features offer complementary
information, Random Similarity Forests can achieve results that are not
attainable by the individual feature types alone. In situations were only
numeric or only complex features are informative, the proposed algorithm
achieves results that are close to the algorithm using only the informative
representation.

The method proposed in this paper opens many new avenues for further research.
Regarding the classifier itself, it would be interesting to verify the impact
of different proportions of numeric and complex features on classifier
performance. To control this effect, one could implement a weighting mechanism
that influences how many features of each type are taken into account at a
given tree node split. Moreover, the effects of other forms of randomization
on the proposed algorithm, such as random splits in the projection space,
could also be inspected. Finally, it will be interesting to observe what types
of features and distance measures are most useful when applied to different
types of multiomics data, which seem to be one of the most promisisng
application domains for Random Similarity Forests.

## References

  * [1] Davies et al., H.: HRDetect is a predictor of BRCA1 and BRCA2 deficiency based on mutational signatures. Nature Medicine 23(4), 517–525 (2017)
  * [2] Breiman, L., Friedman, J.H., Olshen, R.A., Stone, C.J.: Classification and Regression Trees. Wadsworth International Group (1984)
  * [3] Breiman, L.: Random forests. Machine Learning 45(1), 5–32 (2001)
  * [4] Czekalski, S., Morzy, M.: Similarity forests revisited: A swiss army knife for machine learning. In: Advances in Knowledge Discovery and Data Mining. pp. 42–53 (2021)
  * [5] Demsar, J.: Statistical comparisons of classifiers over multiple data sets. Journal of Machine Learning Research 7, 1–30 (2006)
  * [6] Fernández-Delgado, M., Cernadas, E., Barro, S., Amorim, D.: Do we need hundreds of classifiers to solve real world classification problems? J. Mach. Learn. Res. 15(1), 3133–3181 (2014)
  * [7] Geurts, P., Ernst, D., Wehenkel, L.: Extremely randomized trees. Machine Learning 63(1), 3–42 (2006)
  * [8] Huang, J., Ling, C.X.: Using AUC and accuracy in evaluating learning algorithms. IEEE Transactions on Knowledge and Data Engineering 17(3), 299–310 (2005)
  * [9] LeCun, Y., Bengio, Y., Hinton, G.E.: Deep learning. Nat. 521(7553), 436–444 (2015)
  * [10] Liang, J., Liu, Q., Nie, N., Zeng, B., Zhang, Z.: An improved algorithm based on knn and random forest. In: Proceedings of the 3rd International Conference on Computer Science and Application Engineering (2019)
  * [11] McCabe, S., Torres, L., LaRock, T., Haque, S.A., Yang, C.H., Hartle, H., Klein, B.: netrd: A library for network reconstruction and graph distances. Journal of Open Source Software 6(62), 2990 (2021)
  * [12] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E.: Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825–2830 (2011)
  * [13] Piernik, M., Morzy, T.: A study on using data clustering for feature extraction to improve the quality of classification. Knowledge and Information Systems (2021)
  * [14] Reel, P.S., Reel, S., Pearson, E., Trucco, E., Jefferson, E.: Using machine learning approaches for multi-omics data analysis: A review. Biotechnology Advances 49, 107739 (2021)
  * [15] Sathe, S., Aggarwal, C.C.: Similarity forests. In: Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. p. 395–403 (2017)
  * [16] Scornet, E.: Random forests and kernel methods. IEEE Transactions on Information Theory 62(3), 1485–1500 (2016)
  * [17] Tsai, C.F., Lin, W.Y., Hong, Z.F., Hsieh, C.Y.: Distance-based features in pattern classification. EURASIP Journal on Advances in Signal Processing 2011(1), 62 (2011)