# Scalable and Provably Fair Exposure Control
for Large-Scale Recommender Systems

Riku Togashi  CyberAgentTokyoJapan<EMAIL_ADDRESS>,  Kenshi Abe
CyberAgentTokyoJapan<EMAIL_ADDRESS>and  Yuta Saito  Cornell
UniversityNew YorkUnited States<EMAIL_ADDRESS>
(2024)

###### Abstract.

Typical recommendation and ranking methods aim to optimize the satisfaction of
users, but they are often oblivious to their impact on the items (e.g.,
products, jobs, news, video) and their providers. However, there has been a
growing understanding that the latter is crucial to consider for a wide range
of applications, since it determines the utility of those being recommended.
Prior approaches to fairness-aware recommendation optimize a regularized
objective to balance user satisfaction and item fairness based on some notion
such as exposure fairness. These existing methods have been shown to be
effective in controlling fairness, however, most of them are computationally
inefficient, limiting their applications to only unrealistically small-scale
situations. This indeed implies that the literature does not yet provide a
solution to enable a flexible control of exposure in the industry-scale
recommender systems where millions of users and items exist. To enable a
computationally efficient exposure control even for such large-scale systems,
this work develops a scalable, fast, and fair method called _ex posure-aware
ADMM (exADMM)_. exADMM is based on implicit alternating least squares (iALS),
a conventional scalable algorithm for collaborative filtering, but optimizes a
regularized objective to achieve a flexible control of accuracy-fairness
tradeoff. A particular technical challenge in developing exADMM is the fact
that the fairness regularizer destroys the separability of optimization
subproblems for users and items, which is an essential property to ensure the
scalability of iALS. Therefore, we develop a set of optimization tools to
enable yet scalable fairness control with provable convergence guarantees as a
basis of our algorithm. Extensive experiments performed on three
recommendation datasets demonstrate that exADMM enables a far more flexible
fairness control than the vanilla version of iALS, while being much more
computationally efficient than existing fairness-aware recommendation methods.

Recommender systems, Fairness, Collaborative filtering, Scalability

††journalyear: 2024††copyright: rightsretained††conference: Proceedings of the
ACM Web Conference 2024; May 13–17, 2024; Singapore, Singapore††booktitle:
Proceedings of the ACM Web Conference 2024 (WWW ’24), May 13–17, 2024,
Singapore, Singapore††doi: 10.1145/3589334.3645390††isbn:
979-8-4007-0171-9/24/05††ccs: Information systems Collaborative filtering

##  1\. Introduction

Personalized recommender system has been a core function of many online
platforms such as e-commerce, advertising, dating app, and online job markets.
In these systems, the items to be recommended and ranked are products, job
candidates, or other entities that transfer economic benefit, and it is widely
recognized that how they are exposed to users has a crucial influence on their
economic success (Morik et al., 2020; Saito and Joachims, 2022; Singh and
Joachims, 2019). It has also been recognized that recommender systems are
responsible for and should be aware of potential societal concerns in diverse
contexts, such as popularity bias (Park and Tuzhilin, 2008; Abdollahpouri et
al., 2017, 2019; Saito et al., 2020; Saito, 2020), sales concentration in
e-commerce (Fleder and Hosanagar, 2009; Brynjolfsson et al., 2011), filter
bubbles, biased news recommendation in social media sites (Resnick et al.,
2013; Helberger et al., 2018), and item-side fairness in two-sided markets
(Burke, 2017; Abdollahpouri et al., 2020). In essence, these concerns all
demand a form of _exposure control_ to ensure that each item receives “fair”
exposure to relevant users while not greatly sacrificing user satisfaction.
However, implementing exposure control is technically challenging,
particularly in large-scale systems, due to the complex and often intractable
nature of its objective function. This work therefore aims to develop an
effective and scalable approach to control the tradeoff between exposure
allocation and recommendation accuracy.

##### Related Work

In the context of fairness-aware recommendation and ranking, there exist
numerous studies on learning fair probabilistic rankings based on _pre-
trained_ preferences (Biega et al., 2018; Singh and Joachims, 2018; Memarrast
et al., 2021; Do et al., 2021; Do and Usunier, 2022; Saito and Joachims,
2022). The problem is often formulated as a convex optimization on doubly
stochastic matrices with the size of $|\mathcal{V}|\times|\mathcal{V}|$ (where
$\mathcal{V}$ is the item set) for each user. Whereas formulating ranking
optimization through doubly stochastic matrices is advantageous for
differentiability and convexity, this approach may not be applicable to most
industry systems because of the space complexity of
$\mathcal{O}(|\mathcal{U}||\mathcal{V}|K)$ for top-$K$ recommendation (where
$\mathcal{U}$ is the user set). More recent methods (Do et al., 2021; Usunier
et al., 2022; Do and Usunier, 2022) are based on the Frank–Wolfe-type
efficient algorithm (Frank and Wolfe, 1956; Jaggi, 2013), which requires a
top-$K$ sorting of items for each user at each iteration, resulting in a
computational cost of $\mathcal{O}(|\mathcal{U}||\mathcal{V}|\log K)$ per
training epoch, which is still prohibitively high. Patro et al. (2020)
proposed a greedy round-robin algorithm called FairRec, which also does not
scale well because its round-robin scheduling is not parallelizable. Notably,
these _post-processing_ methods require, a priori, a
$|\mathcal{U}|\times|\mathcal{V}|$ (dense) preference matrix (e.g., a
preference matrix estimated by an matrix factorization, MF, model). The
preference matrix is costly to retain in the memory space or even impossible
to materialize due to its cost of $\mathcal{O}(|\mathcal{U}||\mathcal{V}|)$.
Therefore, these post-processing approaches cannot exploit feedback sparsity,
leading to the time and space costs of
$\mathcal{O}(|\mathcal{U}||\mathcal{V}|)$. Note that a similar approach is
also adopted to control popularity bias (Steck, 2018; Abdollahpouri et al.,
2019). In particular, Abdollahpouri et al. (2019) proposed a re-ranking
algorithm based on xQuAD (Santos et al., 2010), which also suffers from a
quadratic computational cost to the number of items for each user, and thus it
is infeasible for large-scale systems.

In contrast to the post-processing approach, various studies have explored an
_in-processing_ counterpart where a single model is trained to jointly
optimize recommendation accuracy and exposure fairness (Kamishima et al.,
2011, 2013; Yao and Huang, 2017; Abdollahpouri et al., 2017; Burke et al.,
2018; Singh and Joachims, 2019; Zehlike and Castillo, 2020; Morik et al.,
2020; Yadav et al., 2021; Oosterhuis, 2021; Wu et al., 2022). To represent a
stochastic ranking policy, several studies (Singh and Joachims, 2019; Yadav et
al., 2021; Oosterhuis, 2021) rely on the Plackett–Luce (PL) model (Plackett,
1975), which has a cost of $\mathcal{O}(|\mathcal{U}||\mathcal{V}|K)$ per
training epoch. When optimizing the PL model and the joint objective of
ranking quality and exposure equality, stochastic gradient descent (SGD) is
often applied. Although SGD allows flexible objectives and reduces the
computational cost per training step, it is often difficult to apply in
practice due to its severely slow convergence, particularly when the item
catalogue is large (Yu et al., 2014a; Chen et al., 2023).

Compared to often inefficient fairness methods, iALS is a conventional
algorithm to enable scalable recommender systems. It has also been shown by
Rendle et al. (2020, 2021) to perform more effectively in terms of
recommendation accuracy than neural collaborative filtering (NCF) (He et al.,
2017) with a proper hyperparameter tuning. Rendle et al. (2021) reported that
MF models with the ALS solver (Hu et al., 2008) can be further improved by
using a customized Tikhonov regularization, which is an extension of the well-
known technique proposed by Zhou et al. (2008). This latest version of iALS
(Rendle et al., 2021) shows competitive accuracy with that of state-of-the-art
methods, such as Mult-VAE (Liang et al., 2018), while substantially improving
computational cost and scalability. This empirical evidence motivates us to
extend iALS to build a first recommendation method that is fairness-aware and
scalable to large-scale systems with over million users and items.

##### Our Contributions

To enable a scalable and provably fair exposure control for large-scale
recommender systems, this work develops a new recommendation algorithm called
_exposure-aware ADMM (exADMM)_ , which is an extension of the celebrated iALS
algorithm to achieve a scalable and flexible exposure-control. This extension
is novel and non-trivial, since any fairness regularizer in the objective
function introduces dependency between all users and items. In particular,
this intrinsic dependency inevitably destroys optimization separability, which
is a crucial property that ensures the scalability of iALS. To overcome this
technical difficulty in building a scalable method to control item fairness,
we develop a set of novel optimization tools based on the alternating
direction method of multipliers (ADMM) (Boyd et al., 2011). Furthermore, we
provide a convergence guarantee for the proposed algorithm in terms of a non-
trivial objective that includes a fairness regularizer, despite the non-convex
and multi-block optimization. Finally, we provide a comprehensive empirical
analysis on three datasets and demonstrate that exADMM outperforms the vanilla
version of iALS in terms of fairness control while maintaining its scalability
and computational efficiency. In addition, exADMM achieves similar
effectiveness in terms of accuracy-fairness tradeoff compared to typical fair
recommendation methods while being much more scalable and computationally
efficient.

Our contributions can be summarized as follows.

  * •

We propose a first scalable method (exADMM) to enable a flexible control of
accuracy-fairness tradeoff for large-scale recommender systems with over
million users and items.

  * •

We develop a set of optimization tools based on ADMM to enable an extension of
iALS to exADMM with provable convergence guarantees.

  * •

We empirically demonstrate that exADMM achieves similar scalability to iALS
and similarly effective accuracy-fairness tradeoff compared to existing
(computationally inefficient) fairness methods.

##  2\. Problem Formulation and iALS

This section formulates the typical recommendation problem and the core
technical details of iALS as a basis of our method.

Given users $\mathcal{U}\coloneqq[|\mathcal{U}|]$ and items
$\mathcal{V}\coloneqq[|\mathcal{V}|]$, let
$\mathbf{R}\in\\{0,1\\}^{|\mathcal{U}|\times|\mathcal{V}|}$ be an implicit
feedback matrix whose $(i,j)$-element has the value of $1$ when user
$i\in\mathcal{U}$ has interacted with item $j\in\mathcal{V}$; otherwise, it
has a value of $0$. We represent the number of observed interactions by that
of non-zero entries in $\mathbf{R}$, which is denoted as
$\mathrm{nz}(\mathbf{R})$. iALS is an MF-based method, and its model
parameters are $d$-dimensional embeddings,
$\mathbf{U}\in\mathbb{R}^{|\mathcal{U}|\times d}$ and
$\mathbf{V}\in\mathbb{R}^{|\mathcal{V}|\times d}$ for users and items,
respectively. These parameters are typically learned by minimizing the
following objective:

$\displaystyle L(\mathbf{V},\mathbf{U})$
$\displaystyle\coloneqq\frac{1}{2}\left\lVert\mathbf{R}\odot(\mathbf{R}-\mathbf{U}\mathbf{V}^{\top})\right\rVert_{F}^{2}+\frac{\alpha_{0}}{2}\left\lVert\mathbf{U}\mathbf{V}^{\top}\right\rVert_{F}^{2}$
(1)
$\displaystyle+\frac{1}{2}\left\lVert\mathbf{\Lambda}_{U}^{1/2}\mathbf{U}\right\rVert_{F}^{2}+\frac{1}{2}\left\lVert\mathbf{\Lambda}_{V}^{1/2}\mathbf{V}\right\rVert_{F}^{2},$

where operator $\odot$ is the Hadamard element-wise product, and the second
term is the L2 norm of the recovered score matrix
$\mathbf{U}\mathbf{V}^{\top}$ (i.e., implicit regularizer (Bayer et al.,
2017)) with a weight parameter $\alpha_{0}>0$.
$\mathbf{\Lambda}_{U}\in\mathbb{R}^{|\mathcal{U}|\times|\mathcal{U}|}$ and
$\mathbf{\Lambda}_{V}\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{V}|}$ are
diagonal matrices (a.k.a. Tikhonov matrices (Zhou et al., 2008)) representing
the weights for L2 regularization. Let $\mathbf{r}_{i,\cdot}$ and
$\mathbf{r}_{\cdot,j}$ be the (column) vectors that correspond to the $i$-th
row and the $j$-th column of $\mathbf{R}$, respectively. The frequency-based
weighting strategy sets the weights with internal hyperparameters
$\lambda_{L2}>0$ and exponent $\eta\geq 0$ as follows:

$\displaystyle(\mathbf{\Lambda}_{U})_{i,i}\coloneqq\lambda_{L2}\left(\left\lVert\mathbf{r}_{i,\cdot}\right\rVert_{1}\\!+\\!\alpha_{0}|\mathcal{V}|\right)^{\eta},\,\,(\mathbf{\Lambda}_{V})_{j,j}\coloneqq\lambda_{L2}\left(\left\lVert\mathbf{r}_{\cdot,j}\right\rVert_{1}\\!+\\!\alpha_{0}|\mathcal{U}|\right)^{\eta}.$

Hereafter, we use $\lambda_{U}^{(i)}\coloneqq(\mathbf{\Lambda}_{U})_{i,i}$ and
$\lambda_{V}^{(j)}\coloneqq(\mathbf{\Lambda}_{V})_{j,j}$.

iALS solves the minimization problem in Eq. 1 by alternating the optimization
of $\mathbf{V}$ and $\mathbf{U}$. Specifically, in the $k$-th step, iALS
updates $\mathbf{U}$ and $\mathbf{V}$ via

$\displaystyle\mathbf{U}^{k+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{U}}\|\mathbf{R}\odot(\mathbf{R}-\mathbf{U}(\mathbf{V}^{k})^{\top})\|_{F}^{2}+\alpha_{0}\|\mathbf{U}(\mathbf{V}^{k})^{\top}\|_{F}^{2}+\|\mathbf{\Lambda}_{U}^{1/2}\mathbf{U}\|_{F}^{2},$
$\displaystyle\mathbf{V}^{k+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{V}}\|\mathbf{R}\odot(\mathbf{R}-{\mathbf{U}^{k+1}}\mathbf{V}^{\top})\|_{F}^{2}+\alpha_{0}\|{\mathbf{U}^{k+1}}\mathbf{V}^{\top}\|_{F}^{2}+\|\mathbf{\Lambda}_{V}^{1/2}\mathbf{V}\|_{F}^{2}.$

Owing to the alternating strategy, the optimization of $\mathbf{U}$ and
$\mathbf{V}$ can be divided into independent convex problems for each row of
$\mathbf{U}$ and $\mathbf{V}$. Let us use $\mathbf{u}_{i}\in\mathbb{R}^{d}$ to
denote the (column) vector that corresponds to the $i$-th row of $\mathbf{U}$.
Then, its update can simply be done via the following row-wise independent
problem:

$\displaystyle\mathbf{u}_{i}^{k+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{u}_{i}}\left\lVert\mathbf{r}_{i}\odot(\mathbf{r}_{i}-\mathbf{V}^{k}\mathbf{u}_{i})\right\rVert_{2}^{2}+\alpha_{0}\left\lVert\mathbf{V}^{k}\mathbf{u}_{i}\right\rVert_{2}^{2}+\lambda_{U}^{(i)}\left\lVert\mathbf{u}_{i}\right\rVert_{2}^{2}$
$\displaystyle=\left(\sum_{j\in\mathcal{V}}r_{i,j}\mathbf{v}_{j}^{k}(\mathbf{v}_{j}^{k})^{\top}+\alpha_{0}\mathbf{G}_{V}^{k}+\lambda_{U}^{(i)}\mathbf{I}\right)^{-1}\sum_{j\in\mathcal{V}}r_{i,j}\mathbf{v}_{j}^{k},$

where
$\mathbf{G}_{V}^{k}\coloneqq\sum_{j\in\mathcal{V}}\mathbf{v}_{j}^{k}({\mathbf{v}_{j}}^{k})^{\top}=(\mathbf{V}^{k})^{\top}\mathbf{V}^{k}$
is the Gramian of the item embeddings in the $k$-th step, where
$\mathbf{v}_{j}^{k}\in\mathbb{R}^{d}$ denotes the column vector that
corresponds to the $j$-th row of $\mathbf{V}^{k}$. When $\mathbf{G}_{V}^{k}$
is pre-computed, the expected computational cost for each subproblem is
reduced to $\mathcal{O}((\mathrm{nz}(\mathbf{R})/|\mathcal{V}|)d^{2}+d^{3})$
(a.k.a. the Gramian trick (Rendle et al., 2021)), which consists of (i)
computation of the Gramian for interacted items
$\sum_{j\in\mathcal{V}}r_{i,j}(\mathbf{v}_{j}^{k})(\mathbf{v}_{j}^{k})^{\top}$
in $\mathcal{O}((\mathrm{nz}(\mathbf{R})/|\mathcal{V}|)d^{2})$ and (ii)
solving a linear system
$\mathbf{H}^{k}\mathbf{u}_{i}^{k+1}=\sum_{i\in\mathcal{U}}r_{i,j}\mathbf{v}_{j}^{k}$,
where
$\mathbf{H}^{k}\coloneqq\sum_{j\in\mathcal{V}}r_{i,j}\mathbf{v}_{j}^{k}(\mathbf{v}_{j}^{k})^{\top}+\alpha_{0}\mathbf{G}_{V}^{k}+\lambda_{U}^{(i)}\mathbf{I}$
in $\mathcal{O}(d^{3})$. Because the update of $\mathbf{V}$ is analogous to
that of $\mathbf{U}$, the total cost of updating $\mathbf{U}$ and $\mathbf{V}$
is
$\mathcal{O}(\mathrm{nz}(\mathbf{R})d^{2}+(|\mathcal{U}|+|\mathcal{V}|)d^{3})$.
This is much lower than
$\mathcal{O}(|\mathcal{U}||\mathcal{V}|d^{2}+(|\mathcal{U}|+|\mathcal{V}|)d^{3})$
due to feedback sparsity, i.e.,
$\mathrm{nz}(\mathbf{R})\ll|\mathcal{U}||\mathcal{V}|$. In summary, iALS
retains scalability, despite its objective involving _all_ user-item pairs.
The crux is that iALS exploits the Gramian trick and feedback sparsity to
avoid the intractable factor $\mathcal{O}(|\mathcal{U}||\mathcal{V}|)$.

##  3\. The exADMM Algorithm

This section develops our proposed algorithm, exADMM, which is an extension of
iALS to enable a scalable exposure control.

###  3.1. A Regularized Objective for Fairness

The aim of this paper is to enable a scalable control of item exposure so that
individual items can receive attention from users more fairly while not
sacrificing recommendation accuracy much. To this end, we consider minimizing
the following regularized objective.

(2)
$\displaystyle\min_{\mathbf{V},\mathbf{U}}\underbrace{L(\mathbf{V},\mathbf{U})}_{\textit{typical
prediction
loss}}+\quad\lambda_{ex}\cdot\underbrace{R_{ex}(\mathbf{V},\mathbf{U})}_{\mathclap{\textit{(item)
fairness regularizer}}},$

where $R_{ex}(\mathbf{V},\mathbf{U})$ is a penalty term to induce exposure
equality, and $\lambda_{ex}$ is the weight hyperparameter to control the
balance between recommendation quality and exposure equality. As already
discussed, the scalability of iALS is due to the simplicity of its objective.
To retain this desirable property, we need to carefully define an exposure
regularizer $R_{ex}(\mathbf{V},\mathbf{U})$ in a still tractable way.

To define our regularizer, let us denote the predicted score for user $i$ and
item $j$ by $\hat{r}_{i,j}=(\mathbf{U}\mathbf{V}^{\top})_{i,j}$, and we
predict the item ranking for user $i$ according to the decreasing order of
$\\{\hat{r}_{i,j}\mid j\in\mathcal{V}\\}$. Evidently, there is a monotonic
relationship between $\hat{r}_{i,j}$ and the amount of exposure that item $j$
will receive in a ranked list. Hence, we can evaluate exposure inequality
under a recommendation model by the variability of items’ scores averaged over
the users, i.e., $\frac{1}{|\mathcal{U}|}\sum_{i\in\mathcal{U}}\hat{r}_{i,j}$
(this is for item $j$). There exist several possible measures of variability
such as Gini indices (Atkinson et al., 1970; Do and Usunier, 2022), standard
deviation (Do et al., 2021), and variance (Wu et al., 2021). In this work, we
consider the following second moment of the predicted scores as the
regularizer $R_{ex}$.

$\displaystyle R_{ex}(\mathbf{V},\mathbf{U})$
$\displaystyle\coloneqq\frac{1}{|\mathcal{V}|}\sum_{j\in\mathcal{V}}\left(\frac{1}{|\mathcal{U}|}\sum_{i\in\mathcal{U}}\hat{r}_{i,j}\right)^{2}=\frac{1}{|\mathcal{V}|}\left\lVert\frac{1}{|\mathcal{U}|}\mathbf{V}\mathbf{U}^{\top}\mathds{1}\right\rVert_{2}^{2},$

where $\mathds{1}$ is the $|\mathcal{U}|\times 1$ column vector of which the
elements are all $1$. The fairness regularizer defined above is the L2 norm of
the average scores predicted for the items. This is considered one of the
reasonable measures of exposure inequality since it takes a large value for
items whose average scores are either extremely large or small. Moreover, we
can draw a clear technical distinction between our fairness regularizer
$R_{ex}(\mathbf{V},\mathbf{U})$ and implicit regularizer
$\|\mathbf{V}\mathbf{U}^{\top}\|_{F}^{2}$ of the vanilla iALS in Eq. 1. That
is, the implicit regularizer penalizes the score
$(\mathbf{U}\mathbf{V}^{\top})_{i,j}$ of each user-item pair
$(i,j)\in\mathcal{U}\times\mathcal{V}$ _independently_ , whereas our penalty
term introduces a _structural dependency_ into the recovered matrix
$\mathbf{U}\mathbf{V}^{\top}$. Unfortunately, this structural dependency
destroys the optimization separability with respect to the rows of
$\mathbf{U}$ due to the averaged user embedding
$(1/|\mathcal{U}|)\mathbf{U}^{\top}\mathds{1}$ that appears in its definition.
Optimizing $R_{ex}$ is thus not straightforward, particularly in large-scale
settings, which motivates us to develop novel tools to handle this fairness
regularizer in a scalable and provable fashion.

###  3.2. Scalable Optimization based on ADMM

To enable parallel optimization of our exposure-controllable objective in Eq.
2, we adopt an approach based on ADMM, which is an optimization framework with
high parallelism (Boyd et al., 2011) and has been adopted to enable scalable
recommendations (Yu et al., 2014b; Cheng et al., 2014; Smith et al., 2017;
Ioannidis et al., 2019; Steck et al., 2020; Steck and Liang, 2021). To
decouple the row- and column-wise dependencies in $\mathbf{U}$ introduced by
our fairness regularizer, we first reformulate the optimization problem by
introducing an auxiliary variable $\mathbf{s}\in\mathbb{R}^{d}$ as

(3)
$\displaystyle\min_{\mathbf{V},\mathbf{U},\mathbf{s}}\;L(\mathbf{V},\mathbf{U})+\frac{\lambda_{ex}}{2}\left\lVert\mathbf{V}\mathbf{s}\right\rVert_{2}^{2},\quad\text{s.t.~{}}\mathbf{s}=\frac{1}{|\mathcal{U}|}\mathbf{U}^{\top}\mathds{1}.$

Here, we replaced $(1/|\mathcal{U}|)\mathbf{U}^{\top}\mathds{1}$ in the
fairness regularizer with $\mathbf{s}$ while introducing an additional linear
equality constraint. This can be further reformulated to the following saddle-
point optimization:

$\displaystyle\min_{\mathbf{V},\mathbf{U},\mathbf{s}}\max_{\mathbf{w}}L_{\rho}(\mathbf{V},\mathbf{U},\mathbf{s},\mathbf{w}),$

where

$\displaystyle L_{\rho}(\mathbf{V},\mathbf{U},\mathbf{s},\mathbf{w})$
$\displaystyle\coloneqq
L(\mathbf{V},\mathbf{U})+\frac{\lambda_{ex}}{2}\left\lVert\mathbf{V}\mathbf{s}\right\rVert_{2}^{2}+\frac{\rho}{2}\left\lVert\frac{1}{|\mathcal{U}|}\mathbf{U}^{\top}\mathds{1}-\mathbf{s}+\mathbf{w}\right\rVert_{2}^{2}-\frac{\rho}{2}\|\mathbf{w}\|_{2}^{2}.$

$L_{\rho}$ is the Lagrangian augmented by the ADMM penalty term with weight
$\rho>0$, and $\mathbf{w}\in\mathbb{R}^{d}$ is the dual variable (i.e.,
Lagrange multipliers) scaled by $1/\rho$. We can perform optimization in the
$(k+1)$-th step by iteratively updating each variable as

$\displaystyle\mathbf{V}^{k+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{V}}L_{\rho}(\mathbf{V},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k}),$
$\displaystyle\mathbf{U}^{k+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{U}}L_{\rho}(\mathbf{V}^{k+1},\mathbf{U},\mathbf{s}^{k},\mathbf{w}^{k}),$
$\displaystyle\mathbf{s}^{k+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{s}}L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s},\mathbf{w}^{k}),$
$\displaystyle\mathbf{w}^{k+1}$
$\displaystyle=\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}^{k+1},$

The update of $\mathbf{w}$ corresponds to the gradient ascent with respect to
the dual problem
$\max_{\mathbf{w}}\min_{\mathbf{V},\mathbf{U},\mathbf{s}}L_{\rho}(\mathbf{V},\mathbf{U},\mathbf{s},\mathbf{w})$
with step size $\rho$ (Boyd et al., 2011).

####  3.2.1. Update of $\mathbf{V}$.

Next, we derive how to update $\mathbf{V}$ in the $(k+1)$-th step, which
comprises independent optimization problems for the rows of $\mathbf{V}$.
Suppose that $\mathbf{v}_{j}^{k+1}\in\mathbb{R}^{d}$ and
$\mathbf{r}_{\cdot,j}\in\\{0,1\\}^{|\mathcal{U}|}$ are the column vectors
indicating the $j$-th row of $\mathbf{V}^{k+1}$ and the $j$-th column of
$\mathbf{R}$, respectively. The update can then be performed by solving the
following linear system.

$\displaystyle\mathbf{v}_{j}^{k+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{v}_{j}}\Bigl{\\{}\frac{1}{2}\left\lVert\mathbf{r}_{\cdot,j}\odot(\mathbf{r}_{\cdot,j}-\mathbf{U}^{k}\mathbf{v}_{j})\right\rVert_{2}^{2}+\frac{\alpha_{0}}{2}\left\lVert\mathbf{U}^{k}\mathbf{v}_{j}\right\rVert_{2}^{2}$
$\displaystyle\qquad\phantom{\operatorname*{argmin}_{\mathbf{v}_{j}}\Bigl{\\{}}+\frac{\lambda_{V}^{(j)}}{2}\left\lVert\mathbf{v}_{j}\right\rVert_{2}^{2}+\frac{\lambda_{ex}}{2}\left(\mathbf{v}_{j}^{\top}\mathbf{s}^{k}\right)^{2}\Bigr{\\}}$
$\displaystyle=\left(\sum_{i\in\mathcal{U}}r_{i,j}\mathbf{u}_{i}^{k}({\mathbf{u}_{i}^{k}})^{\top}\\!+\\!\alpha_{0}\mathbf{G}_{U}^{k}\\!+\\!\lambda_{ex}\mathbf{s}^{k}(\mathbf{s}^{k})^{\top}\\!+\\!\lambda_{V}^{(j)}\mathbf{I}\right)^{-1}\\!\\!\sum_{i\in\mathcal{U}}r_{i,j}\mathbf{u}_{i}^{k},$

where
$\mathbf{G}_{U}^{k}\coloneqq\sum_{i\in\mathcal{U}}\mathbf{u}_{i}^{k}{(\mathbf{u}_{i}^{k})}^{\top}$
is the Gramian of the user embeddings in the $k$-th step. Notably, we can pre-
compute $\mathbf{G}_{U}^{k}$ and $\mathbf{s}^{k}(\mathbf{s}^{k})^{\top}$, and
thus the update of $\mathbf{V}$ achieves the same complexity as that of iALS.

####  3.2.2. Update of $\mathbf{U}$.

Updating $\mathbf{U}$ is the most intricate part of our algorithm. In the
$(k+1)$-th step, our aim is to solve the following

$\displaystyle\operatorname*{argmin}_{\mathbf{U}}$
$\displaystyle\Biggl{\\{}L(\mathbf{V}^{k+1},\mathbf{U})+\frac{\lambda_{ex}}{2}\|\mathbf{V}^{k+1}\mathbf{s}^{k}\|_{2}^{2}+\frac{\rho}{2}\left\lVert\frac{1}{|\mathcal{U}|}\mathbf{U}^{\top}\mathds{1}-\mathbf{s}^{k}+\mathbf{w}^{k}\right\rVert_{2}^{2}\Biggr{\\}}.$

The issue here is that the penalty term of ADMM (the fourth term of RHS)
destroys the separability regarding the rows of $\mathbf{U}$. We resolve this
using a proximal gradient method (Rockafellar, 1976; Duchi and Singer, 2009;
Liu et al., 2019). Specifically, we consider a linear approximation (i.e., the
first-order Taylor expansion around the current estimate $\mathbf{U}^{k}$) of
the objective except for the ADMM penalization. This yields the following
objective:

$\displaystyle\mathbf{U}^{k+1}=\operatorname*{argmin}_{\mathbf{U}}$
$\displaystyle\Bigl{\\{}\langle\mathbf{U}-\mathbf{U}^{k},\nabla_{\mathbf{U}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})\rangle_{F}+\frac{1}{2\gamma}\left\lVert\mathbf{U}-\mathbf{U}^{k}\right\rVert_{F}^{2}$
$\displaystyle\qquad+\frac{\rho}{2}\left\lVert\frac{1}{|\mathcal{U}|}\mathbf{U}^{\top}\mathds{1}-\mathbf{s}^{k}+\mathbf{w}^{k}\right\rVert_{2}^{2}\Bigr{\\}},$

where $g(\mathbf{V},\mathbf{U},\mathbf{s})\coloneqq
L(\mathbf{V},\mathbf{U})+(\nicefrac{{\lambda_{ex}}}{{2}})\left\lVert\mathbf{V}\mathbf{s}\right\rVert_{2}^{2}$.
Here, we introduce a regularization term
$(1/2\gamma)\lVert\mathbf{U}-\mathbf{U}^{k}\rVert_{F}^{2}$, which is referred
to as the proximal term (Rockafellar, 1976). By completing the square, the
above update can be rearranged into the following parallel and non-parallel
steps:

$\displaystyle\mathbf{U}^{k+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{U}}\frac{\rho}{2}\left\lVert\frac{1}{|\mathcal{U}|}\mathbf{U}^{\top}\mathds{1}-\mathbf{s}^{k}+\mathbf{w}^{k}\right\rVert_{2}^{2}+\frac{1}{2\gamma}\left\lVert\mathbf{U}-\left(\mathbf{U}^{k}-\gamma\nabla_{\mathbf{U}}g^{k}\right)\right\rVert_{F}^{2}$
$\displaystyle=\underbrace{\mathrm{prox}_{\gamma}^{k}}_{\text{non-
parallel}}(\underbrace{\vphantom{\mathrm{prox}_{\gamma}^{k}}\mathbf{U}^{k}-\gamma\nabla_{\mathbf{U}}g^{k}}_{\text{parallel}}),$

where

$\displaystyle\mathrm{prox}_{\gamma}^{k}(\widetilde{\mathbf{U}})$
$\displaystyle\coloneqq\operatorname*{argmin}_{\mathbf{U}}\frac{\rho}{2}\left\lVert\frac{1}{|\mathcal{U}|}\mathbf{U}^{\top}\mathds{1}-\mathbf{s}^{k}+\mathbf{w}^{k}\right\rVert_{2}^{2}+\frac{1}{2\gamma}\left\lVert\mathbf{U}-\widetilde{\mathbf{U}}\right\rVert_{F}^{2}$
$\displaystyle=\left(\frac{\rho}{|\mathcal{U}|^{2}}\mathds{1}\mathds{1}^{\top}+\frac{1}{\gamma}\mathbf{I}\right)^{-1}\left(\frac{1}{\gamma}\widetilde{\mathbf{U}}+\frac{\rho}{|\mathcal{U}|}\mathds{1}(\mathbf{s}^{k}-\mathbf{w}^{k})^{\top}\right).$

$\nabla_{\mathbf{U}}g^{k}$ is used to represent
$\nabla_{\mathbf{U}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})$ for
brevity. Note that the term $\mathbf{U}^{k}-\gamma\nabla_{\mathbf{U}}g^{k}$
corresponds to a gradient descent with respect to the iALS objective.111Note
that $\nabla_{\mathbf{U}}g(\mathbf{V},\mathbf{U},\mathbf{s})$ is equivalent to
the gradient of $L(\mathbf{V},\mathbf{U})$ with respect to $\mathbf{U}$
because we can ignore the constant exposure penalty
$(\lambda_{ex}/2)\|\mathbf{V}\mathbf{s}\|_{2}^{2}$. Therefore, we can update
$\mathbf{U}$ in two row-wise parallel and non-parallel steps, that is, (i)
gradient descent
$\widetilde{\mathbf{U}}^{k+1}=\mathbf{U}^{k}-\gamma\nabla_{\mathbf{U}}g^{k}$
and (ii) proximal mapping
$\mathbf{U}^{k+1}=\mathrm{prox}_{\gamma}^{k}(\widetilde{\mathbf{U}}^{k+1})$.

##### Parallel gradient computation

The gradient
$\nabla_{\mathbf{U}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})$ can be
independently computed for each row of $\mathbf{U}$ as follows:

$\displaystyle\nabla_{\mathbf{u}_{i}}g^{k}$
$\displaystyle=\left(\sum_{j\in\mathcal{V}}r_{i,j}{\mathbf{v}_{j}}^{k+1}{\left(\mathbf{v}_{j}^{k+1}\right)}^{\top}+\alpha_{0}\mathbf{G}_{V}^{k+1}+\lambda_{U}^{(i)}\mathbf{I}\right)\mathbf{u}_{i}^{k}-\sum_{j\in\mathcal{V}}r_{i,j}\mathbf{v}_{j}^{k+1}.$

Similar to iALS, we can efficiently compute the gradient by pre-computing the
Gramian $\mathbf{G}_{V}^{k+1}=(\mathbf{V}^{k+1})^{\top}\mathbf{V}^{k+1}$.
Thus, the gradient descent $\mathbf{U}^{k}-\nabla_{\mathbf{U}}g^{k}$ can be
performed in parallel with respect to users while maintaining efficiency by
exploiting the Gramian trick and feedback sparsity. It should be noted that we
can avoid computing the inverse Hessian in $\mathcal{O}(d^{3})$ unlike the
$\mathbf{U}$ step of iALS.

##### Efficient proximal mapping

The proximal mapping step requires a costly inversion of the
$|\mathcal{U}|\times|\mathcal{U}|$ matrix in $\mathcal{O}(|\mathcal{U}|^{3})$
for a naive computation. This is problematic because, in practice, $\rho$ and
$\gamma$ may increase/decrease during iterations (Boyd et al., 2011). However,
we can indeed compute the inverse matrix efficiently by leveraging the
Sherman-Morrison formula (Sherman and Morrison, 1950) (a special case of the
Woodbury matrix identity (Woodbury, 1950)), which yields the following

$\displaystyle\left(\frac{\rho}{|\mathcal{U}|^{2}}\mathds{1}\mathds{1}^{\top}+\frac{1}{\gamma}\mathbf{I}\right)^{-1}$
$\displaystyle=-\frac{(\gamma\mathbf{I})(\rho/|\mathcal{U}|^{2})\mathds{1}\mathds{1}^{\top}(\gamma\mathbf{I})}{1+(\rho/|\mathcal{U}|^{2})\mathds{1}^{\top}(\gamma\mathbf{I})\mathds{1}}+\gamma\mathbf{I}$
$\displaystyle=\gamma\left(-\frac{1}{|\mathcal{U}|^{2}\left(\frac{1}{|\mathcal{U}|}\\!+\\!\frac{1}{\rho\gamma}\right)}\mathds{1}\mathds{1}^{\top}+\mathbf{I}\right).$

Therefore, we can derive the proximal mapping $\mathrm{prox}_{\gamma}^{k}$ as
the following closed-form solution:

(4)
$\displaystyle\mathrm{prox}_{\gamma}^{k}(\mathbf{U})=\left(\frac{-1}{|\mathcal{U}|^{2}\left(\frac{1}{|\mathcal{U}|}\\!+\\!\frac{1}{\rho\gamma}\right)}\mathds{1}\mathds{1}^{\top}+\mathbf{I}\right)\left(\mathbf{U}\\!+\\!\frac{\rho\gamma}{|\mathcal{U}|}\mathds{1}(\mathbf{s}^{k}\\!-\\!\mathbf{w}^{k})^{\top}\right).$

The naive computation of $\mathrm{prox}_{\gamma}^{k}$ is still computationally
costly due to the multiplication of $|\mathcal{U}|\times|\mathcal{U}|$ and
$|\mathcal{U}|\times d$ matrices in $\mathcal{O}(|\mathcal{U}|^{2}d)$. Let us
here define
$\widehat{\mathbf{U}}\coloneqq\mathbf{U}+(\nicefrac{{\rho\gamma}}{{|\mathcal{U}|}})\mathds{1}(\mathbf{s}^{k}-\mathbf{w}^{k})^{\top}$
for simplicity, and we can further rewrite Eq. 4 as follows:

$\displaystyle\mathrm{prox}_{\gamma}^{k}(\mathbf{U})$
$\displaystyle=-\frac{1}{|\mathcal{U}|^{2}\left(\frac{1}{|\mathcal{U}|}+\frac{1}{\rho\gamma}\right)}\cdot\mathds{1}\left(\widehat{\mathbf{U}}^{\top}\mathds{1}\right)^{\top}+\widehat{\mathbf{U}}.$

Thus, we can perform this matrix multiplication efficiently by (i) computing
each row of $\widehat{\mathbf{U}}$ in parallel (i.e.,
$\widehat{\mathbf{u}}_{i}=\mathbf{u}_{i}+(\nicefrac{{\rho\gamma}}{{|\mathcal{U}|}})(\mathbf{s}^{k}-\mathbf{w}^{k})$),
(ii) computing the accumulated user embedding
$\mathbf{t}=\widehat{\mathbf{U}}^{\top}\mathds{1}$, and then (iii) adding
$-|\mathcal{U}|^{-2}\left(\nicefrac{{1}}{{|\mathcal{U}|}}+\nicefrac{{1}}{{\rho\gamma}}\right)^{-1}\cdot\mathbf{t}$
to each row of $\widehat{\mathbf{U}}$. Thus, the cost of the matrix-matrix
multiplication in Eq. 4 is reduced to $\mathcal{O}(|\mathcal{U}|d)$, which is
more efficient than $\mathcal{O}(|\mathcal{U}|^{2}d)$ of the naive
implementation. The computational efficiency is advantageous even when $\rho$
and $\gamma$ are fixed during optimization.

####  3.2.3. Update of $\mathbf{s}$.

We can perform the update of $\mathbf{s}$ by computing the following solution:

$\displaystyle\mathbf{s}^{k+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{s}}\left\\{\frac{\lambda_{ex}}{2}\left\lVert\mathbf{V}^{k+1}\mathbf{s}\right\rVert_{2}^{2}+\frac{\rho}{2}\left\lVert\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}+\mathbf{w}^{k}\right\rVert_{2}^{2}\right\\}$
$\displaystyle=\rho\left(\lambda_{ex}\mathbf{G}_{V}^{k+1}+\rho\mathbf{I}\right)^{-1}\left(\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}+\mathbf{w}^{k}\right).$

We can reuse the Gramian $\mathbf{G}_{V}^{k+1}$ for this step following its
pre-computation in the $\mathbf{U}$ step. The computation cost here is thus
$\mathcal{O}(|\mathcal{U}|d+d^{3})$, which includes (i) the computation of
$(1/|\mathcal{U}|)(\mathbf{U}^{k+1})^{\top}\mathds{1}$ and (ii) the solution
of a linear system of size $d^{2}$.

###  3.3. Complexity Analysis

Algorithm 1 shows the detailed implementation of exADMM. First, the user and
item embeddings are initialized with independent normal noise with a
$\sigma/\sqrt{d}$ standard deviation (Rendle et al., 2021). In line 10 of
Algorithm 1, we pre-compute
$\mathbf{G}_{V}^{k}=\sum_{j\in\mathcal{V}}\mathbf{v}_{j}^{k}(\mathbf{v}_{j}^{k})^{\top}$,
which can be done in parallel for each item and be reused in the update of
$\mathbf{U}$ and $\mathbf{s}$. The averaged user vector
$\mathbf{t}=(1/|\mathcal{U}|)(\mathbf{U})^{\top}\mathds{1}$ can be reused for
the $\mathbf{s}$ and $\mathbf{w}$ steps, hence, we compute this in line 20.
Consequently, the computational costs for updating $\mathbf{V}$, $\mathbf{U}$,
$\mathbf{s}$, and $\mathbf{w}$ are (i)
$\mathcal{O}(\mathrm{nz}(\mathbf{R})d^{2}+|\mathcal{V}|d^{3})$, (ii)
$\mathcal{O}(\mathrm{nz}(\mathbf{R})d^{2}+|\mathcal{U}|d^{2})$, (iii)
$\mathcal{O}(|\mathcal{U}|d+d^{3})$, and (iv) $\mathcal{O}(d)$, respectively.
Therefore, the overall cost is
$\mathcal{O}(\mathrm{nz}(\mathbf{R})d^{2}+|\mathcal{U}|d^{2}+|\mathcal{V}|d^{3})$,
which is indeed even lower than
$\mathcal{O}(\mathrm{nz}(\mathbf{R})d^{2}+(|\mathcal{U}|+|\mathcal{V}|)d^{3})$
of iALS. This is because we avoid solving the linear system when updating
$\mathbf{U}$ by applying the proximal gradient method with the efficient
$\mathrm{prox}_{\gamma}^{k}$.

Algorithm 1 exADMM

1:Implicit feedback matrix $\mathbf{R}$

2:$\forall
i\in\mathcal{U},\mathbf{u}_{i}^{0}\sim\mathcal{N}(0,(\sigma/\sqrt{d})\mathbf{I})$,
$\forall
j\in\mathcal{V},\mathbf{v}_{j}^{0}\sim\mathcal{N}(0,(\sigma/\sqrt{d})\mathbf{I})$,

3:$\mathbf{s}^{0}\leftarrow(1/|\mathcal{U}|)(\mathbf{U}^{0})^{\top}\mathds{1}$,
$\mathbf{w}^{0}\leftarrow\vec{0}$

4:for $k=0,\dots,T-1$ do

5:
$\mathbf{G}_{U}^{k}\leftarrow\sum_{i}\mathbf{u}_{i}^{k}{(\mathbf{u}_{i}^{k})}^{\top}$,
$\mathbf{G}_{s}^{k}\leftarrow\mathbf{s}^{k}{(\mathbf{s}^{k})}^{\top}$ //
$\mathcal{O}(|\mathcal{U}|d^{2})$ and $\mathcal{O}(d^{2})$

6: for $j=1,\dots,|\mathcal{V}|$ do // parallelizable loop

7:
$\mathbf{G}_{j}^{k}\leftarrow\sum_{i}r_{i,j}\mathbf{u}_{i}^{k}{(\mathbf{u}_{i}^{k})}^{\top}$
// $\mathcal{O}((\mathrm{nz}(\mathbf{R})/|\mathcal{V}|)d^{2})$

8:
$\mathbf{v}_{j}^{k+1}\leftarrow\left(\mathbf{G}_{j}^{k}+\alpha_{0}\mathbf{G}_{U}^{k}+\lambda_{ex}\mathbf{G}_{s}^{k}+\lambda_{V}^{(j)}\mathbf{I}\right)^{-1}\sum_{i}r_{i,j}\mathbf{u}_{i}^{k}$
// $\mathcal{O}(d^{3})$

9: end for

10:
$\mathbf{G}_{V}^{k+1}\leftarrow\sum_{j}\mathbf{v}_{j}^{k+1}{(\mathbf{v}_{j}^{k+1})}^{\top}$
// $\mathcal{O}(|\mathcal{V}|d^{2})$

11: for $i=1,\dots,|\mathcal{U}|$ do // parallelizable loop

12:
$\mathbf{G}_{i}^{k+1}\leftarrow\sum_{j}r_{i,j}\mathbf{v}_{j}^{k+1}{(\mathbf{v}_{j}^{k+1})}^{\top}$
// $\mathcal{O}((\mathrm{nz}(\mathbf{R})/|\mathcal{U}|)d^{2})$

13:
$\nabla_{\mathbf{u}_{i}}g^{k+1}\leftarrow\left(\mathbf{G}_{i}^{k+1}+\alpha_{0}\mathbf{G}_{V}^{k+1}+\lambda_{U}^{(i)}\mathbf{I}\right)\mathbf{u}_{i}^{k}-\sum_{j}r_{i,j}\mathbf{v}_{j}^{k+1}$
// $\mathcal{O}(d^{2})$

14:
$\widehat{\mathbf{u}}_{i}^{k+1}\leftarrow\mathbf{u}_{i}^{k}-\gamma\nabla_{\mathbf{u}_{i}}g^{k+1}+\frac{\rho\gamma}{|\mathcal{U}|}(\mathbf{s}^{k}-\mathbf{w}^{k})$
// $\mathcal{O}(d)$

15: end for

16: $\mathbf{t}\leftarrow\sum_{i}\widehat{\mathbf{u}}_{i}^{k+1}$ //
$\mathcal{O}(|\mathcal{U}|d)$

17: for $i=1,\dots,|\mathcal{U}|$ do // parallelizable loop

18:
$\mathbf{u}_{i}^{k+1}\leftarrow\widehat{\mathbf{u}}_{i}^{k+1}-|\mathcal{U}|^{-2}\left(\nicefrac{{1}}{{|\mathcal{U}|}}+\nicefrac{{1}}{{\rho\gamma}}\right)^{-1}\mathbf{t}$
// $\mathcal{O}(d)$

19: end for

20: $\mathbf{t}\leftarrow\frac{1}{|\mathcal{U}|}\sum_{i}\mathbf{u}_{i}^{k+1}$
// $\mathcal{O}(|\mathcal{U}|d)$

21:
$\mathbf{s}^{k+1}\leftarrow\rho\left(\lambda_{ex}\mathbf{G}_{V}^{k+1}+\rho\mathbf{I}\right)^{-1}\left(\mathbf{t}-\mathbf{w}^{k}\right)$
// $\mathcal{O}(d^{3})$

22: $\mathbf{w}^{k+1}\leftarrow\mathbf{w}^{k}+\mathbf{t}-\mathbf{s}^{k+1}$ //
$\mathcal{O}(d)$

23:end for

24:return $\mathbf{U},\mathbf{V}$

###  3.4. Convergence Analysis

The objective defined in Eq. 3 has more than two variables (i.e., three-block
optimization), which are _coupled_ (e.g., $\mathbf{U},\mathbf{V}$ in the iALS
loss function). However, multi-block ADMM does not retain a convergence
guarantee in general (Chen et al., 2016). Various algorithms have been
developed for optimization separability and provable convergence under coupled
variables (Wang et al., 2014; Deng et al., 2017; Liu et al., 2019). For
instance, Liu et al. (2019) proposed a variant of ADMM for non-convex
problems, which completely decouples variables by introducing linear
approximation when updating _all_ the coupled ones, thereby enabling parallel
gradient descent. By contrast, exADMM applies linearization only to the
$\mathbf{U}$ step and works in an alternate way. This strategy enables second-
order acceleration in the update of $\mathbf{V}$ and $\mathbf{s}$, whereas
this partial linearization might impair convergence at first glance.
Nonetheless, the following provides a convergence guarantee for exADMM, which
is our main theoretical contribution.222The proofs of the theorem and related
lemmas are provided in Appendix A.

###### Theorem 3.1.

Assume that there exist constants $C_{V},C_{U},C_{\mathbf{s}}>0$ such that
$\|\mathbf{V}^{k}\|_{F}^{2}\leq C_{V}$, $\|\mathbf{U}^{k}\|_{F}^{2}\leq
C_{U}$, $\|\mathbf{s}^{k}\|_{2}^{2}\leq C_{s}$ for $\forall k\geq 0$. For
$\rho\geq\max\left(\frac{24\lambda_{ex}^{2}C_{V}C_{\mathbf{s}}}{\underline{\lambda}_{V}},\frac{1}{2}\\!+\\!\sqrt{\frac{1}{4}\\!+\\!6\lambda_{ex}^{2}C_{V}^{2}}\right)$
and
$\gamma\leq\frac{1}{\sqrt{|\mathcal{U}|}((1+\alpha_{0})C_{V}+\bar{\lambda}_{U})+1}$,
where $\bar{\lambda}_{U}\coloneqq\max_{i\in\mathcal{U}}\lambda_{U}^{(i)}$ and
$\underline{\lambda}_{V}\coloneqq\min_{j\in\mathcal{V}}\lambda_{V}^{(j)}$, the
augmented Lagrangian
$L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
converges to some value, and residual norms
$\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F},\|\mathbf{U}^{k+1}-\mathbf{U}^{k}\|_{F},\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}$,
and $\|\mathbf{w}^{k+1}-\mathbf{w}^{k}\|_{2}$ converge to $0$. Furthermore,
the gradients of $L_{\rho}$ with respect to $\mathbf{V}$, $\mathbf{U}$,
$\mathbf{s}$, and $\mathbf{w}$ converge to $0$.

Theorem 3.1 illustrates that the sequence
$\\{\mathbf{U}^{k},\mathbf{V}^{k},\mathbf{s}^{k},\mathbf{w}^{k}\\}$ will
converge to the feasible set, in which
$\mathbf{s}=(1/|\mathcal{U}|)\mathbf{U}^{\top}\mathds{1}$ in Eq. 3 holds.
Moreover, the derivative of the augmented Lagrangian with respect to the
primal variables (i.e., $\mathbf{V}$, $\mathbf{U}$, and $\mathbf{s}$) will
converge to zero, which implies that the limit points of
$\\{\mathbf{U}^{k},\mathbf{V}^{k},\mathbf{s}^{k},\mathbf{w}^{k}\\}$ should be
the saddle points, i.e., the KKT points of Eq. 3 if there exist. Notably, the
above convergence relies on the fact that the objective is strongly convex
with respect to each variable when the other variables are fixed. This
property is inherited from iALS, and therefore exADMM takes advantage of iALS
in both scalability and convergence guarantee while enabling flexible control
of accuracy-fairness tradeoff via a (seemingly) hard-to-optimize regularizer
$R_{ex}(\mathbf{V},\mathbf{U})$.

Table 1. Statistics of the datasets.

Dataset | # of Users | # of Items | # of Interactions
---|---|---|---
ML-20M | $136{,}677$ | $20{,}108$ |  $10$M
MSD | $571{,}355$ | $41{,}140$ |  $33.6$M
Epinions | $6{,}287$ | $3{,}999$ |  $0.13$M

Table 2. Comparison of the scalability (time and space complexity) of the
compared methods

Methods | Time Complexity | Space Complexity | Is it parallelizable across users? | Is it parallelizable across items?
---|---|---|---|---
iALS | $\mathcal{O}\left(\mathrm{nz}(\mathbf{R})d^{2}+(|\mathcal{U}|+|\mathcal{V}|)d^{3}\right)$ | $\mathcal{O}\left((|\mathcal{U}|+|\mathcal{V}|)d\right)$ | ✔ | ✔
Mult-VAE | $\mathcal{O}\left(|\mathcal{U}||\mathcal{V}|d\right)$ | $\mathcal{O}\left(|\mathcal{V}|d\right)$ | ✔ | ✗
FairRec | $\mathcal{O}\left(|\mathcal{U}||\mathcal{V}|(d+\log K)\right)$ | $\mathcal{O}\left(|\mathcal{U}||\mathcal{V}|\right)$ | ✗ | ✗
Multi-FR | $\mathcal{O}\left(\mathrm{nz}(\mathbf{R})|\mathcal{V}|+|\mathcal{U}||\mathcal{V}|d\right)$ | $\mathcal{O}\left(|\mathcal{V}|d\right)$ | ✔ | ✗
exADMM (ours) | $\mathcal{O}\left(\mathrm{nz}(\mathbf{R})d^{2}+|\mathcal{V}|d^{3}\right)$ | $\mathcal{O}\left((|\mathcal{U}|+|\mathcal{V}|)d\right)$ | ✔ | ✔

##  4\. Empirical Evaluation

This section empirically compares exADMM with iALS and existing fair
recommendation methods regarding their effectiveness in accuracy-fairness
control and scalability. Our experiment code is available at
https://github.com/riktor/exADMM-recommender.

###  4.1. Experiment Design

##### Datasets

Our experiments use MovieLens 20M (ML-20M) (Harper and Konstan, 2015), Million
Song Dataset (MSD) (Bertin-Mahieux et al., 2011), and the Epinions dataset
(Epinions) (Massa and Avesani, 2007). Following the standard protocol to
evaluate recommendation effectiveness (Liang et al., 2018; Rendle et al.,
2021), we generate implicit feedback datasets by binarizing the raw explicit
feedback data by keeping interactions with ratings of four or five for ML-20M
and Epinions. For MSD, we use all the recorded interactions as implicit
feedback. Note that, for Epinions, we only retain users and items with more
than 20 interactions following conventional studies (Abdollahpouri et al.,
2017, 2019). Table 1 shows the statistics of the resulting implicit feedback
datasets.

Our experiments follow a strong generalization setting where we use all
interactions of 80% of the users for training and consider the remaining two
sets of 10% of the users as holdout splits. In the validation and testing
phases, each model predicts the preference scores of all items for each user
on the validation and test sets.

##### Compared Methods

Since our aim is to develop a scalable method to enable accuracy-fairness
control, we compare our method against an efficient but unfair method and fair
but inefficient methods. Specifically, we first include the vanilla version of
iALS as a scalable baseline, which does not consider item fairness. Therefore,
against this baseline, our aim is to achieve better accuracy-fairness control
and similar scalability. In addition to iALS, we include Mult-VAE (Liang et
al., 2018), a state-of-the-art deep recommendation method, as a reference to
provide the best achievable recommendation accuracy on each dataset. Besides,
we consider fairness-aware recommendation methods such as the MF-based in-
processing method called _Multi-FR_ (Wu et al., 2022). We also consider a
post-processing method called FairRec (Patro et al., 2020) combined with the
vanilla iALS algorithm to pre-train user-item preference matrix. We thus call
this baseline iALS+FairRec. Against these fairness-aware baselines, our aim is
to achieve a similarly flexible and effective accuracy-fairness control by our
algorithm, which is much more scalable and computationally efficient. Table 2
summarizes the time complexity, space complexity, and parallelizablity of each
method where we can see that Mult-VAE, Multi-FR, and FairRec are particularly
not scalable because their complexity depends on the problematic factor of
$|\mathcal{U}||\mathcal{V}|$.

Throughout the experiments, we train iALS and exADMM for $T=50$ training
epochs, Mult-VAE for $T=200$ epochs, and Multi-FR for $T=500$ epochs with a
constant standard deviation $\sigma=0.1$ for initialization. We tune
$\\{\lambda_{L2},\alpha_{0}\\}$ for iALS and iALS+FairRec,
$\\{\lambda_{L2},\alpha_{0},\lambda_{ex},\rho,\gamma\\}$ for exADMM, and
$\\{p,\tau,\rho_{d}\\}$ for Multi-FR (where $p,\tau,\rho_{d}$ are the user
patience, temperature of smooth rank functions, and dropout rate). To ensure
that $\lambda_{ex}$ and $\rho$ are scale independent of $|\mathcal{U}|$, we
reparametrize $\lambda_{ex}$ and $\rho$ by
$\lambda_{ex}=\lambda_{ex}^{*}\cdot|\mathcal{U}|^{2}$ and
$\rho=\rho^{*}\cdot|\mathcal{U}|^{2}$, respectively, and tune
$\lambda_{ex}^{*}$ and $\rho^{*}$ instead of the original ones. We implement
iALS, iALS+FairRec, and exADMM based on the efficient C++ implementation
provided by Rendle et al. (2021)333https://github.com/google-research/google-
research/tree/master/ials. For a fair comparison, we use frequency-based re-
scaling of $\mathbf{\Lambda}_{U}$ and $\mathbf{\Lambda}_{V}$ (Rendle et al.,
2021) for both iALS and exADMM. Note that we adopt Denoising Auto-Encoder
(Liang et al., 2018) as a backbone model of Multi-FR because it can make
predictions for the holdout users without costly SGD iterations in the testing
phase.

Figure 1. Tradeoff between recommendation accuracy (nDCG@K) and exposure
equality (Gini@K) achieved by each method. Figure 2. Distribution of item
exposure achieved by our method with different hyperparameter ($\lambda_{ex}$)
settings.

##### Evaluation Metrics

We use the normalized cumulative gain (nDCG) as the measure of recommendation
accuracy. To formally define this accuracy metric, let
$\mathcal{V}_{i}\subset\mathcal{V}$ be the held-out items that user $i$
interacts with and $\pi_{i}(k)\in\mathcal{V}$ be the $k$-th item in the ranked
list to be evaluated for $i$. Then, we can define nDCG@$K$ as

(5) $\displaystyle\mathrm{nDCG@}K(i,\pi_{i})$
$\displaystyle\coloneqq\frac{\mathrm{DCG@}K(i,\pi_{i})}{\mathrm{DCG@}K(i,\pi_{i}^{*})},$
where
$\displaystyle\qquad\mathrm{DCG@}K(i,\pi_{i})\coloneqq\sum_{k=1}^{K}\frac{\mathbb{I}\\{\pi_{i}(k)\in\mathcal{V}_{i}\\}}{\log_{2}(k+1)},$

and $\pi_{i}^{*}$ is an ideal ranking for user $i$.

In addition to the accuracy metric (DCG@$K$), we use Gini@$K$ to measure the
exposure inequality based on the Gini index (Atkinson et al., 1970), which is
widely used to evaluate exposure inequality in related research (Do and
Usunier, 2022; Wu et al., 2022). Specifically, Gini@$K$ is defined as follows:

(6)
$\displaystyle\text{Gini@}K(\mathbf{o})\coloneqq\frac{1}{2\left\lVert\mathbf{o}\right\rVert_{1}|\mathcal{V}|^{2}}\sum_{j\in\mathcal{V}}\sum_{l\in\mathcal{V}}|o_{j}-o_{l}|,$

where $\mathbf{o}\in\mathbb{R}^{|\mathcal{V}|}$ is an
$|\mathcal{V}|$-dimensional vector, whose $j$-th element $o_{j}$ indicates the
total exposure given to item $j$, i.e.,
$o_{j}\coloneqq\sum_{i}o(i,j;\pi_{i})$. In our experiments, we define
$o(i,j;\pi_{i})\coloneqq\mathbb{I}\\{\pi_{i}^{-1}(j)\leq
K\\}/\log_{2}(\pi_{i}^{-1}(j)+1)$ following the examination model of the DCG
metric. Note that a lower value of Gini@K indicates that recommendations are
more fair towards the items, but to do so without sacrificing accuracy and
scalability is particularly challenging.

###  4.2. Results and Discussion

##### Accuracy-Fairness Tradeoff

First, we evaluate and compare how well each method can control the tradeoff
between recommendation accuracy (nDCG@K) and exposure fairness (Gini@K) in
Figure 1 on the three datasets and three different values of $K$. In the
figures, we present the Pareto frontier for exADMM, Multi-FR, and iALS+FairRec
across various hyperparameter settings obtained through grid
search.444Hyperparameters range from $\lambda_{L2}\in[1\mathrm{e}{-4},1.0]$,
$\alpha_{0}\in[1\mathrm{e}{-4},1.0]$,
$\lambda_{ex}^{*}\in[1\mathrm{e}{-10},1\mathrm{e}{-4}]$,
$\rho^{*}\in[1\mathrm{e}{-10},1\mathrm{e}{-4}]$,
$\gamma=\\{0.001,0.01,0.05\\}$, $\beta\in[0.1,1.0]$, and
$p\in[0.6,1.0],\tau\in[1\mathrm{e}{-3},1.0],\rho_{d}=[0.1,0.9]$. For
iALS+FairRec, we first select the best setting of $\alpha_{0}$ and
$\lambda_{L2}$ in iALS in terms of nDCG@K and then tune the hyperparameters of
FairRec, namely, the length of rankings $K$ and the scale $l\in(0,1]$ of the
minimum allocation constraint for each item
$l\cdot(K\cdot|\mathcal{U}|)/|\mathcal{V}|$. We search each parameter in a
logarithmic scale unless otherwise noted. For all methods, we use $d=32$ for
Epinions, $d=256$ for ML-20M, and $d=512$ for MSD. Note that neither iALS nor
Mult-VAE incorporates a fairness regularizer or a specific hyperparameter for
managing the accuracy-fairness tradeoff in their objectives. Therefore, we
adjust their hyperparameters based on nDCG@K in the validation split and then
report nDCG@K and Gini@K in the test split, represented as (blue and yellow)
dots. This is why these methods achieve better accuracy than fairness-aware
methods including exADMM, but it is also true that they produce substantial
unfairness among items according to their Gini@K. It should also be noted that
Multi-FR is infeasible on ML-20M and MSD, and iALS+FairRec is infeasible on
MSD. This is due to their excessive time and space complexity, which depend on
$|\mathcal{U}||\mathcal{V}|$. As a result, their results are not shown in the
figures corresponding to these datasets.

Table 3. Computation time (seconds) to complete training and inference of the
compared methods. NA means that the method was infeasible on the corresponding
dataset.

| Epinions | ML-20M | MSD
---|---|---|---
| $K=5$ | $K=20$ | $K=10$ | $K=50$ | $K=10$ | $K=50$
Methods | Train | Inference | Train | Inference | Train | Inference | Train | Inference | Train | Inference | Train | Inference
iALS | 0.7 | 0.1 | 0.7 | 0.1 | 67.2 | 4.2 | 67.2 | 4.2 | $857.6$ | 39.4 | $857.6$ | 39.4
Mult-VAE | 72.6 | 0.5 | 72.6 | 0.5 | $1{,}872.6$ | 14.7 | $1{,}872.6$ | 14.7 | $22{,}580.1$ | 145.2 | $22{,}580.1$ | 145.2
FairRec | 0.7 | 0.3 | 0.7 | 0.3 | 67.2 | 23.3 | 67.2 | 30.1 | $857.6$ | NA | $857.6$ | NA
Multi-FR | $4{,}011.3$ | 0.5 | $4{,}011.3$ | 0.5 | NA | NA | NA | NA | NA | NA | NA | NA
exADMM (ours) | 1.7 | 0.4 | 1.7 | 0.4 | 98.0 | 9.5 | 98.0 | 9.5 | $1{,}533.1$ | 89.6 | $1{,}533.1$ | 89.6

From the figures, we can first see that exADMM achieves a similar accuracy-
fairness tradeoff compared to Multi-FR on Epinions even though exADMM is much
more scalable. Next, the comparisons between exADMM and iALS+FairRec (which
has much larger time and space complexity compared to exADMM) are complicated,
and it is not straightforward to determine which method performs better in
terms of accuracy-fairness tradeoff. However, an interesting trend emerges
from our empirical results on Epinions and ML-20M (FairRec is infeasible on
MSD). That is, exADMM is likely to achieve a better accuracy-fairness tradeoff
than iALS+FairRec when we need to achieve high recommendation accuracy, while
iALS+FairRec is likely to perform better in terms of the tradeoff when
enforcing a strong fairness requirement. This interesting difference can be
attributed to the fact that exADMM is an in-processing method while FairRec is
a post-processing method. Overall, it is remarkable that exADMM, which is much
more scalable, achieves a competitive effectiveness in terms of accuracy-
fairness control compared to Multi-FR and iALS+FairRec, which do not consider
scalability.

To visualize how flexibly exADMM can control exposure distribution via its
hyperparameter ($\lambda_{ex}$), Figure 2 illustrates the cumulative item
exposure (called the Lorenz curve) induced by exADMM with several different
values of $\lambda_{ex}$ on Epinions (top row), ML-20M (mid row), and MSD
(bottom row). Each curve in the figures shows the cumulative item exposure
relative to the case when $\lambda_{ex}=0.0$ (i.e., the vanilla version of
iALS), and this is why $\lambda_{ex}=0.0$ (iALS) always has flat lines. In
addition to the curve for $\lambda_{ex}=0.0$, we present curves induced by
three other values of $\lambda_{ex}$ (small, medium, and large) for each
dataset. Specifically, we use $\lambda_{ex}=2\mathrm{e}{-4}$ (small),
$\lambda_{ex}=3\mathrm{e}{-3}$ (medium), and $\lambda_{ex}=9\mathrm{e}{-2}$
(large) for Epinions, $\lambda_{ex}=1\mathrm{e}{-10}$ (small),
$\lambda_{ex}=1\mathrm{e}{-7}$ (medium), and $\lambda_{ex}=3\mathrm{e}{-6}$
(large) for ML-20M, and $\lambda_{ex}=1\mathrm{e}{-12}$ (small),
$\lambda_{ex}=1\mathrm{e}{-10}$ (medium), and $\lambda_{ex}=1\mathrm{e}{-3}$
(large) for MSD. Figure 2 demonstrates that exADMM clearly achieves fairer
exposure distribution compared to iALS ($\lambda_{ex}=0.0$). We can also see
that we can flexibly and accurately control the item exposure distribution via
the hyperparameter $\lambda_{ex}$ of our method. That is, we can see a
monotonic relationship between $\lambda_{ex}$ and fairness of exposure
distribution, suggesting that $\lambda_{ex}$ of exADMM works as an appropriate
controller of item fairness.

##### Computational Complexity

Finally, we empirically evaluate the computational efficiency of the methods
in terms of both training and inference. Table 3 reports the average elapsed
time to complete training and inference of each method on the three datasets
and two values of $K$. Note that the “NA” values that we see for Multi-FR and
FairRec indicate that their training and/or inference are infeasible.

From the table, it is evident that exADMM achieves computational efficiency
close to iALS across all datasets and all values of $K$, which implies that
our method is able to control the accuracy-fairness tradeoff while retaining
scalability. We also observe that Mult-VAE needs approximately 15-40 times
longer training time compared to exADMM. This issue of Mult-VAE will be
exacerbated as the item space grows, given its time complexity of
$\mathcal{O}\left(|\mathcal{U}||\mathcal{V}|)d\right)$. The training procedure
of Multi-FR is about 2,350 times slower than that of exADMM on Epinions.
Besides, Multi-FR is infeasible on ML-20M and MSD. FairRec suffers from longer
inference time compared to iALS and exADMM, and on MSD, its inference becomes
infeasible. Therefore, even though Multi-FR and FairRec show their usefulness
in terms of exposure control on datasets with a limited size, due to their
inefficiency and large complexity (as in Table 2), they are impractical for
most industry-scale systems, which can even be larger than MSD. This
demonstrates that exADMM is the first method that enables an effective control
of accuracy-fairness tradeoff and is scalable to systems of practical size.

##  5\. Conclusion

The feasibility of exposure-controllable recommendation is indispensable for
solving accuracy-fairness tradeoff in practice, however, it has been
disregarded in academic research. Therefore, this work develops a novel
scalable method called exADMM, to enable flexible exposure control. Despite
the technical difficulty in handling the exposure regularizer in parallel,
exADMM achieves this while maintaining scalability with yet provable
convergence guarantees. Empirical evaluations are promising and demonstrate
that our method is the first to achieve flexible control of accuracy-fairness
tradeoff in a scalable and computationally efficient way. Our work also raises
several intriguing questions for future studies such as extensions to more
refined fairness regularizers beyond the mere second moment, a scalable post-
processing approach to control fairness, and a scalable control of impact-
based fairness (Saito and Joachims, 2022).

## References

  * (1)
  * Abdollahpouri et al. (2020) Himan Abdollahpouri, Gediminas Adomavicius, Robin Burke, Ido Guy, Dietmar Jannach, Toshihiro Kamishima, Jan Krasnodebski, and Luiz Pizzato. 2020.  Multistakeholder recommendation: Survey and research directions. In _UMUAI_.
  * Abdollahpouri et al. (2017) Himan Abdollahpouri, Robin Burke, and Bamshad Mobasher. 2017.  Controlling popularity bias in learning-to-rank recommendation. In _RecSys_.
  * Abdollahpouri et al. (2019) Himan Abdollahpouri, Robin Burke, and Bamshad Mobasher. 2019.  Managing popularity bias in recommender systems with personalized re-ranking. In _International FLAIRS Conference_.
  * Atkinson et al. (1970) Anthony B Atkinson et al. 1970\.  On the measurement of inequality. In _Journal of Economic Theory_.
  * Bayer et al. (2017) Immanuel Bayer, Xiangnan He, Bhargav Kanagal, and Steffen Rendle. 2017.  A generic coordinate descent framework for learning from implicit feedback. In _WWW_.
  * Bertin-Mahieux et al. (2011) Thierry Bertin-Mahieux, Daniel PW Ellis, Brian Whitman, and Paul Lamere. 2011.  The million song dataset. In _ISMIR_.
  * Biega et al. (2018) Asia J Biega, Krishna P Gummadi, and Gerhard Weikum. 2018.  Equity of attention: Amortizing individual fairness in rankings. In _SIGIR_.
  * Boyd et al. (2011) Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al. 2011\.  Distributed optimization and statistical learning via the alternating direction method of multipliers. In _Foundations and Trends® in Machine learning_.
  * Brynjolfsson et al. (2011) Erik Brynjolfsson, Yu Hu, and Duncan Simester. 2011.  Goodbye pareto principle, hello long tail: The effect of search costs on the concentration of product sales. In _Management Science_.
  * Burke (2017) Robin Burke. 2017.  Multisided fairness for recommendation. In _arXiv preprint arXiv:1707.00093_.
  * Burke et al. (2018) Robin Burke, Nasim Sonboli, and Aldo Ordonez-Gauger. 2018.  Balanced neighborhoods for multi-sided fairness in recommendation. In _FAccT_.
  * Chen et al. (2016) Caihua Chen, Bingsheng He, Yinyu Ye, and Xiaoming Yuan. 2016.  The direct extension of ADMM for multi-block convex minimization problems is not necessarily convergent. In _Mathematical Programming_.
  * Chen et al. (2023) Chong Chen, Weizhi Ma, Min Zhang, Chenyang Wang, Yiqun Liu, and Shaoping Ma. 2023.  Revisiting negative sampling vs. non-sampling in implicit recommendation. In _TOIS_.
  * Cheng et al. (2014) Yao Cheng, Liang Yin, and Yong Yu. 2014.  Lorslim: Low rank sparse linear methods for top-n recommendations. In _ICDM_.
  * Deng et al. (2017) Wei Deng, Ming-Jun Lai, Zhimin Peng, and Wotao Yin. 2017.  Parallel multi-block ADMM with O (1/k) convergence. In _Journal of Scientific Computing_.
  * Do et al. (2021) Virginie Do, Sam Corbett-Davies, Jamal Atif, and Nicolas Usunier. 2021.  Two-sided fairness in rankings via Lorenz dominance. In _NIPS_.
  * Do and Usunier (2022) Virginie Do and Nicolas Usunier. 2022.  Optimizing generalized Gini indices for fairness in rankings. In _SIGIR_.
  * Duchi and Singer (2009) John Duchi and Yoram Singer. 2009.  Efficient online and batch learning using forward backward splitting. In _JMLR_.
  * Fleder and Hosanagar (2009) Daniel Fleder and Kartik Hosanagar. 2009.  Blockbuster culture’s next rise or fall: The impact of recommender systems on sales diversity. In _Management Science_.
  * Frank and Wolfe (1956) Marguerite Frank and Philip Wolfe. 1956.  An algorithm for quadratic programming. In _Naval Research Logistics Quarterly_.
  * Harper and Konstan (2015) F Maxwell Harper and Joseph A Konstan. 2015.  The movielens datasets: History and context. In _TiiS_.
  * He et al. (2017) Xiangnan He, Lizi Liao, Hanwang Zhang, Liqiang Nie, Xia Hu, and Tat-Seng Chua. 2017.  Neural collaborative filtering. In _WWW_.
  * Helberger et al. (2018) Natali Helberger, Kari Karppinen, and Lucia D’acunto. 2018.  Exposure diversity as a design principle for recommender systems. In _Information, Communication & Society_.
  * Hu et al. (2008) Yifan Hu, Yehuda Koren, and Chris Volinsky. 2008.  Collaborative filtering for implicit feedback datasets. In _ICDM_.
  * Ioannidis et al. (2019) Vassilis N Ioannidis, Ahmed S Zamzam, Georgios B Giannakis, and Nicholas D Sidiropoulos. 2019.  Coupled graphs and tensor factorization for recommender systems and community detection. In _TKDE_.
  * Jaggi (2013) Martin Jaggi. 2013.  Revisiting Frank-Wolfe: Projection-free sparse convex optimization. In _ICML_.
  * Kamishima et al. (2013) Toshihiro Kamishima, Shotaro Akaho, Hideki Asoh, and Jun Sakuma. 2013.  Efficiency Improvement of Neutrality-Enhanced Recommendation. In _Human Decision Making in Recommender Systems (Decisions@ RecSys’ 13)_.
  * Kamishima et al. (2011) Toshihiro Kamishima, Shotaro Akaho, and Jun Sakuma. 2011.  Fairness-aware learning through regularization approach. In _ICDM_.
  * Liang et al. (2018) Dawen Liang, Rahul G Krishnan, Matthew D Hoffman, and Tony Jebara. 2018.  Variational autoencoders for collaborative filtering. In _WWW_.
  * Liu et al. (2019) Qinghua Liu, Xinyue Shen, and Yuantao Gu. 2019.  Linearized ADMM for nonconvex nonsmooth optimization with convergence analysis. In _IEEE Access_.
  * Massa and Avesani (2007) Paolo Massa and Paolo Avesani. 2007.  Trust-aware recommender systems. In _RecSys_.
  * Memarrast et al. (2021) Omid Memarrast, Ashkan Rezaei, Rizal Fathony, and Brian Ziebart. 2021.  Fairness for Robust Learning to Rank. In _arXiv preprint arXiv:2112.06288_.
  * Morik et al. (2020) Marco Morik, Ashudeep Singh, Jessica Hong, and Thorsten Joachims. 2020.  Controlling fairness and bias in dynamic learning-to-rank. In _SIGIR_.
  * Oosterhuis (2021) Harrie Oosterhuis. 2021.  Computationally efficient optimization of plackett-luce ranking models for relevance and fairness. In _SIGIR_.
  * Park and Tuzhilin (2008) Yoon-Joo Park and Alexander Tuzhilin. 2008.  The long tail of recommender systems and how to leverage it. In _RecSys_.
  * Patro et al. (2020) Gourab K Patro, Arpita Biswas, Niloy Ganguly, Krishna P Gummadi, and Abhijnan Chakraborty. 2020.  Fairrec: Two-sided fairness for personalized recommendations in two-sided platforms. In _WWW_.
  * Plackett (1975) Robin L Plackett. 1975.  The analysis of permutations. In _Journal of the Royal Statistical Society: Series C (Applied Statistics)_.
  * Rendle et al. (2020) Steffen Rendle, Walid Krichene, Li Zhang, and John Anderson. 2020.  Neural collaborative filtering vs. matrix factorization revisited. In _RecSys_.
  * Rendle et al. (2021) Steffen Rendle, Walid Krichene, Li Zhang, and Yehuda Koren. 2021.  Revisiting the Performance of iALS on Item Recommendation Benchmarks. In _arXiv preprint arXiv:2110.14037_.
  * Resnick et al. (2013) Paul Resnick, R Kelly Garrett, Travis Kriplean, Sean A Munson, and Natalie Jomini Stroud. 2013.  Bursting your (filter) bubble: strategies for promoting diverse exposure. In _CSCW_.
  * Rockafellar (1976) R Tyrrell Rockafellar. 1976.  Monotone operators and the proximal point algorithm. In _SIAM Journal on Control and Optimization_.
  * Saito (2020) Yuta Saito. 2020.  Unbiased pairwise learning from biased implicit feedback. In _Proceedings of the 2020 ACM SIGIR on International Conference on Theory of Information Retrieval_. 5–12.
  * Saito and Joachims (2022) Yuta Saito and Thorsten Joachims. 2022.  Fair Ranking as Fair Division: Impact-Based Individual Fairness in Ranking. In _KDD_.
  * Saito et al. (2020) Yuta Saito, Suguru Yaginuma, Yuta Nishino, Hayato Sakata, and Kazuhide Nakata. 2020.  Unbiased recommender learning from missing-not-at-random implicit feedback. In _Proceedings of the 13th International Conference on Web Search and Data Mining_. 501–509.
  * Santos et al. (2010) Rodrygo LT Santos, Craig Macdonald, and Iadh Ounis. 2010.  Exploiting query reformulations for web search result diversification. In _WWW_.
  * Sherman and Morrison (1950) Jack Sherman and Winifred J Morrison. 1950.  Adjustment of an inverse matrix corresponding to a change in one element of a given matrix. In _The Annals of Mathematical Statistics_.
  * Singh and Joachims (2018) Ashudeep Singh and Thorsten Joachims. 2018.  Fairness of exposure in rankings. In _KDD_.
  * Singh and Joachims (2019) Ashudeep Singh and Thorsten Joachims. 2019.  Policy learning for fairness in ranking. In _NIPS_.
  * Smith et al. (2017) Shaden Smith, Alec Beri, and George Karypis. 2017.  Constrained tensor factorization with accelerated AO-ADMM. In _ICPP_.
  * Steck (2018) Harald Steck. 2018.  Calibrated recommendations. In _RecSys_.
  * Steck et al. (2020) Harald Steck, Maria Dimakopoulou, Nickolai Riabov, and Tony Jebara. 2020.  Admm slim: Sparse recommendations for many users. In _WSDM_.
  * Steck and Liang (2021) Harald Steck and Dawen Liang. 2021.  Negative Interactions for Improved Collaborative Filtering: Don’t go Deeper, go Higher. In _RecSys_.
  * Usunier et al. (2022) Nicolas Usunier, Virginie Do, and Elvis Dohmatob. 2022.  Fast online ranking with fairness of exposure. In _FAccT_.
  * Wang et al. (2014) Huahua Wang, Arindam Banerjee, and Zhi-Quan Luo. 2014.  Parallel direction method of multipliers. In _NIPS_.
  * Woodbury (1950) Max A Woodbury. 1950.  _Inverting modified matrices_.  Statistical Research Group.
  * Wu et al. (2022) Haolun Wu, Chen Ma, Bhaskar Mitra, Fernando Diaz, and Xue Liu. 2022.  A Multi-objective Optimization Framework for Multi-stakeholder Fairness-aware Recommendation. In _TOIS_.
  * Wu et al. (2021) Yao Wu, Jian Cao, Guandong Xu, and Yudong Tan. 2021.  TFROM: A two-sided fairness-aware recommendation model for both customers and providers. In _SIGIR_.
  * Yadav et al. (2021) Himank Yadav, Zhengxiao Du, and Thorsten Joachims. 2021.  Policy-gradient training of fair and unbiased ranking functions. In _SIGIR_.
  * Yao and Huang (2017) Sirui Yao and Bert Huang. 2017.  Beyond parity: Fairness objectives for collaborative filtering. In _NIPS_.
  * Yu et al. (2014a) Hsiang-Fu Yu, Cho-Jui Hsieh, Si Si, and Inderjit S Dhillon. 2014a.  Parallel matrix factorization for recommender systems. In _Knowledge and Information Systems_.
  * Yu et al. (2014b) Zhi-Qin Yu, Xing-Jian Shi, Ling Yan, and Wu-Jun Li. 2014b.  Distributed stochastic ADMM for matrix factorization. In _CIKM_.
  * Zehlike and Castillo (2020) Meike Zehlike and Carlos Castillo. 2020.  Reducing disparate exposure in ranking: A learning to rank approach. In _WWW_.
  * Zhou et al. (2008) Yunhong Zhou, Dennis Wilkinson, Robert Schreiber, and Rong Pan. 2008.  Large-scale parallel collaborative filtering for the netflix prize. In _AAIM_.

##  Appendix A Proofs of Convergence Guarantee

###  A.1. Proof of 3.1

###### Proof of Theorem 3.1.

In the proof, we use the following lemma on the smoothness of $g$:

###### Lemma A.1.

For any $\mathbf{V},\mathbf{V}^{\prime}\in\mathbb{R}^{|\mathcal{V}|\times d}$,
$\mathbf{s},\mathbf{s}^{\prime}\in\mathbb{R}^{d}$, and
$\mathbf{U},\mathbf{U}^{\prime}\in\mathbb{R}^{|\mathcal{U}|\times d}$,
function $g$ satisfies the following inequalities:

$\displaystyle\|\nabla_{\mathbf{V}}g(\mathbf{V},\mathbf{U},\mathbf{s})-\nabla_{\mathbf{V}}g(\mathbf{V}^{\prime},\mathbf{U},\mathbf{s})\|_{F}\leq\sqrt{|\mathcal{V}|}\left(\left(1+\alpha_{0}\right)\|\mathbf{U}\|_{F}^{2}+\lambda_{ex}\|\mathbf{s}\|_{2}^{2}+\bar{\lambda}_{V}\right)\|\mathbf{V}-\mathbf{V}^{\prime}\|_{F},$
$\displaystyle\|\nabla_{\mathbf{U}}g(\mathbf{V},\mathbf{U},\mathbf{s})-\nabla_{\mathbf{U}}g(\mathbf{V},\mathbf{U}^{\prime},\mathbf{s}^{\prime})\|_{F}\leq\sqrt{|\mathcal{U}|}\left((1+\alpha_{0})\|\mathbf{V}\|_{F}^{2}+\bar{\lambda}_{U}\right)\|\mathbf{U}-\mathbf{U}^{\prime}\|_{F},$
$\displaystyle\|\nabla_{\mathbf{s}}g(\mathbf{V},\mathbf{U},\mathbf{s})-\nabla_{\mathbf{s}}g(\mathbf{V},\mathbf{U}^{\prime},\mathbf{s}^{\prime})\|_{2}\leq\lambda_{ex}\|\mathbf{V}\|_{F}^{2}\|\mathbf{s}-\mathbf{s}^{\prime}\|_{2},$
$\displaystyle\|\nabla_{\mathbf{s}}g(\mathbf{V},\mathbf{U},\mathbf{s})-\nabla_{\mathbf{s}}g(\mathbf{V}^{\prime},\mathbf{U},\mathbf{s})\|_{2}\leq\lambda_{ex}(\|\mathbf{V}\|_{F}+\|\mathbf{V}^{\prime}\|_{F})\|\mathbf{s}\|_{2}\|\mathbf{V}-\mathbf{V}^{\prime}\|_{F},$

where $\bar{\lambda}_{U}\coloneqq\max_{i\in\mathcal{U}}\lambda_{U}^{(i)}$ and
$\bar{\lambda}_{V}\coloneqq\max_{j\in\mathcal{V}}\lambda_{V}^{(j)}$.

We prove the first part of the theorem. We decompose the difference of
$L_{\rho}$ before and after a single epoch update into that before and after
each alternating step.

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k+1})-L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
$\displaystyle=\left(L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k},\mathbf{w}^{k})-L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})\right)$
$\displaystyle\phantom{=}+\left(L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k})-L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k},\mathbf{w}^{k})\right)$
(7)
$\displaystyle\phantom{=}+\left(L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k+1})-L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k})\right).$

Lemma A.1 implies the upper bound on each term in the RHS:

###### Lemma A.2.

The update of $\mathbf{V}$ and $\mathbf{U}$ in the $(k+1)$-step satisfies

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k},\mathbf{w}^{k})-L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
$\displaystyle\leq\frac{\sqrt{|\mathcal{U}|}((1+\alpha_{0})C_{V}+\bar{\lambda}_{U})-1/\gamma}{2}\|\mathbf{U}^{k+1}-\mathbf{U}^{k}\|_{F}^{2}-\frac{\underline{\lambda}_{V}}{2}\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2}.$

###### Lemma A.3.

The update of $\mathbf{s}$ in the $(k+1)$-th step satisfies

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k})-L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k},\mathbf{w}^{k})\leq-\frac{\rho}{2}\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}^{2}.$

###### Lemma A.4.

The update of $\mathbf{w}$ in the $(k+1)$-th step satisfies

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k+1})-L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k})$
$\displaystyle\leq\frac{3\lambda_{ex}^{2}C_{V}^{2}}{\rho}\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}^{2}+\frac{6\lambda_{ex}^{2}C_{V}C_{s}}{\rho}\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2}.$

By using Section A.1 and Lemmas A.2 to A.4, under the assumptions
$\rho\geq\max\left(\frac{24\lambda_{ex}^{2}C_{V}C_{\mathbf{s}}}{\underline{\lambda}_{V}},\frac{1}{2}+\sqrt{\frac{1}{4}+6\lambda_{ex}^{2}C_{V}^{2}}\right)$
and
$\gamma\leq\frac{1}{\sqrt{|\mathcal{U}|}((1+\alpha_{0})C_{V}+\bar{\lambda}_{U})+1}$,
we have

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k+1})-L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
$\displaystyle\leq\frac{\sqrt{|\mathcal{U}|}((1+\alpha_{0})C_{V}+\bar{\lambda}_{U})-1/\gamma}{2}\|\mathbf{U}^{k+1}-\mathbf{U}^{k}\|_{F}^{2}-\frac{\underline{\lambda}_{V}}{2}\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2}-\frac{\rho}{2}\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}^{2}$
$\displaystyle\phantom{\leq}+\frac{3\lambda_{ex}^{2}C_{V}^{2}}{\rho}\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}^{2}+\frac{6\lambda_{ex}^{2}C_{V}C_{s}}{\rho}\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2}$
$\displaystyle=\frac{\sqrt{|\mathcal{U}|}((1+\alpha_{0})C_{V}+\bar{\lambda}_{U})-1/\gamma}{2}\|\mathbf{U}^{k+1}-\mathbf{U}^{k}\|_{F}^{2}+\left(-\frac{\underline{\lambda}_{V}}{2}+\frac{6\lambda_{ex}^{2}C_{V}C_{s}}{\rho}\right)\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2}+\left(-\frac{\rho}{2}+\frac{3\lambda_{ex}^{2}C_{V}^{2}}{\rho}\right)\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}^{2}$
(8)
$\displaystyle\leq-\frac{1}{2}\|\mathbf{U}^{k+1}-\mathbf{U}^{k}\|_{F}^{2}-\frac{\underline{\lambda}_{V}}{4}\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2}-\frac{1}{2}\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}^{2}\leq
0.$

Therefore,
$L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$ is
monotonically decreasing.

Here, we obtain the following lower bound on
$L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$:

###### Lemma A.5.

$\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k}$, and $\mathbf{w}^{k}$ updated
by exADMM satisfy

$\displaystyle
L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})\geq\frac{\rho-\lambda_{ex}C_{V}}{2}\left\|\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}-\mathbf{s}^{k}\right\|_{2}^{2}.$

Thus, when $\rho\geq\lambda_{ex}C_{V}$ holds,
$L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$ is
lower bounded by $0$. Therefore, owing to its monotonic decrease,
$L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
converges to some constant value, and
$L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k+1})-L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
converges to $0$. From Section A.1 and the fact that
$L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k+1})-L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
converges to $0$, $\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}$,
$\|\mathbf{U}^{k+1}-\mathbf{U}^{k}\|_{F}$, and
$\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}$ also converge to $0$. Finally, from
Lemma A.4, we have:

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k+1})-L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k})$
$\displaystyle=\rho\|\mathbf{w}^{k+1}-\mathbf{w}^{k}\|_{2}^{2}$
$\displaystyle\leq\frac{3\lambda_{ex}^{2}C_{V}^{2}}{\rho}\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}^{2}+\frac{6\lambda_{ex}^{2}C_{V}C_{s}}{\rho}\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2},$

and then we can also state that $\|\mathbf{w}^{k+1}-\mathbf{w}^{k}\|_{2}$
converges to $0$.

We next prove the second part of the theorem. Since $\mathbf{V}^{k+1}$
minimizes $L_{\rho}(\mathbf{V},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$,
it holds that
$\nabla_{\mathbf{V}}L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})=0$,
and we obtain the following inequality from Lemma A.1:

$\displaystyle\|\nabla_{\mathbf{V}}L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})\|_{F}$
$\displaystyle=\|\nabla_{\mathbf{V}}L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})-\nabla_{\mathbf{V}}L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})\|_{F}$
$\displaystyle=\|\nabla_{\mathbf{V}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})-\nabla_{\mathbf{V}}g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})\|_{F}$
$\displaystyle\leq\sqrt{|\mathcal{V}|}\left(\left(1+\alpha_{0}\right)C_{U}+\lambda_{ex}C_{s}+\bar{\lambda}_{V}\right)\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}.$

Since $\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}$ converges to $0$,
$\nabla_{\mathbf{V}}L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
converge to $0$. Similarly, we have:

$\displaystyle\|\nabla_{\mathbf{U}}L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})\|_{F}$
$\displaystyle=\left\|\nabla_{\mathbf{U}}g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})+\frac{\rho}{|\mathcal{U}|^{2}}\mathds{1}\mathds{1}^{\top}\mathbf{U}^{k}+\frac{\rho}{|\mathcal{U}|}\mathds{1}(\mathbf{w}^{k}-\mathbf{s}^{k})^{\top}\right\|_{F}$
$\displaystyle=\left\|\nabla_{\mathbf{U}}g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})-\nabla_{\mathbf{U}}g(\mathbf{V}^{k},\mathbf{U}^{k-1},\mathbf{s}^{k-1})+\frac{\rho}{|\mathcal{U}|}\mathds{1}(\mathbf{w}^{k}-\mathbf{w}^{k-1})^{\top}-\frac{\rho}{|\mathcal{U}|}\mathds{1}(\mathbf{s}^{k}-\mathbf{s}^{k-1})^{\top}-\frac{1}{\gamma}(\mathbf{U}^{k}-\mathbf{U}^{k-1})\right\|_{F}$
$\displaystyle\leq\|\nabla_{\mathbf{U}}g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})-\nabla_{\mathbf{U}}g(\mathbf{V}^{k},\mathbf{U}^{k-1},\mathbf{s}^{k-1})\|_{F}+\frac{1}{\gamma}\|\mathbf{U}^{k}-\mathbf{U}^{k-1}\|_{F}+\frac{\rho}{|\mathcal{U}|}\|\mathds{1}(\mathbf{w}^{k}-\mathbf{w}^{k-1})^{\top}\|_{F}+\frac{\rho}{|\mathcal{U}|}\|\mathds{1}(\mathbf{s}^{k}-\mathbf{s}^{k-1})^{\top}\|_{F}$
$\displaystyle\leq\sqrt{|\mathcal{U}|}\left((1+\alpha_{0})C_{V}+\bar{\lambda}_{U}\right)\|\mathbf{U}^{k}-\mathbf{U}^{k-1}\|_{F}+\frac{1}{\gamma}\|\mathbf{U}^{k}-\mathbf{U}^{k-1}\|_{F}+\frac{\rho}{|\mathcal{U}|}\|\mathds{1}(\mathbf{w}^{k}-\mathbf{w}^{k-1})^{\top}\|_{F}+\frac{\rho}{|\mathcal{U}|}\|\mathds{1}(\mathbf{s}^{k}-\mathbf{s}^{k-1})^{\top}\|_{F},$

where the second equality follows from the fact that $\mathbf{U}^{k}$
minimizes
$\frac{\rho}{2}\|\frac{1}{|\mathcal{U}|}\mathbf{U}^{\top}\mathds{1}+\mathbf{w}^{k-1}-\mathbf{s}^{k-1}\|_{2}^{2}-\frac{\rho}{2}\|\mathbf{w}^{k-1}\|_{2}^{2}+\frac{1}{2\gamma}\|\mathbf{U}-\mathbf{U}^{k-1}\|_{F}^{2}+\langle\mathbf{U}-\mathbf{U}^{k-1},\nabla_{\mathbf{U}}g(\mathbf{V}^{k},\mathbf{U}^{k-1},\mathbf{s}^{k-1})\rangle_{F}$.
Since $\|\mathbf{U}^{k}-\mathbf{U}^{k-1}\|_{F}$,
$\|\mathbf{w}^{k}-\mathbf{w}^{k-1}\|_{2}$ converge to $0$, this inequality
implies that
$\nabla_{\mathbf{U}}L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
converges to $0$.

Because $\mathbf{s}^{k}$ minimizes
$L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s},\mathbf{w}^{k-1})$, we also
have

$\displaystyle\|\nabla_{\mathbf{s}}L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})\|_{2}$
$\displaystyle=\|\nabla_{\mathbf{s}}L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k-1})-\nabla_{\mathbf{s}}L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})\|_{2}$
$\displaystyle=\rho\|\mathbf{w}^{k}-\mathbf{w}^{k-1}\|_{2}.$

Thus, since $\|\mathbf{w}^{k}-\mathbf{w}^{k-1}\|_{2}$ converges to $0$,
$\nabla_{\mathbf{s}}L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
converges to $0$.

Finally, it holds that

$\displaystyle\|\nabla_{\mathbf{w}}L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})\|_{2}$
$\displaystyle=\rho\left\|\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}-\mathbf{s}^{k}\right\|_{2}$
$\displaystyle=\rho\|\mathbf{w}^{k}-\mathbf{w}^{k-1}\|_{2},$

and hence
$\nabla_{\mathbf{w}}L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
converges to $0$. ∎

###  A.2. Proof of Lemma A.2

###### Proof.

From the definition of
$L_{\rho}(\mathbf{V},\mathbf{U},\mathbf{s}^{k},\mathbf{w}^{k})$, we have:

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k},\mathbf{w}^{k})-L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
$\displaystyle=g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k})+\frac{\rho}{2}\left\lVert\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}^{k}\right\rVert_{2}^{2}-\frac{\rho}{2}\|\mathbf{w}^{k}\|_{2}^{2}-g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})-\frac{\rho}{2}\left\lVert\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}-\mathbf{s}^{k}\right\rVert_{2}^{2}+\frac{\rho}{2}\|\mathbf{w}^{k}\|_{2}^{2}$
(9)
$\displaystyle=g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k})-g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})+g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})-g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})+\frac{\rho}{2}\left\lVert\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}^{k}\right\rVert_{2}^{2}-\frac{\rho}{2}\left\lVert\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}-\mathbf{s}^{k}\right\rVert_{2}^{2}.$

Denoting the Gramian of $\mathbf{U}$ by
$\mathbf{G}_{U}\coloneqq\mathbf{U}^{\top}\mathbf{U}$, we have

$\displaystyle\langle\nabla_{\mathbf{V}}g(\mathbf{V},\mathbf{U},\mathbf{s})-\nabla_{\mathbf{V}}g(\mathbf{V}^{\prime},\mathbf{U},\mathbf{s}),\mathbf{V}-\mathbf{V}^{\prime}\rangle_{F}$
$\displaystyle=\sum_{j\in\mathcal{V}}\Biggl{\langle}\left(\sum_{i\in\mathcal{U}}r_{i,j}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}+\alpha_{0}\mathbf{G}_{U}+\lambda_{ex}\mathbf{s}\mathbf{s}^{\top}+\lambda_{V}^{(j)}\mathbf{I}\right)\mathbf{v}_{j}-\sum_{i\in\mathcal{U}}r_{i,j}\mathbf{u}_{i}-\left(\sum_{i\in\mathcal{U}}r_{i,j}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}+\alpha_{0}\mathbf{G}_{U}+\lambda_{ex}\mathbf{s}\mathbf{s}^{\top}+\lambda_{V}^{(j)}\mathbf{I}\right)\mathbf{v}_{j}^{\prime}+\sum_{i\in\mathcal{U}}r_{i,j}\mathbf{u}_{i},\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\Biggr{\rangle}$
$\displaystyle=\sum_{j\in\mathcal{V}}\left\langle\left(\sum_{i\in\mathcal{U}}r_{i,j}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}\right)(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}),\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\right\rangle+\sum_{j\in\mathcal{V}}\left\langle\alpha_{0}\mathbf{G}_{U}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}),\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\right\rangle$
$\displaystyle+\sum_{j\in\mathcal{V}}\left\langle\lambda_{ex}\mathbf{s}\mathbf{s}^{\top}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}),\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\right\rangle+\sum_{j\in\mathcal{V}}\left\langle\lambda_{V}^{(j)}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}),\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\right\rangle$
$\displaystyle=\sum_{j\in\mathcal{V}}\sum_{i\in\mathcal{U}}r_{i,j}(\mathbf{u}_{i}^{\top}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}))^{2}+\alpha_{0}\sum_{j\in\mathcal{V}}\sum_{i\in\mathcal{U}}(\mathbf{u}_{i}^{\top}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}))^{2}+\lambda_{ex}\sum_{j\in\mathcal{V}}(\mathbf{s}^{\top}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}))^{2}+\sum_{j\in\mathcal{V}}\|\lambda_{V}^{(j)}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime})\|_{2}^{2}$
$\displaystyle\geq\underline{\lambda}_{V}\sum_{j\in\mathcal{V}}\|\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\|_{2}^{2}=\underline{\lambda}_{V}\|\mathbf{V}-\mathbf{V}^{\prime}\|_{F}^{2},$

where
$\underline{\lambda}_{V}\coloneqq\min_{j\in\mathcal{V}}\lambda_{V}^{(j)}$.
Thus, the function $g$ is $\underline{\lambda}_{V}$-strongly convex with
respect to $\mathbf{V}$. We also have

$\displaystyle
g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})-g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})$
$\displaystyle\leq\langle\nabla_{\mathbf{V}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k}),\mathbf{V}^{k+1}-\mathbf{V}^{k}\rangle_{F}-\frac{\underline{\lambda}_{V}}{2}\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2}$
(10)
$\displaystyle=-\frac{\underline{\lambda}_{V}}{2}\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2},$

where the last equality follows from the fact that $\mathbf{V}^{k+1}$
minimizes $g(\mathbf{V},\mathbf{U}^{k},\mathbf{s}^{k})$; hence
$\nabla_{\mathbf{V}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})=0$
holds. Moreover, since $\mathbf{U}^{k+1}$ minimizes
$\frac{\rho}{2}\left\|\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U})^{\top}\mathds{1}-\mathbf{s}^{k}\right\|_{2}^{2}-\frac{\rho}{2}\|\mathbf{w}^{k}\|_{2}^{2}+\frac{1}{2\gamma}\|\mathbf{U}-\mathbf{U}^{k}\|_{F}^{2}+\langle\mathbf{U}-\mathbf{U}^{k},\nabla_{\mathbf{U}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})\rangle_{F}$,
we have:

$\displaystyle\frac{\rho}{2}\left\|\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}^{k}\right\|_{2}^{2}-\frac{\rho}{2}\|\mathbf{w}^{k}\|_{2}^{2}+\frac{1}{2\gamma}\|\mathbf{U}^{k+1}-\mathbf{U}^{k}\|_{F}^{2}+\langle\mathbf{U}^{k+1}-\mathbf{U}^{k},\nabla_{\mathbf{U}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})\rangle_{F}$
(11)
$\displaystyle\leq\frac{\rho}{2}\left\|\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}-\mathbf{s}^{k}\right\|_{2}^{2}-\frac{\rho}{2}\|\mathbf{w}^{k}\|_{2}^{2}.$

By combining Sections A.2, A.2 and A.2, we obtain:

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k},\mathbf{w}^{k})-L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
(12) $\displaystyle\leq
g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k})-g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})-\langle\mathbf{U}^{k+1}-\mathbf{U}^{k},\nabla_{\mathbf{U}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})\rangle_{F}-\frac{1}{2\gamma}\|\mathbf{U}^{k+1}-\mathbf{U}^{k}\|_{F}^{2}-\frac{\underline{\lambda}_{V}}{2}\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2}.$

On the other hand, under the assumption in Theorem 3.1, from Lemma A.1, the
function $g(\mathbf{V}^{k+1},\mathbf{U},\mathbf{s}^{k})$ is
$\sqrt{|\mathcal{U}|}\left((1+\alpha_{0})C_{V}+\bar{\lambda}_{U}\right)$-smooth
with respect to $\mathbf{U}$. Then, for any $\mathbf{U},\mathbf{U}^{\prime}$,
we have:

(13) $\displaystyle
g(\mathbf{V}^{k+1},\mathbf{U}^{\prime},\mathbf{s}^{k})-g(\mathbf{V}^{k+1},\mathbf{U},\mathbf{s}^{k})-\langle\nabla_{\mathbf{U}}g(\mathbf{V}^{k+1},\mathbf{U},\mathbf{s}^{k}),\mathbf{U}^{\prime}-\mathbf{U}\rangle_{F}$
$\displaystyle\leq\frac{\sqrt{|\mathcal{U}|}((1+\alpha_{0})C_{V}+\bar{\lambda}_{U})}{2}\|\mathbf{U}-\mathbf{U}^{\prime}\|_{F}^{2}.$

By combining Section A.2 and Eq. 13, we obtain the following inequality:

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k},\mathbf{w}^{k})-L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
$\displaystyle\leq\frac{\sqrt{|\mathcal{U}|}((1+\alpha_{0})C_{V}+\bar{\lambda}_{U})-1/\gamma}{2}\|\mathbf{U}^{k+1}-\mathbf{U}^{k}\|_{F}^{2}-\frac{\underline{\lambda}_{V}}{2}\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2}.$

∎

###  A.3. Proof of Lemma A.3

###### Proof.

Let us define $h^{k}(\mathbf{s})\coloneqq
g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s})+\frac{\rho}{2}\left\|\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}\right\|_{2}^{2}-\frac{\rho}{2}\|\mathbf{w}^{k}\|_{2}^{2}$.
We have:

$\displaystyle\langle\nabla h^{k}(\mathbf{s})-\nabla
h^{k}(\mathbf{s}^{\prime}),\mathbf{s}-\mathbf{s}^{\prime}\rangle$
$\displaystyle=\left\langle\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s})-\rho\left(\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}\right)-\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{\prime})+\rho\left(\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}^{\prime}\right),\mathbf{s}-\mathbf{s}^{\prime}\right\rangle$
$\displaystyle=\left\langle\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s})-\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{\prime})+\rho(\mathbf{s}-\mathbf{s}^{\prime}),\mathbf{s}-\mathbf{s}^{\prime}\right\rangle$
$\displaystyle\geq\rho\|\mathbf{s}-\mathbf{s}^{\prime}\|_{2}^{2},$

where the inequality follows from the convexity of
$g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\cdot)$. Consequently, $h^{k}$ is a
$\rho$-strongly convex function. Therefore, we have:

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k})-L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k},\mathbf{w}^{k})$
$\displaystyle=h^{k}(\mathbf{s}^{k+1})-h^{k}(\mathbf{s}^{k})\leq\langle\nabla
h^{k}(\mathbf{s}^{k+1}),\mathbf{s}^{k+1}-\mathbf{s}^{k}\rangle-\frac{\rho}{2}\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}^{2}=-\frac{\rho}{2}\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}^{2},$

where the last equality follows from that $\mathbf{s}^{k+1}$ minimizes
$h^{k}(\mathbf{s})$, i.e., $\nabla h^{k}(\mathbf{s}^{k+1})=0$. ∎

###  A.4. Proof of Lemma A.4

###### Proof.

From the definition of
$L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w})$ and
the update rule of $\mathbf{w}^{k}$, we have:

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k+1})-L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k})$
$\displaystyle=\frac{\rho}{2}\left\|\mathbf{w}^{k+1}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}^{k+1}\right\|_{2}^{2}-\frac{\rho}{2}\|\mathbf{w}^{k+1}\|_{2}^{2}-\frac{\rho}{2}\left\|\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}^{k+1}\right\|_{2}^{2}+\frac{\rho}{2}\|\mathbf{w}^{k}\|_{2}^{2}$
(14)
$\displaystyle=\rho\langle\mathbf{w}^{k+1}-\mathbf{w}^{k},\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}^{k+1}\rangle=\rho\|\mathbf{w}^{k+1}-\mathbf{w}^{k}\|_{2}^{2}.$

On the other hand, since $\mathbf{s}^{k+1}$ minimizes the convex function
$h^{k}(\mathbf{s})$, the first-order optimality condition implies:

$\displaystyle\nabla h^{k}(\mathbf{s}^{k+1})$
$\displaystyle=\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1})-\rho\left(\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k+1})^{\top}\mathds{1}-\mathbf{s}^{k+1}\right)$
$\displaystyle=\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1})-\rho\mathbf{w}^{k+1}=0.$

Thus, it holds that

(15)
$\displaystyle\mathbf{w}^{k+1}=\frac{1}{\rho}\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1}).$

By combining Section A.4, Eq. 15, and Lemma A.1, we obtain:

$\displaystyle
L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k+1})-L_{\rho}(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1},\mathbf{w}^{k})$
$\displaystyle=\frac{1}{\rho}\|\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1})-\nabla_{\mathbf{s}}g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})\|_{2}^{2}$
$\displaystyle=\frac{1}{\rho}\|\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1})-\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})+\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})-\nabla_{\mathbf{s}}g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})\|_{2}^{2}$
$\displaystyle\leq\frac{1}{\rho}\Bigl{(}\|\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k+1},\mathbf{s}^{k+1})-\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})\|_{2}+\|\nabla_{\mathbf{s}}g(\mathbf{V}^{k+1},\mathbf{U}^{k},\mathbf{s}^{k})-\nabla_{\mathbf{s}}g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})\|_{2}\Bigr{)}^{2}$
$\displaystyle\leq\frac{1}{\rho}\left(\lambda_{ex}\|\mathbf{V}^{k+1}\|_{F}^{2}\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}+\lambda_{ex}\left(\|\mathbf{V}^{k+1}\|_{F}+\|\mathbf{V}^{k}\|_{F}\right)\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}\|\mathbf{s}^{k}\|_{2}\right)^{2}$
$\displaystyle\leq\frac{3}{\rho}\left(\lambda_{ex}^{2}\|\mathbf{V}^{k+1}\|_{F}^{4}\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}^{2}+\lambda_{ex}^{2}\left(\|\mathbf{V}^{k+1}\|_{F}^{2}+\|\mathbf{V}^{k}\|_{F}^{2}\right)\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2}\|\mathbf{s}^{k}\|_{2}^{2}\right)$
$\displaystyle\leq\frac{3\lambda_{ex}^{2}}{\rho}\left(C_{V}^{2}\|\mathbf{s}^{k+1}-\mathbf{s}^{k}\|_{2}^{2}+2C_{V}C_{s}\|\mathbf{V}^{k+1}-\mathbf{V}^{k}\|_{F}^{2}\right),$

where the third inequality follows from $(a+b+c)^{2}\leq 3(a^{2}+b^{2}+c^{2})$
for $a,b,c\in\mathbb{R}$. ∎

###  A.5. Proof of Lemma A.5

###### Proof of Lemma A.5.

Under the assumption in Theorem 3.1, the function
$g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s})$ is $\lambda_{ex}C_{V}$-smooth
with respect to $\mathbf{s}$ from Lemma A.1, and then we have, for any
$\mathbf{s},\mathbf{s}^{\prime}$,

(16) $\displaystyle
g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{\prime})-g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s})-\langle\nabla_{\mathbf{s}}g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}),\mathbf{s}^{\prime}-\mathbf{s}\rangle\leq\frac{\lambda_{ex}C_{V}}{2}\|\mathbf{s}-\mathbf{s}^{\prime}\|_{2}^{2}.$

By combining Eq. 15 and Eq. 16, we obtain:

$\displaystyle
L_{\rho}(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k},\mathbf{w}^{k})$
$\displaystyle=g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})+\frac{\rho}{2}\left\|\mathbf{w}^{k}+\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}-\mathbf{s}^{k}\right\|_{2}^{2}-\frac{\rho}{2}\|\mathbf{w}^{k}\|_{2}^{2}$
$\displaystyle=g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})+\rho\langle\mathbf{w}^{k},\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}-\mathbf{s}^{k}\rangle+\frac{\rho}{2}\left\|\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}-\mathbf{s}^{k}\right\|_{2}^{2}$
$\displaystyle=g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k})-\langle\nabla_{\mathbf{s}}g(\mathbf{V}^{k},\mathbf{U}^{k},\mathbf{s}^{k}),\mathbf{s}^{k}-\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}\rangle+\frac{\rho}{2}\left\|\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}-\mathbf{s}^{k}\right\|_{2}^{2}$
$\displaystyle\geq
g(\mathbf{V}^{k},\mathbf{U}^{k},\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1})+\frac{\rho-\lambda_{ex}C_{V}}{2}\left\|\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}-\mathbf{s}^{k}\right\|_{2}^{2}$
$\displaystyle\geq\frac{\rho-\lambda_{ex}C_{V}}{2}\left\|\frac{1}{|\mathcal{U}|}(\mathbf{U}^{k})^{\top}\mathds{1}-\mathbf{s}^{k}\right\|_{2}^{2},$

where the last inequality follows from
$g(\mathbf{V},\mathbf{U},\mathbf{s})\geq 0$ for any $\mathbf{v}$,
$\mathbf{U}$, and $\mathbf{s}$. ∎

###  A.6. Proof for Lemma A.1

###### Proof.

For fixed $\mathbf{U},\mathbf{s}$, for all $\mathbf{V},\mathbf{V}^{\prime}$,
we have the following

$\displaystyle\|\nabla_{\mathbf{V}}g(\mathbf{V},\mathbf{U},\mathbf{s})-\nabla_{\mathbf{V}}g(\mathbf{V}^{\prime},\mathbf{U},\mathbf{s})\|_{F}$
$\displaystyle=\sum_{j\in\mathcal{V}}\left\|\left(\sum_{i\in\mathcal{U}}r_{i,j}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}+\alpha_{0}\mathbf{G}_{U}+\lambda_{ex}\mathbf{s}\mathbf{s}^{\top}+\lambda_{V}^{(j)}\mathbf{I}\right)(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime})\right\|_{2}$
$\displaystyle\leq\sum_{j\in\mathcal{V}}\sum_{i\in\mathcal{U}}\|r_{i,j}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime})\|_{2}+\sum_{j\in\mathcal{V}}\sum_{i\in\mathcal{U}}\|\alpha_{0}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime})\|_{2}+\sum_{j\in\mathcal{V}}\|\lambda_{ex}\mathbf{s}\mathbf{s}^{\top}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime})\|_{2}+\sum_{j\in\mathcal{V}}\|\lambda_{V}^{(j)}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime})\|_{2}$
$\displaystyle=\sum_{j\in\mathcal{V}}\sum_{i\in\mathcal{U}}r_{i,j}|\mathbf{u}_{i}^{\top}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime})|\cdot\|\mathbf{u}_{i}\|_{2}+\alpha_{0}\sum_{j\in\mathcal{V}}\sum_{i\in\mathcal{U}}|\mathbf{u}_{i}^{\top}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime})|\cdot\|\mathbf{u}_{i}\|_{2}+\lambda_{ex}\sum_{j\in\mathcal{V}}|\mathbf{s}^{\top}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime})|\cdot\|\mathbf{s}\|_{2}+\sum_{j\in\mathcal{V}}\|\lambda_{V}^{(j)}(\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime})\|_{2}$
$\displaystyle\leq\sum_{j\in\mathcal{V}}\sum_{i\in\mathcal{U}}r_{i,j}\|\mathbf{u}_{i}\|_{2}^{2}\|\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\|_{2}+\alpha_{0}\sum_{j\in\mathcal{V}}\sum_{i\in\mathcal{U}}\|\mathbf{u}_{i}\|_{2}^{2}\|\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\|_{2}+\lambda_{ex}\sum_{j\in\mathcal{V}}\|\mathbf{s}\|_{2}^{2}\|\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\|_{2}+\bar{\lambda}_{V}\sum_{j\in\mathcal{V}}\|\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\|_{2}$
$\displaystyle\leq\left(\left(1+\alpha_{0}\right)\sum_{i\in\mathcal{U}}\|\mathbf{u}_{i}\|_{2}^{2}+\lambda_{ex}\|\mathbf{s}\|_{2}^{2}+\bar{\lambda}_{V}\right)\sum_{j\in\mathcal{V}}\|\mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\|_{2}$
$\displaystyle\leq\sqrt{|\mathcal{V}|}\left(\left(1+\alpha_{0}\right)\|\mathbf{U}\|_{F}^{2}+\lambda_{ex}\|\mathbf{s}\|_{2}^{2}+\bar{\lambda}_{V}\right)\|\mathbf{V}-\mathbf{V}^{\prime}\|_{F},$

where $\bar{\lambda}_{V}\coloneqq\max_{j\in\mathcal{V}}\lambda_{V}^{(j)}$.
Note that the second inequality follows from the Cauchy–Schwarz inequality.

In addition, we have, for a fixed $\mathbf{V}$, for all
$\mathbf{s},\mathbf{s}^{\prime}$ and $\mathbf{U},\mathbf{U}^{\prime}$,

$\displaystyle\|\nabla_{\mathbf{U}}g(\mathbf{V},\mathbf{U},\mathbf{s})-\nabla_{\mathbf{U}}g(\mathbf{V},\mathbf{U}^{\prime},\mathbf{s}^{\prime})\|_{F}$
$\displaystyle=\sum_{i\in\mathcal{U}}\Biggl{\|}\left(\sum_{j\in\mathcal{V}}r_{i,j}\mathbf{v}_{j}\mathbf{v}_{j}^{\top}+\alpha_{0}\sum_{j\in\mathcal{V}}\mathbf{v}_{j}\mathbf{v}_{j}^{\top}+\lambda_{U}^{(i)}\mathbf{I}\right)\mathbf{u}_{i}-\sum_{j\in\mathcal{V}}r_{i,j}\mathbf{v}_{j}-\left(\sum_{j\in\mathcal{V}}r_{i,j}\mathbf{v}_{j}\mathbf{v}_{j}^{\top}+\alpha_{0}\sum_{j\in\mathcal{V}}\mathbf{v}_{j}\mathbf{v}_{j}^{\top}+\lambda_{U}^{(i)}\mathbf{I}\right)\mathbf{u}_{i}^{\prime}+\sum_{j\in\mathcal{V}}r_{i,j}\mathbf{v}_{j}\Biggr{\|}_{2}$
$\displaystyle=\sum_{i\in\mathcal{U}}\left\|\left(\sum_{j\in\mathcal{V}}r_{i,j}\mathbf{v}_{j}\mathbf{v}_{j}^{\top}+\alpha_{0}\sum_{j\in\mathcal{V}}\mathbf{v}_{j}\mathbf{v}_{j}^{\top}+\lambda_{U}^{(i)}\mathbf{I}\right)(\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime})\right\|_{2}$
$\displaystyle\leq\sum_{i\in\mathcal{U}}\sum_{j\in\mathcal{V}}\left\|\left(r_{i,j}\mathbf{v}_{j}\mathbf{v}_{j}^{\top}\right)(\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime})\right\|_{2}+\sum_{i\in\mathcal{U}}\sum_{j\in\mathcal{V}}\left\|\left(\alpha_{0}\mathbf{v}_{j}\mathbf{v}_{j}^{\top}\right)(\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime})\right\|+\sum_{i\in\mathcal{U}}\left\|\lambda_{U}^{(i)}(\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime})\right\|_{2}$
$\displaystyle=\sum_{i\in\mathcal{U}}\sum_{j\in\mathcal{V}}r_{i,j}\left\|\mathbf{v}_{j}\mathbf{v}_{j}^{\top}(\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime})\right\|_{2}+\alpha_{0}\sum_{i\in\mathcal{U}}\sum_{j\in\mathcal{V}}\left\|\mathbf{v}_{j}\mathbf{v}_{j}^{\top}(\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime})\right\|_{2}+\sum_{i\in\mathcal{U}}\left\|\lambda_{U}^{(i)}(\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime})\right\|_{2}$
$\displaystyle=\sum_{i\in\mathcal{U}}\sum_{j\in\mathcal{V}}r_{i,j}|\mathbf{v}_{j}^{\top}(\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime})|\left\|\mathbf{v}_{j}\right\|_{2}+\alpha_{0}\sum_{i\in\mathcal{U}}\sum_{j\in\mathcal{V}}|\mathbf{v}_{j}^{\top}(\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime})|\left\|\mathbf{v}_{j}\right\|_{2}+\sum_{i\in\mathcal{U}}\left\|\lambda_{U}^{(i)}(\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime})\right\|_{2}$
$\displaystyle\leq\sum_{i\in\mathcal{U}}\sum_{j\in\mathcal{V}}r_{i,j}\left\|\mathbf{v}_{j}\right\|_{2}^{2}\left\|\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime}\right\|_{2}+\alpha_{0}\sum_{i\in\mathcal{U}}\sum_{j\in\mathcal{V}}\left\|\mathbf{v}_{j}\right\|^{2}\left\|\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime}\right\|_{2}+\sum_{i\in\mathcal{U}}\left\|\lambda_{U}^{(i)}(\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime})\right\|_{2}$
$\displaystyle\leq(1+\alpha_{0})\left(\sum_{j\in\mathcal{V}}\|\mathbf{v}_{j}\|_{2}^{2}\right)\sum_{i\in\mathcal{U}}\|\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime}\|_{2}+\bar{\lambda}_{U}\sum_{i\in\mathcal{U}}\|\mathbf{u}_{i}-\mathbf{u}_{i}^{\prime}\|_{2}$
$\displaystyle\leq\sqrt{|\mathcal{U}|}\left((1+\alpha_{0})\|\mathbf{V}\|_{F}^{2}+\bar{\lambda}_{U}\right)\|\mathbf{U}-\mathbf{U}^{\prime}\|_{F},$

where the second inequality follows from the Cauchy–Schwarz inequality.

Here, for a fixed $\mathbf{V}$, for all $\mathbf{s},\mathbf{s}^{\prime}$ and
$\mathbf{U},\mathbf{U}^{\prime}$, we have:

$\displaystyle\|\nabla_{\mathbf{s}}g(\mathbf{V},\mathbf{U},\mathbf{s})-\nabla_{\mathbf{s}}g(\mathbf{V},\mathbf{U}^{\prime},\mathbf{s}^{\prime})\|_{2}$
$\displaystyle=\left\|\lambda_{ex}\mathbf{V}^{\top}\mathbf{V}\mathbf{s}-\lambda_{ex}\mathbf{V}^{\top}\mathbf{V}\mathbf{s}^{\prime}\right\|_{2}$
$\displaystyle=\lambda_{ex}\|\mathbf{V}^{\top}\mathbf{V}(\mathbf{s}-\mathbf{s}^{\prime})\|_{2}$
$\displaystyle\leq\lambda_{ex}\|\mathbf{V}^{\top}\mathbf{V}\|_{F}\|\mathbf{s}-\mathbf{s}^{\prime}\|_{2}$
$\displaystyle\leq\lambda_{ex}\|\mathbf{V}\|_{F}^{2}\|\mathbf{s}-\mathbf{s}^{\prime}\|_{2},$

where the first/second inequality follows from the Cauchy–Schwarz inequality.

Finally, for fixed $\mathbf{U}$ and $\mathbf{s}$, for all
$\mathbf{V},\mathbf{V}^{\prime}$, we have:

$\displaystyle\|\nabla_{\mathbf{s}}g(\mathbf{V},\mathbf{U},\mathbf{s})-\nabla_{\mathbf{s}}g(\mathbf{V}^{\prime},\mathbf{U},\mathbf{s})\|_{2}$
$\displaystyle=\lambda_{ex}\|\mathbf{V}^{\top}\mathbf{V}\mathbf{s}-\mathbf{V}^{\prime\top}\mathbf{V}^{\prime}\mathbf{s}\|_{2}$
$\displaystyle=\lambda_{ex}\|(\mathbf{V}^{\top}(\mathbf{V}-\mathbf{V}^{\prime})+(\mathbf{V}-\mathbf{V}^{\prime})^{\top}\mathbf{V}^{\prime})\mathbf{s}\|_{2}$
$\displaystyle\leq\lambda_{ex}\|\mathbf{V}^{\top}(\mathbf{V}-\mathbf{V}^{\prime})\mathbf{s}\|_{2}+\lambda_{ex}\|(\mathbf{V}-\mathbf{V}^{\prime})^{\top}\mathbf{V}^{\prime}\mathbf{s}\|_{2}$
$\displaystyle\leq\lambda_{ex}\|\mathbf{V}^{\top}(\mathbf{V}-\mathbf{V}^{\prime})\|_{F}\cdot\|\mathbf{s}\|_{2}+\lambda_{ex}\|(\mathbf{V}-\mathbf{V}^{\prime})^{\top}\mathbf{V}^{\prime}\|_{F}\cdot\|\mathbf{s}\|_{2}$
$\displaystyle\leq\lambda_{ex}(\|\mathbf{V}\|_{F}+\|\mathbf{V}^{\prime}\|_{F})\cdot\|\mathbf{s}\|_{2}\|\mathbf{V}-\mathbf{V}^{\prime}\|_{F},$

where the second/third inequality follows from the Cauchy–Schwarz inequality.
∎