Towards Bursting Filter Bubble\\via Contextual Risks and Uncertainties
Rikiya Takahashi
SmartNews
Inc.
6-25-16
Jingumae
Shibuya-ku
Tokyo 150-0001
Japan
Rikiya.Takahashi@gmail.com
A rising topic in computational journalism is how to enhance the diversity in news served to subscribers to foster exploration behavior in news reading. Despite the success of preference learning in personalized news recommendation, their over-exploitation causes filter bubble that isolates readers from opposing viewpoints and hurts long-term user experiences with lack of serendipity. Since news providers can recommend neither opposite nor diversified opinions if unpopularity of these articles is surely predicted, they can only bet on the articles whose forecasts of click-through rate involve high variability (risks) or high estimation errors (uncertainties). We propose a novel Bayesian model of uncertainty-aware scoring and ranking for news articles. The Bayesian binary classifier models probability of success (defined as a news click) as a Beta-distributed random variable conditional on a vector of the context (user features, article features, and other contextual features). The posterior of the contextual coefficients can be computed efficiently using a low-rank version of Laplace’s method via thin Singular Value Decomposition. Efficiencies in personalized targeting of exceptional articles, which are chosen by each subscriber in test period, are evaluated on real-world news datasets. The proposed estimator slightly outperformed existing training and scoring algorithms, in terms of efficiency in identifying successful outliers. %Valuable news often accompanies serendipity and serving mind-changing articles %to news subscribers is a rising and essential issue in computational journalism. %Despite the success of preference learning in personalized news recommendation, %their over-exploitation in terms of exploration-exploitation trade-off %causes filter bubble that isolates readers from opposing viewpoints. %Since news providers can recommend neither opposite nor diversified opinions %as long as unpopularity of these articles is surely predicted, %they can bet only on the articles whose forecasts of click-through rate involve high fundamental variabilities (risks) and/or %high estimation errors (uncertainties). %Being inspired by an insight that %estimation error of variability statistic is higher for some contexts associated with higher fundamental variability, %the authors newly propose a Bayesian binary classifier whose probability of success is %a Beta-distributed random variable conditional on vector of context. %A low-rank version of Laplace's method, via thin Singular Value Decomposition, %is introduced for yielding both efficient Bayesian training and test recommendation, %where the recommendation is done via novel uncertainty-aware scoring measures derived from the approximate posterior. %Efficiencies in personalized targeting of exceptional articles, which are chosen by each subscriber in test period %despite the dissimilarity to the chosen ones in training period, %are evaluated for real-world news datasets. %The proposed estimator experimentally outperformed existing training and scoring algorithms, %in terms of efficiency in identifying the successful outliers.
anti_filter_bubble.tex
\documentclass{article}
\pdfoutput=1

% if you need to pass options to natbib, use, e.g.:
% \PassOptionsToPackage{numbers, compress}{natbib}
% before loading nips_2017
%
% to avoid loading the natbib package, add option nonatbib:
% \usepackage[nonatbib]{nips_2017}

\usepackage[final,nonatbib]{arxiv_2017}

% to compile a camera-ready version, add the [final] option, e.g.:
%\usepackage[final]{nips_2017}

\usepackage[utf8]{inputenc} % allow utf-8 input
\usepackage[T1]{fontenc}    % use 8-bit T1 fonts
\usepackage{hyperref}       % hyperlinks
\usepackage{url}            % simple URL typesetting
\usepackage{booktabs}       % professional-quality tables
\usepackage{amsfonts}       % blackboard math symbols
\usepackage{nicefrac}       % compact symbols for 1/2, etc.
\usepackage{microtype}      % microtypography
\usepackage{url,amsmath,amssymb,amsthm,braket,algorithmic,algorithm}
%\usepackage[dvipdfmx]{graphicx}
\usepackage{graphicx}
%\usepackage{authblk}

\title{Towards Bursting Filter Bubble\\via Contextual Risks and Uncertainties}

% The \author macro works with any number of authors. There are two
% commands used to separate the names and addresses of multiple
% authors: \And and \AND.
%
% Using \And between authors leaves it to LaTeX to determine where to
% break the lines. Using \AND forces a line break at that point. So,
% if LaTeX puts 3 of 4 authors names on the first line, and the last
% on the second line, try using \AND instead of \And before the third
% author name.

%\author[1]{Rikiya Takahashi}
%\author[2]{Shunan Zhang}
%\affil[1]{SmartNews, Inc.}
%\affil[2]{Smartnews International, Inc.}
\author{
  Rikiya Takahashi\\
  SmartNews, Inc.\\
  6-25-16, Jingumae, Shibuya-ku\\
  Tokyo 150-0001, Japan\\
  \texttt{Rikiya.Takahashi@gmail.com} \\
  %% examples of more authors
 \And
  Shunan Zhang\\
  Smartnews International, Inc.\\
  144 2nd St, San Francisco\\
  CA 94105, USA\\
  \texttt{shunan.z@gmail.com} \\
}
\begin{document}
% \nipsfinalcopy is no longer used

\maketitle



\begin{abstract}
A rising topic in computational journalism is how to enhance the diversity in news served to subscribers to foster exploration behavior in news reading. Despite the success of preference learning in personalized news recommendation, their over-exploitation causes filter bubble that isolates readers from opposing viewpoints and hurts long-term user experiences with lack of serendipity. Since news providers can recommend neither opposite nor diversified opinions if unpopularity of these articles is surely predicted, they can only bet on the articles whose forecasts of click-through rate involve high variability (risks) or high estimation errors (uncertainties). We propose a novel Bayesian model of 
uncertainty-aware scoring and ranking for news articles. The Bayesian binary classifier models probability of success (defined as a news click) as a Beta-distributed random variable conditional on a vector of the context (user features, article features, and other contextual features). The posterior of the contextual coefficients can be computed efficiently using a low-rank version of Laplace’s method via thin Singular Value Decomposition. Efficiencies in personalized targeting of exceptional articles, which are chosen by each subscriber in test period, are evaluated on real-world news datasets. The proposed estimator slightly outperformed existing training and scoring algorithms, in terms of efficiency in identifying successful outliers.
%Valuable news often accompanies serendipity and serving mind-changing articles 
%to news subscribers is a rising and essential issue in computational journalism.
%Despite the success of preference learning in personalized news recommendation,
%their over-exploitation in terms of exploration-exploitation trade-off
%causes filter bubble that isolates readers from opposing viewpoints.
%Since news providers can recommend neither opposite nor diversified opinions 
%as long as unpopularity of these articles is surely predicted,
%they can bet only on the articles whose forecasts of click-through rate involve high fundamental variabilities (risks) and/or %high estimation errors (uncertainties).
%Being inspired by an insight that
%estimation error of variability statistic is higher for some contexts associated with higher fundamental variability,
%the authors newly propose a Bayesian binary classifier whose probability of success is
%a Beta-distributed random variable conditional on vector of context.
%A low-rank version of Laplace's method, via thin Singular Value Decomposition,
%is introduced for yielding both efficient Bayesian training and test recommendation,
%where the recommendation is done via novel uncertainty-aware scoring measures derived from the approximate posterior.
%Efficiencies in personalized targeting of exceptional articles, which are chosen by each subscriber in test period
%despite the dissimilarity to the chosen ones in training period,
%are evaluated for real-world news datasets.
%The proposed estimator experimentally outperformed existing training and scoring algorithms,
%in terms of efficiency in identifying the successful outliers.
\end{abstract}


\section{Introduction}
\label{sec:intro}

% Msg: personalized news recommendation has been successful
Personalized news recommendation has been a successful application of machine learning
although it would also endanger people's open-mindedness towards opposing viewpoints.
Targeted news delivery has been reinforced by collections of click response logs, 
with preference learning algorithms such as
bilinear models \cite{chu2009personalized,sharma2015feature},
topic models \cite{wu2010topic},
matrix factorization \cite{das2007google,agarwal2010matrix,gunasekar2015consistent}, learning to rank \cite{dali2010learning},
or their hybrid 
%\cite{claypool1999combining,liu2010personalized,zheng2013penetrate}.
\cite{claypool1999combining,zheng2013penetrate}.
In terms of exploration-exploitation trade-off, however, the exploitative nature of data-oriented personalization
has been criticized as a cause of the filter bubble \cite{pariser2011filter},
which is a compound problem fueled both by narrow-minded human's willing choice \cite{bakshy2015exposure}
and by machines over-confident about their scoring.
In this paper, we address the problem of over-confident machines while psychologically opening people's mind is out of the current scope.

% Msg: it is in vain to recommend articles whose unpopularity is statistically significant.
Even by diversifying recommendations or by exploring around uncertain preferences,
online news providers cannot deliver articles whose unpopularity is surely predicted.
Submodular maximization algorithms \cite{yue2011linear} or 
Determinantal Point Processes (DPPs \cite{affandi2012markov,gillenwater2014em,ren2015summarizing}; see \cite{kulesza2013determinantal} for good summary)
generate a subset of input items associated with univariate scores (e.g., estimated click-through rates).
The output subset consists of some high-score items that are well diversified.
Unfortunately, due to the exclusion of low-score items,
diversification algorithms never select items whose unpopularity is statistically significant.
Another direction towards serendipity is to use not unbiased but optimistically-biased estimates
as in contextual multi-armed bandit algorithms
\cite{liu2010personalized,li2010contextual,yue2011linear,may2012optimistic,tang2014ensemble}.
Bandit algorithms, however, also fail to select surely unpopular items,
because expected cumulative reward is their main objective in repeated trials.
As a consequence, as long as a large number of total clicks must be retained,
news providers cannot bet against readers' preference evidenced by large samples.

% Msg: philosophy and underlying deliberation of our model
We instead recommend articles whose popularity involves high variabilities (risks) and/or high estimation errors (uncertainties),
by introducing a new Bayesian classifier that explicitly formalizes the dependence of variabilities on recommendation contexts.
It is worth observing that Maximum A Posteriori (MAP) point estimation and logistic-loss models (e.g., factorization machines
\cite{rendle2012factorization,juan2016field} and deep neural networks \cite{covington2016deep}) are broadly used in real-world systems.
We hypothesize that content providers too much rely on the logistic-loss MAP estimates 
because of their successes in retrieving popular items, 
and they often forget the fact that substance of their model
is an erroneous estimate deviated from the true preference.
Furthermore, we gaze at a statistical insight that estimation error of variability statistic is particularly high when the true variability is high.
By performing Bayesian interval estimation and assuming the true variability to be a function of context,
we robustly quantify the total stochasticity with clear separation between the essential variability and estimation error.
These deliberated philosophies behind our model are expected to increase the chance of discovering exceptional items,
whose preference uncertainty has been underestimated by the existing models.

% Msg: Technical components of our model
Our key ideas to accurately quantify the context-dependent risks and uncertainties 
are conditional Beta distribution and low-rank Laplace approximation with affordable computational costs.
We model the probability of success as  a random variable to obey a Beta distribution,
whose two parameters are both functions of input vector.
In training, we approximate the posterior by a low-rank Gaussian distribution
whose estimates of the variance-covariance matrices are given by thin Singular Value Decomposition (SVD).
Thanks to this thin SVD,
we can perform both training and recommendation
efficiently by sparse high-dimensional matrix libraries.
We also derive a variety of uncertainty-aware scoring measures
over the closed forms of the approximate posterior, and evaluate the performance of every measure by using real-world news-reading datasets.
While we compare linear models in the experiments,
since our model is a robust extension of logistic regression with additional dispersion parameters,
the principle of combining conditional Beta distribution with low-rank Laplace's method
is broadly pluggable into many logistic-loss non-linear models.

The remainder of this paper is organized as follows.
Our Bayesian binary classifier and its approximate posterior inference are introduced in Section \ref{sec:training}.
We then derive several types of posterior-based scoring methods in Section \ref{sec:recommendation}.
Section \ref{sec:related_work} discusses the related work on robust classification
and human's exploration behavior though the latter theme is out of our scope.
In Section \ref{sec:experiments}, we numerically evaluate the proposed algorithm and scoring measures
for the news datasets, with introducing a novel performance indicator for serendipity-oriented targeting.
Section \ref{sec:conclusion} concludes the paper.


\section{Bayesian training of a binary classifier with contextual variabilities}
\label{sec:training}

Let us introduce a Bayesian binary classifier in which
the probability of success  is a random variable conditional on vector of context.
%Thanks to the modeling of context-dependent variabilities,
%the observable rates of success often widely vary across the contexts
%even when these contexts share the same values of input vector.
As shown in Section \ref{subsec:beta_binomial_model},
the frequency of success obeys a Beta-binomial distribution
and we place a Gaussian prior.
Though the exact posterior is intractable, 
low-rank Laplace's method introduced in Section \ref{subsec:laplace_approximation}
provides a closed-form estimate with affordable costs in matrix computation,
and also leads derivation of an approximate marginal likelihood in Section \ref{subsec:hyperparameter_section}.
The low-rank approximate posterior on the maximum-marginal-likelihood prior hyperparameters
is finally used for recommendation whose details are later provided in Section \ref{sec:recommendation}.

\subsection{Beta-binomial model and Gaussian prior}
\label{subsec:beta_binomial_model}

Let $\theta_{ij}\!\in\![0, 1]$ be probability with which user $i$ chooses item $j$, and 
let $\boldsymbol{x}_{ij} \in {\mathbb R}^d$ be vector of context 
when item $j$ is shown to user $i$.
In our news example, each item is an article and the vector $\boldsymbol{x}_{ij}$
is defined on a conjunction of news subscriber $i$'s characteristics and text content of article $j$.
We assume that probability $\theta_{ij}$ 
obeys a conditional Beta distribution on vector $\boldsymbol{x}_{ij}$ as
\begin{equation}
\label{equation:beta_model}
p(\theta_{ij}\vert\boldsymbol{x}_{ij}, \boldsymbol{w})\!=\!Be\!\left(\theta_{ij}; %
\exp((\boldsymbol{\beta}+\boldsymbol{\rho})^\top\boldsymbol{x}_{ij}), %
\exp(\boldsymbol{\rho}^\top\boldsymbol{x}_{ij})\right)\!,
\end{equation}
where $Be\!\left(\theta; \alpha_1, \alpha_2\right)\!\triangleq\!%
\frac{\Gamma(\alpha_1\!+\!\alpha_2)}{\Gamma(\alpha_1)\Gamma(\alpha_2)}%
\theta^{\alpha_1-1}(1\!-\!\theta)^{\alpha_2-1}$
and $\boldsymbol{w}\!\triangleq\!(\boldsymbol{\beta}^\top,\boldsymbol{\rho}^\top)^\top$
is a vector of regression coefficients. Let us define expectation operator ${\mathbb E}_{p}$ as 
${\mathbb E}_{p(\theta)}[f(\theta)]\!\triangleq\!\int_{\theta}f(\theta)p(\theta)d\theta$. 
Since ${\mathbb E}_{p(\theta_{ij}\vert\boldsymbol{x}_{ij}, \boldsymbol{w})}[\theta_{ij}]\!=\!1/(1+\exp(-\boldsymbol{\beta}^\top\boldsymbol{x}_{ij}))$,
Eq.~(\ref{equation:beta_model}) is interpreted as a logistic regression
with additional variabilities introduced by unobservable variables that are not contained in vector $\boldsymbol{x}_{ij}$.

We estimate posterior of $\boldsymbol{w}$.
Let ${\mathcal N}(\cdot; \boldsymbol{\mu}, \boldsymbol{\Sigma})$
be the probability density function of multivariate Gaussian distribution
whose mean is $\boldsymbol{\mu}$ and whose variance-covariance matrix 
is $\boldsymbol{\Sigma}$.
We place an isotropic Gaussian prior $p(\boldsymbol{w})={\mathcal N}(\boldsymbol{\beta}; \boldsymbol{0}_{d}, c_{\boldsymbol{\beta}}^{-1}\boldsymbol{I}_{d}){\mathcal N}(\boldsymbol{\rho}; \boldsymbol{0}_{d}, c_{\boldsymbol{\rho}}^{-1}\boldsymbol{I}_{d})$ where $\boldsymbol{0}_d$ and $\boldsymbol{I}_d$ are the $d$-dimensional zero vector and identity matrix, and $\boldsymbol{c}\!\triangleq\!(c_{\boldsymbol{\beta}}, c_{\boldsymbol{\rho}})^\top$ is a vector of prior hyperparameters.

Because we observe only choice frequencies,
the data likelihood is given by probability mass function of Beta-binomial distribution.
Let $n_{ij}$ and $v_{ij}$ be the frequencies with which user $i$ is exposed to item $j$
and with which user $i$ chooses item $j$, respectively. 
In online news recommendation, $n_{ij}$ and $v_{ij}$ are called 
the numbers of impressions and page views, respectively.
Since $v_{ij}$ is distributed from binomial distribution whose number of trials is $n_{ij}$
and whose probability of success is $\theta_{ij}$, we can obtain the data likelihood 
marginal over the latent choice probabilities as
\begin{equation}
%\label{equation:beta_binomial_model}
p(v_{ij}\vert n_{ij},\boldsymbol{x}_{ij},\boldsymbol{w})\!\triangleq\!%
{\mathbb E}_{p(\theta_{ij}\vert\boldsymbol{x}_{ij}, \boldsymbol{w})}\!\!\left[\theta_{ij}^{v_{ij}}\!(1\!-\!\theta_{ij})^{n_{ij}\!-\!v_{ij}}\!\right]%
\!\!=\!BB\!\left(v_{ij}; n_{ij},\exp(\!(\boldsymbol{\beta}\!+\!\boldsymbol{\rho})\!^\top\!\boldsymbol{x}_{ij}\!),%
\exp(\boldsymbol{\rho}\!^\top\!\!\boldsymbol{x}_{ij}\!)\right)\!\nonumber,
%&\equiv&BB\!\left(v_{ij}; n_{ij}, \exp\left(\boldsymbol{w}^\top\!%
%\begin{pmatrix}%
%S\boldsymbol{x}_{ij}\\%
%\boldsymbol{x}_{ij}%
%\end{pmatrix}%
%\right),%
%\exp\left(\boldsymbol{w}^\top\!%
%\begin{pmatrix}%
%\boldsymbol{0}_d\\%
%\boldsymbol{x}_{ij}%
%\end{pmatrix}%
%\right)\right),\nonumber
\end{equation}
where $BB\!\left(v; n, \alpha_1, \alpha_2\right)$ represents the
unnormalized probability mass function of Beta-binomial distribution
such that $BB\!\left(v; n, \alpha_1, \alpha_2\right)\!\triangleq\!%
\frac{\Gamma(\alpha_1\!+\!\alpha_2)}{\Gamma(\alpha_1\!+\!\alpha_2\!+\!n)}%
\frac{\Gamma(\alpha_1\!+\!v)}{\Gamma(\alpha_1)}%
\frac{\Gamma(\alpha_2\!+\!n\!-\!v)}{\Gamma(\alpha_2)}$.
Since Beta-binomial and Gaussian distributions are not conjugate, the true posterior is analytically intractable.
% and needs to be approximated.

\subsection{Block low-rank Laplace approximation}
\label{subsec:laplace_approximation}

Laplace's method is useful in approximate Bayesian inference particularly when the true posterior is
close to a Gaussian distribution. The approximate Gaussian posterior's 
mean is given by the Maximum A Posteriori (MAP) estimate and 
its variance-covariance matrix is obtained through a quadratic approximation around the MAP estimate. 
Our negative joint log-likelihood is given as ${\mathcal L}_{{\mathcal D},\boldsymbol{c}}(\boldsymbol{w})\!\triangleq\!%
-\sum_{(i,j)\in{\mathcal D}} \log BB\left(v_{ij}; n_{ij}, \exp((\boldsymbol{\beta}+\boldsymbol{\rho})^\top\boldsymbol{x}_{ij}),%
\exp(\boldsymbol{\rho}^\top\boldsymbol{x}_{ij})\right)%
+\frac{c_{\boldsymbol{\beta}}}{2}\Vert\boldsymbol{\beta}\Vert^2%
+\frac{c_{\boldsymbol{\rho}}}{2}\Vert\boldsymbol{\rho}\Vert^2$,
where  $\Vert\cdot\Vert$ denotes the $L_2$ norm
and ${\mathcal D}=\lbrace (i,j) \rbrace$ is the training set of user-item pairs.
%For short notations, 
Let
\begin{eqnarray}
\alpha_{ij}^{+}&\!\triangleq\!&\exp((\boldsymbol{\beta}\!+\!\boldsymbol{\rho})\!^\top\!\boldsymbol{x}_{ij}),~%
\alpha_{ij}^{-}\!\triangleq\!\exp(\boldsymbol{\rho}\!^\top\!\boldsymbol{x}_{ij}),~%
\alpha_{ij}\!\triangleq\!\alpha_{ij}^{+}+\alpha_{ij}^{-},\nonumber\\
g_{ij}^{(\boldsymbol{\beta})}&\!\triangleq\!&
-\alpha_{ij}^{+}\left[\Psi(\alpha_{ij}^{+}+v_{ij})\!-\!\Psi(\alpha_{ij}^{+})%
\!-\!\Psi(\alpha_{ij}+n_{ij})\!+\!\Psi(\alpha_{ij})\right]\nonumber,\\
g_{ij}^{(\boldsymbol{\rho})}&\!\triangleq\!&%
-\alpha_{ij}^{+}\left[\Psi(\alpha_{ij}^{+}+v_{ij})\!-\!\Psi(\alpha_{ij}^{+})\right]-%
\alpha_{ij}^{-}\left[\Psi(\alpha_{ij}^{-}+n_{ij}-v_{ij})\!-\!\Psi(\alpha_{ij}^{-})\right]%
\nonumber\\
&&+\alpha_{ij}\left[\Psi(\alpha_{ij}+n_{ij})\!-\!\Psi(\alpha_{ij})\right]\nonumber,\\
h_{ij}^{(\boldsymbol{\beta})}&\!\triangleq\!&%
-(\alpha_{ij}^{+})^2\left[\Psi'(\alpha_{ij}^{+}+v_{ij})\!-\!\Psi'(\alpha_{ij}^{+})%
\!-\!\Psi'(\alpha_{ij}+n_{ij})\!+\!\Psi'(\alpha_{ij})\right]+g_{ij}^{(\boldsymbol{\beta})}\nonumber,\\
h_{ij}^{(\boldsymbol{\beta},\boldsymbol{\rho})}&\!\triangleq\!&%
-(\alpha_{ij}^{+})^2\left[\Psi'(\alpha_{ij}^{+}+v_{ij})\!-\!\Psi'(\alpha_{ij}^{+})\right]%
+\alpha_{ij}\alpha_{ij}^{+}\left[\Psi'(\alpha_{ij}+n_{ij})\!-\!\Psi'(\alpha_{ij})\right]+g_{ij}^{(\boldsymbol{\beta})}\nonumber\mbox{, and}\\
h_{ij}^{(\boldsymbol{\beta})}&\!\triangleq\!&%
-(\alpha_{ij}^{+})^2\left[\Psi'(\alpha_{ij}^{+}+v_{ij})\!-\!\Psi'(\alpha_{ij}^{+})\right]%
-(\alpha_{ij}^{-})^2\left[\Psi'(\alpha_{ij}^{-}+n_{ij}-v_{ij})\!-\!\Psi'(\alpha_{ij}^{-})\right]\nonumber\\
&&+\alpha_{ij}^2\left[\Psi'(\alpha_{ij}+n_{ij})\!-\!\Psi'(\alpha_{ij})\right]+g_{ij}^{(\boldsymbol{\rho})}\nonumber,
\end{eqnarray}
where $\Psi(\cdot)$ and $\Psi'(\cdot)$ are the digamma and trigamma functions such that 
$\Psi(u)\!\triangleq\!\frac{\partial\log\Gamma(u)}{\partial u}$ and 
$\Psi'(u)\!\triangleq\!\frac{\partial\Psi(u)}{\partial u}$.
The MAP estimate 
$\widehat{\boldsymbol{w}}\!=\!\arg\min_{\boldsymbol{w}}{\mathcal L}_{{\mathcal D},\boldsymbol{c}}(\boldsymbol{w})$
is attained by descending the gradient 
$\frac{\partial{\mathcal L}_{{\mathcal D},\boldsymbol{c}}(\boldsymbol{w})}{\partial\boldsymbol{w}}%
\!=\!%
\sum_{(i,j)\in{\mathcal D}}%
(g_{ij}^{(\boldsymbol{\beta})}\boldsymbol{x}_{ij}^\top,%
g_{ij}^{(\boldsymbol{\rho})}\boldsymbol{x}_{ij}^\top)^\top%
\!+\!%
(c_{\boldsymbol{\beta}}\boldsymbol{\beta}^\top,%
c_{\boldsymbol{\rho}}\boldsymbol{\rho}^\top)^\top$
with the Hessian matrix
\iffalse
\begin{eqnarray}
\dfrac{\partial{\mathcal L}_{{\mathcal D},\boldsymbol{c}}(\boldsymbol{w})}{\partial\boldsymbol{w}}%
\!&=&\!%
\sum_{(i,j)\in{\mathcal D}}
\begin{pmatrix}%
g_{ij}^{(\boldsymbol{\beta})}\boldsymbol{x}_{ij}\\
g_{ij}^{(\boldsymbol{\rho})}\boldsymbol{x}_{ij}
\end{pmatrix}+%
\begin{pmatrix}%
c_{\boldsymbol{\beta}}\boldsymbol{\beta}\\
c_{\boldsymbol{\rho}}\boldsymbol{\rho}
\end{pmatrix}%
\mbox{ and }\nonumber\\
\label{equation:data_hessian}
\dfrac{\partial^2{\mathcal L}_{{\mathcal D},\boldsymbol{c}}(\boldsymbol{w})}{\partial\boldsymbol{w}\partial\boldsymbol{w}^\top}%
\!&=&\!\sum_{(i,j)\in{\mathcal D}}%
\begin{pmatrix}%
h_{ij}^{(\boldsymbol{\beta})}\boldsymbol{x}_{ij}\boldsymbol{x}_{ij}^\top &
h_{ij}^{(\boldsymbol{\beta},\boldsymbol{\rho})}\boldsymbol{x}_{ij}\boldsymbol{x}_{ij}^\top \\
h_{ij}^{(\boldsymbol{\beta},\boldsymbol{\rho})}\boldsymbol{x}_{ij}\boldsymbol{x}_{ij}^\top &
h_{ij}^{(\boldsymbol{\rho})}\boldsymbol{x}_{ij}\boldsymbol{x}_{ij}^\top
\end{pmatrix}+%
\begin{pmatrix}%
c_{\boldsymbol\beta}\boldsymbol{I}_d & \boldsymbol{O}_d\\
\boldsymbol{O}_d & c_{\boldsymbol\rho}\boldsymbol{I}_d
\end{pmatrix},%
\end{eqnarray}
\fi
\begin{equation}
\label{equation:data_hessian}
\dfrac{\partial^2{\mathcal L}_{{\mathcal D},\boldsymbol{c}}(\boldsymbol{w})}{\partial\boldsymbol{w}\partial\boldsymbol{w}^\top}%
~=\!\sum_{(i,j)\in{\mathcal D}}%
\!\begin{pmatrix}%
h_{ij}^{(\boldsymbol{\beta})}\boldsymbol{x}_{ij}\boldsymbol{x}_{ij}^\top &
h_{ij}^{(\boldsymbol{\beta},\boldsymbol{\rho})}\boldsymbol{x}_{ij}\boldsymbol{x}_{ij}^\top \\
h_{ij}^{(\boldsymbol{\beta},\boldsymbol{\rho})}\boldsymbol{x}_{ij}\boldsymbol{x}_{ij}^\top &
h_{ij}^{(\boldsymbol{\rho})}\boldsymbol{x}_{ij}\boldsymbol{x}_{ij}^\top
\end{pmatrix}+%
\begin{pmatrix}%
c_{\boldsymbol\beta}\boldsymbol{I}_d & \boldsymbol{O}_d\\
\boldsymbol{O}_d & c_{\boldsymbol\rho}\boldsymbol{I}_d
\end{pmatrix},%
\end{equation}
where $\boldsymbol{O}_d$ is the $d$-by-$d$ zero matrix.
After convergence,
we approximate the posterior by a factorial form
$q(\boldsymbol{w}\vert {\mathcal D},\boldsymbol{c})\!\triangleq\!%
{\mathcal N}\!\left(\!\boldsymbol{w}; \widehat{\boldsymbol{w}},\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{w}}\!\right)\!\equiv\!%
{\mathcal N}\!\left(\!\boldsymbol{\beta}; \widehat{\boldsymbol{\beta}},\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{\beta}}\!\right)%
{\mathcal N}\!\left(\!\boldsymbol{\rho}; \widehat{\boldsymbol{\rho}},\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{\rho}}\!\right)$%
%such that $\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{w}}\!\equiv\!
%\begin{pmatrix}%
%\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{\beta}}&\boldsymbol{O}_d\\
%\boldsymbol{O}_d&\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{\rho}}%
%\end{pmatrix}$
\footnote{
We avoid to parametrize $\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{w}}$ by authentic full variance-covariance matrix in this paper,
due to the too lengthy closed-form expression stemming from the complex correlation between $\boldsymbol{\beta}$ and $\boldsymbol{\rho}$.}.
%The matrix $\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{w}}$ can be computed from
%$\left(\left.\dfrac{\partial^2{\mathcal L}_{{\mathcal D},\boldsymbol{c}}(\boldsymbol{w})}{\partial\boldsymbol{w}\partial\boldsymbol{w}^\top}\right\vert_{\boldsymbol{w}\!=\!\widehat{\boldsymbol{w}}}\right)^{-1}$ in principle.

The high-dimensionality of context vector $\boldsymbol{x}_{ij}$, however,
does not allow us for materializing each block's full variance-covariance matrix,
whose further diagonalization is neither acceptable due to the poor quantification of uncertainty
when ignoring multi-collinearity among different words in text.

We instead estimate a rank-$k$ ($\ll\!d$) approximation of the inverse variance-covariance matrix as
\begin{equation}
\label{equation:approximate_inverse_covmatrix}
\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{\beta}}^{-1}=%
\boldsymbol{V}_{\boldsymbol{\beta}}%
\boldsymbol{\Lambda}_{\boldsymbol{\beta}}%
\boldsymbol{V}_{\boldsymbol{\beta}}^\top%
+c_{\boldsymbol{\beta}}\boldsymbol{I}_d\mbox{ and }%
\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{\rho}}^{-1}=%
\boldsymbol{V}_{\boldsymbol{\rho}}%
\boldsymbol{\Lambda}_{\boldsymbol{\rho}}%
\boldsymbol{V}_{\boldsymbol{\rho}}^\top%
+c_{\boldsymbol{\rho}}\boldsymbol{I}_d,%
%\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{w}}^{-1}=%
%\begin{pmatrix}%
%\boldsymbol{V}_{0}%
%\boldsymbol{\Lambda}_{0}%
%\boldsymbol{V}_{0}^\top%
%+c_{\boldsymbol{\beta}}\boldsymbol{I}_d
% & %
%\boldsymbol{V}_{1}%
%\boldsymbol{\Lambda}_{1}%
%\boldsymbol{V}_{1}^\top%
% \\
%\boldsymbol{V}_{1}%
%\boldsymbol{\Lambda}_{1}%
%\boldsymbol{V}_{1}^\top%
% & %
%\boldsymbol{V}_{2}%
%\boldsymbol{\Lambda}_{2}%
%\boldsymbol{V}_{2}^\top%
%+c_{\boldsymbol{\rho}}\boldsymbol{I}_d%
%\end{pmatrix}%,
\end{equation}
where 
$\boldsymbol{V}_{\boldsymbol{\beta}},\boldsymbol{V}_{\boldsymbol{\rho}}\!\in\!{\mathbb R}^{d\times k}$,
$\boldsymbol{\Lambda}_{\boldsymbol{\beta}}\!\equiv\!diag(\lambda_{\boldsymbol{\beta} 1},\ldots,\lambda_{\boldsymbol{\beta} k})$, and
$\boldsymbol{\Lambda}_{\boldsymbol{\rho}}\!\equiv\!diag(\lambda_{\boldsymbol{\rho} 1},\ldots,\lambda_{\boldsymbol{\rho} k})$.
By carefully watching Eq.~(\ref{equation:data_hessian}),
one can find that each matrix $\boldsymbol{V}_{(\cdot)}$ is obtained by
thin Singular Value Decomposition (SVD).
Let us assume that users and items are indexed from $i\!=\!1$ to $i\!=\!I$
and from $j\!=\!1$ to $j\!=\!J$, respectively.
For each block, we perform thin SVD of a weighted data matrix as
\begin{equation}
\boldsymbol{X}_{(\cdot)}\!\triangleq\!%
(({h_{11}^{(\cdot)}\vert_{\boldsymbol{w}\!=\!\widehat{\boldsymbol{w}}}})^{1/2}\boldsymbol{x}_{11},%
\ldots,%
({h_{IJ}^{(\cdot)}\vert_{\boldsymbol{w}\!=\!\widehat{\boldsymbol{w}}}})^{1/2}\boldsymbol{x}_{IJ})%
^{\!\!\top}%
\!\!\in{\mathbb R}^{\vert {\mathcal D}\vert\times d}
\mbox{ and }
\boldsymbol{X}_{(\cdot)}\simeq\boldsymbol{U}_{(\cdot)}%
\boldsymbol{\Lambda}_{(\cdot)}^{\frac{1}{2}}%
\boldsymbol{V}_{(\cdot)}^\top\nonumber,
\end{equation}
where $\boldsymbol{U}_{(\cdot)}\!\in\!{\mathbb R}^{\vert{\mathcal D}\vert\times k}$
and matrix $\boldsymbol{\Lambda}_{(\cdot)}^{\frac{1}{2}}$
is given by the $k$-largest singular values of $\boldsymbol{X}_{(\cdot)}$.
%Implementation of the thin SVD in Eq.~(\ref{equation:thin_svd}) should be 
%done on consideration of which of the sample size $\vert{\mathcal D}\vert$ and the input dimensionality $d$ is larger.
Based on the rank-$k$ approximation,
the closed form of each variance-covariance matrix is consequently given as  
\begin{equation}
\label{equation:closed_form_covariance}
\widehat{\boldsymbol{\Sigma}}_{(\cdot)}\!\triangleq\!\left(c_{(\cdot)}\boldsymbol{I}_d\!+\!\boldsymbol{V}_{(\cdot)}\boldsymbol{\Lambda}_{(\cdot)}\boldsymbol{V}_{(\cdot)}^\top\right)^{-1}\!\equiv\!%
c_{(\cdot)}^{-1}\boldsymbol{I}_d\!-\!c_{(\cdot)}^{-1}\boldsymbol{V}_{(\cdot)}(\boldsymbol{I}_d+c\boldsymbol{\Lambda}_{(\cdot)}^{-1})^{-1}\boldsymbol{V}_{(\cdot)}^\top.
\end{equation}

\iffalse
\begin{eqnarray}
\boldsymbol{\Xi}\!&\triangleq&\!%
diag \left(%
\frac{\lambda_{01}}{c_{\boldsymbol{\beta}}+\lambda_{01}},\ldots,%
\frac{\lambda_{0k}}{c_{\boldsymbol{\beta}}+\lambda_{0k}}
\right)\nonumber\\
\boldsymbol{A}^{-1}\!&\triangleq&\!%
c_{\boldsymbol{\beta}}^{-1}\boldsymbol{I}_d\!-\!%
\boldsymbol{V}_0(\boldsymbol{I}_d+c_{\boldsymbol{\beta}}\boldsymbol{\Lambda}_0^{-1})\boldsymbol{V}_0^\top
\nonumber\\
\boldsymbol{D}^{-1}\!&\triangleq&\!%
c_{\boldsymbol{\rho}}^{-1}\boldsymbol{I}_d\!-\!%
\boldsymbol{V}_2(\boldsymbol{I}_d+c_{\boldsymbol{\rho}}\boldsymbol{\Lambda}_2^{-1})\boldsymbol{V}_2^\top
\nonumber\\
\boldsymbol{T}\!&\triangleq&\!%
\boldsymbol{V}_1\boldsymbol{\Lambda}_1\boldsymbol{V}_1^\top\boldsymbol{V}_2(\boldsymbol{I}_d+c_{\boldsymbol{\rho}}\boldsymbol{\Lambda}_2^{-1})^{1/2}\nonumber\\
\boldsymbol{S}^{-1}\!&\triangleq&\!%
\boldsymbol{D}^{-1}-
\boldsymbol{D}^{-1}\boldsymbol{V}_1\boldsymbol{\Lambda}_1\boldsymbol{V}_1^\top%
\left(c_{\boldsymbol{\beta}}\boldsymbol{I}_d+\boldsymbol{V}_0\boldsymbol{\Lambda}_0\boldsymbol{V}_0^\top%
+\boldsymbol{V}_1\boldsymbol{\Lambda}_1\boldsymbol{V}_1^\top\boldsymbol{D}^{-1}\boldsymbol{V}_1\boldsymbol{\Lambda}_1\boldsymbol{V}_1^\top
\right)^{-1}\boldsymbol{V}_1\boldsymbol{\Lambda}_1\boldsymbol{V}_1^\top\boldsymbol{D}^{-1}%
\nonumber\\
\boldsymbol{G}\!&\triangleq&\!%
\boldsymbol{\Lambda}_1(\boldsymbol{V}_1^\top\boldsymbol{V}_0)\boldsymbol{\Xi}%
(\boldsymbol{V}_0^\top\boldsymbol{V}_1)\boldsymbol{\Lambda}_1\in{\mathbb R}^{k\times k}\nonumber\\
\boldsymbol{S}\!&\triangleq&\!%
c_{\boldsymbol{\rho}}\boldsymbol{I}_d+%
\boldsymbol{V}_2\boldsymbol{\Lambda}_2\boldsymbol{V}_2^\top%
+c_{\boldsymbol{\beta}}^{-1}%
\boldsymbol{V}_1(\boldsymbol{G}-\boldsymbol{\Lambda}_1^2)\boldsymbol{V}_1^\top\nonumber\\%
\boldsymbol{A}\!&\triangleq&\!%
c_{\boldsymbol{\beta}}^{-1}\boldsymbol{I}_d-%
c_{\boldsymbol{\beta}}^{-1}\boldsymbol{V}_0%
\boldsymbol{\Xi}
\boldsymbol{V}_0^\top\nonumber\\
\boldsymbol{B}\!&\triangleq&\!%
\boldsymbol{\Xi}
\left(\boldsymbol{V}_0^\top\boldsymbol{V}_1\right)\boldsymbol{\Lambda}\in{\mathbb R}^{k\times k}\nonumber\\
\boldsymbol{C}\!&=&\!%
c_{\boldsymbol{\beta}}^{-1}\boldsymbol{V}_1\boldsymbol{\Lambda}_1\boldsymbol{V}_1^\top-
c_{\boldsymbol{\beta}}^{-1}\boldsymbol{V}_0\boldsymbol{B}\boldsymbol{V}_1^\top%
\end{eqnarray}
\begin{equation}
\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{w}}=%
\begin{pmatrix}%
c_{\boldsymbol{\beta}}^{-1}\boldsymbol{V}_1\boldsymbol{\Lambda}_1\boldsymbol{V}_1^\top-
c_{\boldsymbol{\beta}}^{-1}\boldsymbol{V}_0\boldsymbol{\Xi}\boldsymbol{V}_1^\top%
 & %
\boldsymbol{V}_{(\boldsymbol{\beta},\boldsymbol{\rho})}%
\boldsymbol{V}_{(\boldsymbol{\beta},\boldsymbol{\rho})}^\top%
 \\
\boldsymbol{V}_{(\boldsymbol{\beta},\boldsymbol{\rho})}%
\boldsymbol{V}_{(\boldsymbol{\beta},\boldsymbol{\rho})}^\top%
 & %
\boldsymbol{V}_{(\boldsymbol{\rho},\boldsymbol{\rho})}%
\boldsymbol{V}_{(\boldsymbol{\rho},\boldsymbol{\rho})}^\top%
+c_{\boldsymbol{\rho}}\boldsymbol{I}_d%
\end{pmatrix}%
.
\end{equation}
\fi

\subsection{Hyperparameter optimization with approximate marginal likelihood}
\label{subsec:hyperparameter_section}

Laplace's method also provides a closed-form approximation of the marginal likelihood and 
enables to optimize the hyperparameter vector $\boldsymbol{c}$ in a Bayesian manner.
Our loss is quadratically approximated around the MAP estimate $\widehat{\boldsymbol{w}}$ as
${\mathcal L}_{{\mathcal D},\boldsymbol{c}}(\boldsymbol{w})\!\simeq\!%
{\mathcal L}_{{\mathcal D},\boldsymbol{c}}(\widehat{\boldsymbol{w}})\!+\!%
\frac{1}{2}(\boldsymbol{w}\!-\!\widehat{\boldsymbol{w}})^\top%
\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{w}}^{-1}(\boldsymbol{w}\!-\!\widehat{\boldsymbol{w}})$.
The marginal negative log-likelihood
$\widetilde{{\mathcal L}}_{{\mathcal D},\boldsymbol{c}}\!\triangleq\!%
-\log\int d\boldsymbol{w}\prod_{(i,j)\in{\mathcal D}}p(v_{ij}\vert n_{ij}, \boldsymbol{x}_{ij}, \boldsymbol{w})p(\boldsymbol{w})$ is hence approximated as
\begin{eqnarray}
\exp(-\widetilde{{\mathcal L}}_{{\mathcal D},\boldsymbol{c}})\!&\simeq&\!%
p(\widehat{\boldsymbol{w}})\!\!\prod_{(i,j)\in{\mathcal D}}\!\!\!p(v_{ij}\vert n_{ij}, \boldsymbol{x}_{ij}, \widehat{\boldsymbol{w}})%
\int\exp\left(-\frac{1}{2}(\boldsymbol{w}\!-\!\widehat{\boldsymbol{w}})^\top%
\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{w}}^{-1}(\boldsymbol{w}\!-\!\widehat{\boldsymbol{w}})\right)d\boldsymbol{w}%
\nonumber\\
\!&\equiv&\!
\exp\left(-{\mathcal L}_{{\mathcal D},\boldsymbol{c}}(\widehat{\boldsymbol{w}})\right)
\vert\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{\beta}}\vert^{1/2}
\vert\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{\rho}}\vert^{1/2}
c_{\boldsymbol{\beta}}^{d/2}c_{\boldsymbol{\rho}}^{d/2}\nonumber.
\end{eqnarray}
%which is a product of the joint likelihood at the MAP estimate
%and the square root of determinant-ratio between posterior and prior variance-covariance matrices.
Because each determinant is the product of diagonal elements of the matrix $\boldsymbol{\Lambda}_{(\cdot)}$ in (\ref{equation:approximate_inverse_covmatrix}),
the approximate marginal negative log-likelihood to select the hyperparameter $\boldsymbol{c}$ is 
\begin{equation}
\label{equation:marginal_loss}
\widetilde{{\mathcal L}}_{{\mathcal D},\boldsymbol{c}}
\simeq
{\mathcal L}_{{\mathcal D},\boldsymbol{c}}(\widehat{\boldsymbol{w}})
+\frac{1}{2}\sum_{\ell=1}^k
\log\left(1+\frac{\lambda_{\boldsymbol{\beta}\ell}}{c_{\boldsymbol{\beta}}}\right)%
\!+\frac{1}{2}\sum_{\ell=1}^k\log\left(1+\frac{\lambda_{\boldsymbol{\rho}\ell}}{c_{\boldsymbol{\rho}}}\right)\nonumber.
\end{equation}
Based on the $L_2$-regularization terms in the loss function ${\mathcal L}_{{\mathcal D},\boldsymbol{c}}$ and Eq.~(\ref{equation:marginal_loss}),
we maximize the marginal likelihood by iterating between 
the update of $(\widehat{\boldsymbol{w}}, \widehat{\boldsymbol{\Sigma}}_{\boldsymbol{w}})$ and that of $\boldsymbol{c}$ as
\begin{equation}
c_{\boldsymbol{\beta}}\!=\!\min_{c}\left[c\Vert\widehat{\boldsymbol{\beta}}\Vert^2\!+\!\sum_{\ell=1}^k\log\left(1\!+\!\frac{\lambda_{\boldsymbol{\beta}\ell}}{c}\right)\right]\mbox{ and }%
c_{\boldsymbol{\rho}}\!=\!\min_{c}\left[c\Vert\widehat{\boldsymbol{\rho}}\Vert^2\!+\!\sum_{\ell=1}^k\log\left(1\!+\!\frac{\lambda_{\boldsymbol{\rho}\ell}}{c}\right)\right]\nonumber.
\end{equation}

\section{Uncertainty-aware recommendation}
\label{sec:recommendation}

In order to rank test items, this section introduces several scoring measures 
which are all derived from the approximate posterior or predictive distribution.
While the exact predictive distribution is intractable
due to the non-conjugacy between Beta and Gaussian distributions,
its Monte Carlo approximation is easily obtained as we show 
in Section \ref{subsec:monte_carlo_predictive}. The Monte Carlo approach
eases computation of the scoring measures, whose varieties and characteristics are described in
\ref{subsec:varieties_scores}. 
Dependence of recommendation results on the choice of scoring methods is evaluated in Section \ref{sec:experiments}. 

\subsection{Monte-Carlo predictive distribution}
\label{subsec:monte_carlo_predictive}

Let $\theta_{*}\!\in\![0,1]$  be probability of success in test context associated with
vector $\boldsymbol{x}_{*}\!\in\!{\mathbb R}^d$.
Integrating over the approximate posterior of $\boldsymbol{w}$,
we obtain the predictive distribution of $\theta_{*}$ conditional on $\boldsymbol{x}_{*}$.
The approximate predictive distribution is defined as
$q(\theta_{*}\vert\boldsymbol{x}_{*}, {\mathcal D},\boldsymbol{c})\!\triangleq\!%
{\mathbb E}_{q(\boldsymbol{w}\vert{\mathcal D},\boldsymbol{c})}\left[p(\theta_{*}\vert\boldsymbol{x}_{*}, \boldsymbol{w})\right]\!=\!%
{\mathbb E}_{{\mathcal N}\left(\boldsymbol{\beta}; \widehat{\boldsymbol{\beta}},\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{\beta}}\right)%
{\mathcal N}\left(\boldsymbol{\rho}; \widehat{\boldsymbol{\rho}},\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{\rho}}\right)}
\!\left[Be\!\left(\theta_{*}; %
\exp((\boldsymbol{\beta}\!+\!\boldsymbol{\rho})\!^\top\!\boldsymbol{x}_{*}), %
\exp(\boldsymbol{\rho}\!^\top\!\boldsymbol{x}_{*})\right)\right]$, whose inherent integral is analytically intractable
while is numerically well-approximated by a bidimensional Monte Carlo integration.
Since each block of the variance-covariance matrix
has the common form (\ref{equation:closed_form_covariance}),
variances of $\zeta_{*}\!\triangleq\!\boldsymbol{\beta}^\top\boldsymbol{x}_{*}$ and $\eta_{*}\!\triangleq\!\boldsymbol{\rho}^\top\boldsymbol{x}_{*}$ are
both given as
$\boldsymbol{x}_{*}^\top\widehat{\boldsymbol{\Sigma}}_{(\cdot)}\boldsymbol{x}_{*}\!\equiv\!c_{(\cdot)}^{-1}\boldsymbol{x}_{*}^\top\boldsymbol{x}_{*}\!-\!c_{(\cdot)}^{-1}\boldsymbol{x}_{*}^\top\boldsymbol{V}_{(\cdot)}(\boldsymbol{I}_d+c_{(\cdot)}\boldsymbol{\Lambda}_{(\cdot)}^{-1})^{-1}\boldsymbol{V}_{(\cdot)}^\top\boldsymbol{x}_{*}$. Therefore, an $m$-sample Monte Carlo approximation of the predictive distribution is given as
\begin{eqnarray}
\sigma_{(\cdot)}^2(\boldsymbol{x}_{*})\!\!\!\!&\triangleq&\!\!\!\!%
c_{(\cdot)}^{-1}\!\left(\Vert\boldsymbol{x}_{*}\Vert^2%
\!-\!\sum_{\ell=1}^k \frac{\lambda_{(\cdot)\ell}}{\lambda_{(\cdot)\ell}+c_{(\cdot)}}y_{*(\cdot)\ell}^2\right)%
\mbox{where }(y_{(\cdot)*1},\ldots,y_{(\cdot)*k})^\top\!\triangleq\!\boldsymbol{V}_{(\cdot)}^\top\boldsymbol{x}_{*}
\nonumber,\\
%\sigma_{\boldsymbol{\rho}}^2(\boldsymbol{x}_{*})\!&\triangleq&\!%
%c_{\boldsymbol{\rho}}^{-1}\!\left(\Vert\boldsymbol{x}_{*}\Vert^2%
%\!-\!\sum_{\ell=1}^k \frac{\lambda_{\boldsymbol{\rho}\ell}}{\lambda_{\boldsymbol{\rho}\ell}+c_{\boldsymbol{\rho}}}z_{*\ell}%^2\right)
%\mbox{ where }(z_{*1},\ldots,z_{*k})^\top\!\triangleq\!\boldsymbol{V}_{\boldsymbol{\rho}}^\top\boldsymbol{x}_{*}
%\nonumber,\\
\label{equation:monte_carlo_scores}
\begin{pmatrix}%
\zeta_{*}^{(1)}\\
\eta_{*}^{(1)}%
\end{pmatrix}%
\!,\ldots\!,%
\begin{pmatrix}%
\zeta_{*}^{(m)}\\
\eta_{*}^{(m)}%
\end{pmatrix}%
\!\!\!\!&\sim&\!\!\!\!%
{\mathcal N}\!\left(%
\!\!\begin{pmatrix}%
\widehat{\boldsymbol{\beta}}^\top\boldsymbol{x}_{*}\\
\widehat{\boldsymbol{\rho}}^\top\boldsymbol{x}_{*}
\end{pmatrix}\!, %
\begin{pmatrix}%
\sigma_{\boldsymbol{\beta}}^2(\boldsymbol{x}_{*}) & 0\\
0 & \sigma_{\boldsymbol{\rho}}^2(\boldsymbol{x}_{*})
\end{pmatrix}\right)\mbox{, and }\\
\label{equation:approximate_predictive_distribution}
q(\theta_{*}\vert\boldsymbol{x}_{*}, {\mathcal D},\boldsymbol{c})\!\!\!\!&\simeq&\!\!\!\!%
\frac{1}{m}\sum_{t=1}^m
Be\!\left(\theta_{*}; %
\exp(\zeta_{*}^{(t)}+\eta_{*}^{(t)}), %
\exp(\zeta_{*}^{(t)})\right).%
\end{eqnarray}
The total complexity in (\ref{equation:monte_carlo_scores}) is ${\mathcal O}(\max\lbrace dk, m\rbrace)$ while it is
much lower for sparse $\boldsymbol{x}_{*}$.
Thanks to the low-dimensionality of the integral,
Quasi-Monte Carlo method makes the integration more accurate.

\subsection{Varieties of recommendation scores}
\label{subsec:varieties_scores}

By using Eqs.~(\ref{equation:monte_carlo_scores}) and (\ref{equation:approximate_predictive_distribution}), we can
derive several types of scores used in the final recommendation.
With regarding probability of success as our target variable, 
we introduce a variety of measures based on Upper Confidence Bounds (UCBs) or quantiles. While
many of our measures stem from existing approaches to handle exploration-exploitation trade-off in multi-armed bandit or Bayesian optimization,
our recommendation experiments in Section \ref{sec:experiments} are not repetitious but one-time.

\paragraph{Upper Confidence Bound of Expectation (UCBE)}

UCB in multi-armed bandit \cite{yue2011linear} or Bayesian optimization 
\cite{contal2014gaussian} is usually defined as
an optimistic estimate of the expected reward under uncertainty. In our model,
the expected reward under no uncertainty of $\boldsymbol{w}$
is ${\mathbb E}_{p(\theta_{*}\vert\boldsymbol{x}_{*}, \boldsymbol{w})}[\theta_{*}]\!=\!1/(1+\exp(-\boldsymbol{\beta}^\top\!\boldsymbol{x}_{*}))$, 
which does not depend on $\boldsymbol{\rho}$.
By using the Monte Carlo samples $\zeta_{*}^{(1)},\ldots,\zeta_{*}^{(m)}$ in (\ref{equation:monte_carlo_scores}),
we obtain the mean and standard deviation about the probability of success as
\begin{equation}
\mu_{*}\!\simeq\!\frac{1}{m}\sum_{t=1}^m \frac{1}{1+\exp\left(-\zeta_{*}^{(t)}\right)}\mbox{ and }%
s_{*}^2\!\simeq\!\frac{1}{m\!-\!1}\sum_{t=1}^m \left(\frac{1}{(1+\exp\left(-\zeta_{*}^{(t)}\right)}\!-\!\mu_{*}\right)^2\nonumber,
\end{equation}
respectively. Then the $\nu$-UCBE is finally computed as
\begin{equation}
\label{equation:our_ucb}
{UCBE}_{q(\theta_{*}\vert\boldsymbol{x}_{*}, {\mathcal D},\boldsymbol{c})}^{(\nu)}[\theta_{*}]\!\simeq\!%
\mu_{*}+\Phi^{-1}(\nu)s_{*}\nonumber,
\end{equation}
where $\Phi^{-1}(\cdot)$ is the inverse cumulative distribution function of the standard normal distribution
and $\nu\!>\!0.5$ is a hyperparameter to determine the exploration-exploitation trade-off.
While Eq.~(\ref{equation:our_ucb}) does not explicitly depend on $\boldsymbol{\rho}$, 
the additional complexity by $\boldsymbol{\rho}$ affects the estimate of $\boldsymbol{\beta}$.
The explicit handling of over-dispersion more robustifies UCB than the standard logistic regression.

\paragraph{Upper Confidence Quantile of Expectation (UCQE)}

Another optimistic statistic easily obtained from the posterior is upper quantile of 
the expected reward. Specifically, we compute $100\nu$-percentile 
of ${\mathbb E}_{p(\theta_{*}\vert\boldsymbol{x}_{*}, \boldsymbol{w})}\left[\theta_{*}\right]$ over the posterior 
$q(\theta_{*}\vert\boldsymbol{x}_{*}, {\mathcal D},\boldsymbol{c})$. Because the sigmoid function is a monotonic transformation, 
quantile of the sigmoid is the sigmoid of the quantile, whose computation does not require Monte Carlo samples.
Hence our $\nu$-UCQE measure is given as
\begin{equation}
{UCQE}_{q(\theta_{*}\vert\boldsymbol{x}_{*}, {\mathcal D},\boldsymbol{c})}^{(\nu)}[\theta_{*}]\!=\!%
\frac{1}{1+\exp(-\zeta_{*\nu})}\mbox{ where }
\zeta_{*\nu}=\widehat{\boldsymbol{\beta}}^\top\boldsymbol{x}_{*}+\Phi^{-1}(\nu)\sigma_{\boldsymbol{\beta}}(\boldsymbol{x}_{*})\nonumber.
\end{equation}
Also for UCQE, modeling the over-dispersion introduced by $\boldsymbol{\rho}$
supplies robustification.

\paragraph{Expected Upper Quantile (EUQ)}

Upper quantile of the reward is a valuable statistic when we are interested in outliers
whose probability $\theta_{*}$
is not perfectly predictable from context vector $\boldsymbol{x}_{*}$.
In our case, the quantile of reward under no uncertainty of $\boldsymbol{w}$
is given by the inverse cumulative distribution function of beta distribution 
$F^{-1}(\cdot; \exp((\boldsymbol{\beta}+\boldsymbol{\rho})^\top\boldsymbol{x}_{*}), \exp(\boldsymbol{\rho}^\top\boldsymbol{x}_{*}))$. 
On the Monte Carlo samples 
$(\zeta_{*}^{(1)},\eta_{*}^{(1)}),\ldots,(\zeta_{*}^{(m)},\eta_{*}^{(m)})$ in (\ref{equation:monte_carlo_scores}),
the $\nu$-EUQ measure is empirically computed as
\begin{equation}
{EUQ}_{q(\theta_{*}\vert\boldsymbol{x}_{*}, {\mathcal D},\boldsymbol{c})}^{(\nu)}[\theta_{*}]\!=\!%
\frac{1}{m}\sum_{t=1}^mF^{-1}\!\left(\nu; %
\exp(\zeta_{*}^{(t)}+\eta_{*}^{(t)}), %
\exp(\zeta_{*}^{(t)})\right)\nonumber.%
\end{equation}

\paragraph{Upper Confidence Quantile of Upper Quantile (UCQUQ)}

A further optimistic statistic
is obtained through replacing the expectation in EUQ by upper quantile.
The $(\nu_1,\nu_2)$-UCQUQ measure is the $100\nu_2$-percentile of
the empirical posterior distribution about the $100\nu_1$-percentile of reward, as
\begin{eqnarray}
{UCQUQ}_{q(\theta_{*}\vert\boldsymbol{x}_{*}, {\mathcal D},\boldsymbol{c})}^{(\nu_1,\nu_2)}[\theta_{*}]\!&=&\!G_{m,\nu_1}^{-1}(\nu_2)\mbox{ and }\nonumber\\
G_{m,\nu_1}(\theta)\!&\triangleq&\!%
\frac{1}{m}\sum_{t=1}^m U\left(\theta-F^{-1}\!\left(\nu_1; %
\exp(\zeta_{*}^{(t)}+\eta_{*}^{(t)}), %
\exp(\zeta_{*}^{(t)})\right)\right)\nonumber,%
\end{eqnarray}
where $U(\cdot)$ is the unit step function
and $G_{m,\nu_1}^{-1}(\cdot)$ is the inverse function of $G_{m,\nu_1}(\cdot)$.

\paragraph{Upper Quantile of Predictive distribution (UQP)}

Predictive distribution can also supply an optimistic statistic
that has a similar principle to EUQ. Here the quantile is taken after the marginalization
over the posterior.
Upper quantile of the empirical predictive distribution is given as
\begin{equation}
\label{equation:upper_quantile_predictive}
{UQP}_{q(\theta_{*}\vert\boldsymbol{x}_{*}, {\mathcal D},\boldsymbol{c})}^{(\nu)}[\theta_{*}]\!=\!\theta_{\nu}\mbox{ such that }%
\nu=\frac{1}{m}\sum_{t=1}^m%
F\!\left(\theta_{\nu}; %
\exp(\zeta_{*}^{(t)}+\eta_{*}^{(t)}), %
\exp(\zeta_{*}^{(t)})\right).%
\end{equation}
While there is no closed-form formula of $\theta_{\nu}$
in (\ref{equation:upper_quantile_predictive}),
we can numerically compute the quantile
by applying bi-section method with the cumulative distribution function of Beta distribution
$F(\cdot; \cdot, \cdot)$.

\iffalse
\subsection{On the tail of Beta distribution}

Beta-binomial distribution has a finite variance but a heavier tail than binomial distribution,
how posterior shape is skewed and tailed.
One could be interested how the mean and dispersion parameters are correlated.
Figure \ref{} shows a simple example when the input dimensionality is $d\!=\!1$.

For more rigorous analysis,
point estimate that only uses the posterior mode is expected to provide severely wrong estimate
for the marginal distribution of $\theta_{ij}$.

In other words, 
posterior has a {\em plateau} around the mode.
When we explicitly consider the fundamental unpredictability of $\theta_{ij}$ only by $\boldsymbol{x}_{ij}$.
Ignoring this plateau hugely misleads the 

\fi


\section{Related work}
\label{sec:related_work}

% Beta-binomial distribution, robust confidence interval
Beta-binomial-logit models have been used for robust classification
whereas their over-dispersion has been supplied by not a regression formula but by a scholar hyperparameter (e.g., \cite{tak2017data}).
Contextual-risk models have been used in regression tasks
such as Gaussian Process (GP) regression with input-dependent variances \cite{goldberg1998regression},
while have been uncommon in classification tasks.
Our derivation of the custom Laplace approximation and marginal likelihood
are based on the techniques in GP classification \cite{rasmussen2005gaussian},
for which Expectation Propagation (EP; \cite{minka2001expectation}) is also applicable whereas we avoided too complicated formulas of EP.
Nonparametric conditional density estimation (e.g., \cite{shahbaba2009nonlinear,sugiyama2010conditional})
naturally introduces input-dependent noises while simpler forms are preferable in our task.
$T$-logistic regression \cite{ding2010logistic} is another robust classifier and
Bayesian estimation of robust classifiers
produces robust credible intervals (e.g., \cite{pericchi1991robust}),
while their risks do not depend on inputs.
Overall, to the best of our knowledge, we provide the most parsimonious Bayesian classifier
that suits recommendation of outlying items based on input-dependent variabilities and uncertainties.

% Human behavior on exploration
Humans exhibit systematically predictable behaviors in the face of uncertainty.
One reliable observation is that desire to avoid monetary loss is a strong incentive for exploration
(e.g., win-stay lose-shift algorithm \cite{robbins1952some}, prospect theory \cite{kahneman1979prospect},
loss aversion \cite{usher2004loss}, and regulatory focus theory
\cite{higgins2000making,avnet2006how}) and 
compensating money works as an incentive \cite{frazier2014incentivizing}.
Users do not lose money, however, in online news service when they read a narrow range of articles.
Even in mental level, it is uncertain whether reading merely one-sided opinions lets users feel pains.
Another observation is that diminishing return for the same type of stimulus
naturally leads exploration (e.g., \cite{smith1999diagnosing}, variety-seeking behavior in marketing \cite{kumar2006estimation,ho2014towards}) 
to maintain the Optimum Stimulation Level (OSL) \cite{raju1980optimum}.
Diminishing return has already been exploited in diversified recommendation (e.g., linear submodular bandits \cite{yue2011linear}),
while we already discussed the insufficient power of diversification for surely unpopular articles.

\section{Experimental evaluations}
\label{sec:experiments}

We experimentally evaluate the varieties of our scoring methods.
In Section \ref{subsec:performance_measure},
we define a performance indicator of early detectability for special articles 
that are chosen in the test period despite their dissimilarity to the positive training samples.
%we evaluate the test-set likelihood for synthetic datasets in Section \ref{subsec:synthetic_datasets}.
%The main purpose for using synthetic dataset
%is to assess whether the additional complexity to handle the over-dispersion
%is spoiled by datasets without over-dispersion.
For real-world news-reading datasets and reference models introduced in Section \ref{subsec:experimental_setting},
we compare the performances among the proposed and reference models in Section \ref{subsec:performance_comparisons}.
The proposed models outperform the reference models in terms of the 
serendipity-oriented indicator, while achieving competitive levels of test-set likelihood.

\subsection{Serendipity-oriented performance measure}
\label{subsec:performance_measure}

Our main indicator is a variation of Area Under Curve
with prioritization of articles whose popularity is hard to be early detected.
Let $\boldsymbol{S}_i
\!\triangleq\!(S_{ijj'}\!\equiv\!\boldsymbol{x}_{ij}^\top\boldsymbol{x}_{ij'})\!\in\!{\mathbb R}^{N[i]\times N[i]}$ be
a  similarity matrix among all of the items that user $i$ chose in the training period, where
$N[i]\!\triangleq\!\vert\lbrace j; v_{ij}>0\wedge (i,j)\!\in\!{\mathcal D}\rbrace\vert$.
For test vector of context $\boldsymbol{x}_{i*}$,
we define an augmented $(N[i]\!+\!1)$-by-$(N[i]\!+\!1)$ matrix
$\boldsymbol{S}_{i*}$, via taking the inner product between the test vector
$\boldsymbol{x}_{i*}$ and every training vector assigned with a positive label.
By borrowing the common formula between the partition function in DPPs \cite{kulesza2013determinantal,kulesza2011learning}
and exponentiation of mutual information of GP \cite{contal2014gaussian},
we define reward variable $r_{i*}\!\geq\!1$ for this test context as
\begin{equation}
\label{diversity_gain}
r_{i*}\!\triangleq\!\dfrac{\det (\boldsymbol{I}_{N[i]+1}\!+\!\boldsymbol{S}_{i*})}{\det (\boldsymbol{I}_{N[i]}\!+\!\boldsymbol{S}_{i})}%
\!\equiv\!1\!+\!\boldsymbol{x}_{i*}^\top\boldsymbol{x}_{i*}%
\!-\!\left((\boldsymbol{X}_i\boldsymbol{x}_{i*})^\top%
(\boldsymbol{I}_{N[i]}+\boldsymbol{X}_{i}\boldsymbol{X}_i^\top)^{-1}%
(\boldsymbol{X}_i\boldsymbol{x}_{i*})%
\right)\!,%
\end{equation}
where $\boldsymbol{X}_i\!\in\!{\mathbb R}^{N[i]\times d}$ is user $i$-specific design matrix 
that lines up all of the positively-labeled context vectors.
Eq.~(\ref{diversity_gain}) represents a monotonically non-decreasing gain of diversity among the articles chosen by user $i$, 
when the test article is added into the existing selection. 
By setting each horizontal length $r_{i*}$ and unifying the test samples by all of the users,
we draw one Receiver Operator Characteristic curve.
We name the resulting performance measure the Serendipity-oriented Area Under Curve (SAUC).
While other reward variables are considerable,
we regard Eq.~(\ref{diversity_gain})
as a good starting point for further studies about serendipity,
because of its direction connection with diversity and entropy,


\subsection{Real-world news datasets and reference models}
\label{subsec:experimental_setting}

\begin{table}[t]
\centering
\caption{Basic characteristics of the four news datasets in 2016.
Each of the {\tt Brexit}, {\tt USElect}, and {\tt FBIMail} datasets covers one-week 
records of individual-level impressions and pageviews for selected users,
where a big political news occurred around the boundary between
the training and test periods. The {\tt Normal} dataset is introduced for 
evaluation when there is no big political news.}
\label{table:characteristics_datasets}
\centering
\begin{tabular}{r|r|r|r|r|r|r}
Dataset Name&Training Period&Test Period&\#users&\#items&$\sum_{(i,j)\in{\mathcal D}}v_{ij}$&$\sum_{(i,j)\in{\mathcal D}}n_{ij}$\\\hline
{\tt Brexit}&Jun 20-23&Jun 24-26&9,908&9,026&221,781&2,388,617\\\hline
{\tt USElect}&Nov 5-8&Nov 9-11&9,841&2,419&266,331&2,792,707\\\hline
{\tt FBIMail}&Aug 29-Sep 1&Sep 2-4&9,871&2,123&184,492&2,060,274\\\hline
{\tt Normal}&Apr 1-4&Apr 5-7&9,896&6,618&186,557&1,995,538\\\hline
\end{tabular}
\end{table}

Our news datasets consist of individual-user-level 
page view and impression logs of a mobile-app news service in the United States.
Table \ref{table:characteristics_datasets} shows the characteristics of our 4 datasets.

The vector of context $\boldsymbol{x}_{ij}$ is given by text content of article
and clustering of users.
For each user $i$, we have a binary vector 
to indicate the items that user $i$ chose before the training period.
With standardizing the $L_2$-norm of every binary vector,
we performed clustering of users by spherical $64$-means algorithm.
We assume that cluster label represents each user's preference
that does not change during every one-week period.
Let $\boldsymbol{x}_j\!\in\!{\mathbb R}^{d_0}$ be unit-norm TF-IDF vector 
of article $j$. 
For multi-task learning \cite{chapelle2014simple}, we set
$\boldsymbol{x}_{ij}\!=\!%
(\boldsymbol{x}_j^\top,\boldsymbol{0}_{d_0}^\top,\ldots,\boldsymbol{0}_{d_0}^\top,\boldsymbol{x}_j^\top,%
\boldsymbol{0}_{d_0}^\top,\ldots,\boldsymbol{0}_{d_0}^\top)^\top\!\in\!{\mathbb R}^{(64\!+\!1)d_0}$ where the positions of non-zero elements are determined based on the cluster label of user $i$.
This feature design is used for all of the models.
Design of more sophisticated feature vector is out of this paper's scope.

Because our model's structural advantage comes from the context-dependence of variabilities and uncertainties,
the reference methods are {\tt M-Log}: ordinary $L_2$-regularized MAP logistic regression,
{\tt M-BBL}: MAP estimate of a Beta-binomial-logit model with input-independent over-dispersion \cite{tak2017data},
{\tt L-Log} and  {\tt L-BBL}: Bayesian extensions of {\tt M-Log} and  {\tt M-BBL} by Laplace's method, respectively.
{\tt L-Log} is a linear GP classifier and see \cite{rasmussen2005gaussian} for the concrete formulas.
{\tt M-BBL} and {\tt L-BBL} are the estimates when we replace $\boldsymbol{\rho}\!^\top\!\boldsymbol{x}_{ij}$ in (\ref{equation:beta_model}) by scholar parameter $\rho_0$.
The proposed model and its MAP counterpart are named {\tt L-Prop} and {\tt M-Prop}, respectively.
All of the models are fitted with empirical-Bayes method,
where every $L_2$-regularization hyperparameter $c_{(\cdot)}$
is optimized by gradient descent with initialization such that 
$c_{(\cdot)}\!:=\!\sum_{(i,j)\in{\mathcal D}}n_{ij}$.


\subsection{Performance comparisons}
\label{subsec:performance_comparisons}

Figure \ref{figure:performance} exhibits several performance comparisons
by using SAUC and test-set log-likelihood per impression. AUC is also shown
for comparison. The proposed model achieved the highest likelihood for all of the 4 datasets
and highest SAUC for 2 datasets.
While the differences among models and scoring measures are currently marginal,
the advantage of the proposed model is significant
because of the negligibly small standard deviation among the 5 bootstrap folds.

\begin{figure*}[t]
  \centering
\begin{minipage}{0.495\textwidth}%
  \centering
  \includegraphics[width=0.95\textwidth]{figure/SAUC.eps}\\
\vspace*{-14pt}%
SAUC by MAP estimate or 0.95-UCQE
\vspace*{14pt}%
%  \\
%  Average log-likelihood
\end{minipage}%
\begin{minipage}{0.495\textwidth}%
  \centering
  \includegraphics[width=0.95\textwidth]{figure/AUC.eps}\\
\vspace*{-14pt}%
AUC by MAP estimate or 0.95-UCQE
\vspace*{14pt}%
%  \\
%  Accuracy
\end{minipage}\\
\centering
\begin{minipage}{0.495\textwidth}%
  \centering
  \includegraphics[width=0.95\textwidth]{figure/measures.eps}\\
\vspace*{-14pt}%
SAUC by various scoring\\
within the proposed model
\vspace*{14pt}%
%  \\
%  Average log-likelihood
\end{minipage}%
\begin{minipage}{0.495\textwidth}%
  \centering
  \includegraphics[width=0.95\textwidth]{figure/AvgLL.eps}\\
\vspace*{-14pt}%
Average log-likelihood\\per impression
\vspace*{14pt}%
%  \\
%  Accuracy
\end{minipage}% 
\caption{Performances measures when $\nu_1\!\equiv\!\nu_2\!\equiv\!\nu\!\equiv\!0.95$.
Every score is 5-fold bootstrap mean for test dataset,
where the deviation among the 5 folds was found to be negligibly small.
For fair comparison, the test likelihood for the Beta-binomial models
is evaluated independently for each impression, while the training likelihood is not independent.
The proposed model {\tt L-Prop} achieved the highest likelihood for all of the 4 datasets
while its SAUC and AUC are highest only for 2 datasets {\tt USElect} and {\tt Normal}.
Among the scoring methods, expectation-oriented
UCBE and UCQE measures or predictive-distribution-based UQP are shown to be promising.
}
\label{figure:performance}
\end{figure*}

\section{Conclusion}
\label{sec:conclusion}
This paper introduced a new contextual Bayesian binary classifier
whose risks and/or uncertainties on probabilities of success have been underestimated by existing models.
The proposed model consists of a conditional Beta distribution on input vector of context,
and its estimation is done with its low-rank Laplace approximation using thin SVD.
The closed forms for our approximate posterior yield several uncertainty-aware scoring measures,
and parts of them were experimentally successful in early detecting the exceptional articles
that are clicked despite the dissimilarity to the positive examples in training period.

In the future, we will apply our philosophy on context-dependent risks and uncertainties
for structured non-linear models. We will also qualitatively investigate
article examples associated with high uncertainties.
We hope our work to become a catalyst for further studies in computational journalism.

\newpage
\small{
\bibliographystyle{unsrt}
\bibliography{anti_filter_bubble}
}
\end{document}

anti_filter_bubble.bbl
\begin{thebibliography}{10}

\bibitem{chu2009personalized}
W.~Chu and S.-T. Park.
\newblock Personalized recommendation on dynamic content using predictive
  bilinear models.
\newblock In {\em Proceedings of the 18th International Conference on World
  Wide Web (WWW 2009)}, pages 691--700, New York, NY, USA, 2009. ACM.

\bibitem{sharma2015feature}
M.~Sharma, J.~Zhou, J.~Hu, and G.~Karypis.
\newblock Feature-based factorized bilinear similarity model for cold-start
  top-n item recommendation.
\newblock In {\em Proceedings of the 15th SIAM International Conference on Data
  Mining (SDM 2015)}, pages 190--198, 2015.

\bibitem{wu2010topic}
Y.~Wu, Y.~Ding, X.~Wang, and J.~Xu.
\newblock Topic based automatic news recommendation using topic model and
  affinity propagation.
\newblock In {\em Proceedings of the 27th international conference on Machine
  Learning and Cybernetics (ICMLC 2010)}, pages 1299--1304, 2010.

\bibitem{das2007google}
A.~Das, M.~Datar, A.~Garg, and S.~Rajaram.
\newblock {G}oogle news personalization: Scalable online collaborative
  filtering.
\newblock In {\em Proceedings of the 16th international conference on World
  wide web (WWW 2007)}, pages 271--280, New York, NY, USA, 2010. ACM.

\bibitem{agarwal2010matrix}
D.~Agarwal and B.-C. Chen.
\newblock f{L}{D}{A}: Matrix factorization through latent {D}irichlet
  allocation.
\newblock In {\em Proceedings of the 3rd ACM International Conference on Web
  Search and Data Mining (WSDM 2010)}, pages 91--100, New York, NY, USA, 2010.
  ACM.

\bibitem{gunasekar2015consistent}
S.~Gunasekar, M.~Yamada, D.~Yin, and Y.~Chang.
\newblock Consistent collective matrix completion under joint low rank
  structure.
\newblock In {\em Proceedings of the 18th international conference on
  Artificial Intelligence and Statistics (AISTATS 2015)}, volume~38, 2015.

\bibitem{dali2010learning}
L.~Dali, B.~Fortuna, and J.~Rupnik.
\newblock Learning to rank for personalized news article retrieval.
\newblock In {\em Proceedings of Workshop on Applications of Pattern Analysis
  (WAPA)}, pages 152--159, 2010.

\bibitem{claypool1999combining}
M.~Claypool, A.~Gokhale, T.~Miranda, P.~Murnikov, D.~Netes, and M.~Sartin.
\newblock Combining content-based and collaborative filters in an online
  newspaper.
\newblock In {\em Proceedings of ACM SIGIR Workshop on Recommender Systems},
  1999.

\bibitem{zheng2013penetrate}
L.~Zheng, L.~Li, W.~Hong, and T.~Li.
\newblock {P}{E}{N}{E}{T}{R}{A}{T}{E}: Personalized news recommendation using
  ensemble hierarchical clustering.
\newblock {\em Expert Systems with Applications}, 40:2127--2136, 2013.

\bibitem{pariser2011filter}
Eli Pariser.
\newblock {\em The {F}ilter {B}ubble: How the New Personalized Web Is Changing
  What We Read and How We Think}.
\newblock The Penguin Press, 2011.

\bibitem{bakshy2015exposure}
E.~Bakshy, S.~Messing, and L.~Adamic.
\newblock Exposure to ideologically diverse news and opinion on {F}acebook.
\newblock {\em Science 5}, 348(6239):1130--1132, 2015.

\bibitem{yue2011linear}
Y.~Yue and C.~Guestrin.
\newblock Linear submodular bandits and their application to diversified
  retrieval.
\newblock In {\em Advances in Neural Information Processing Systems 24}, pages
  2483--2491. 2011.

\bibitem{affandi2012markov}
R.H. Affandi, A.~Kulesza, and E.~Fox.
\newblock {M}arkov determinantal point processes.
\newblock In {\em Proceedings of the 28th Conference on Uncertainty in
  Artificial Intelligence (UAI 2012)}, 2012.

\bibitem{gillenwater2014em}
J.A. Gillenwater, A.~Kulesza, E.~Fox, and B.~Taskar.
\newblock {E}xpectation-{M}aximization for learning determinantal point
  processes.
\newblock In {\em Advances in Neural Information Processing Systems 27}, pages
  3149--3157. 2014.

\bibitem{ren2015summarizing}
Z.~Ren and M.~de~Rijke.
\newblock Summarizing contrastive themes via hierarchical non-parametric
  processes.
\newblock In {\em Proceedings of the 38th International ACM SIGIR Conference on
  Research and Development in Information Retrieval (SIGIR 2015)}, pages
  93--102, New York, NY, USA, 2015. ACM.

\bibitem{kulesza2013determinantal}
A.~Kulesza and B.~Taskar.
\newblock Determinantal point processes for machine learning.
\newblock http://arxiv.org/abs/1207.6083, 2013.

\bibitem{liu2010personalized}
J.~Liu, P.~Dolan, and E.R. Pedersen.
\newblock Personalized news recommendation based on click behavior.
\newblock In {\em Proceedings of the 15th International Conference on
  Intelligent User Interfaces (IUI 2010)}, pages 31--40, New York, NY, USA,
  2010. ACM.

\bibitem{li2010contextual}
L.~Li, W.~Chu, J.~Langford, and R.E. Schapire.
\newblock A contextual-bandit approach to personalized news article
  recommendation.
\newblock In {\em Proceedings of the 19th International Conference on World
  Wide Web (WWW 2010)}, pages 661--670, New York, NY, USA, 2010. ACM.

\bibitem{may2012optimistic}
B.C. May, N.~Korda, A.~Lee, and D.S. Leslie.
\newblock Optimistic {B}ayesian sampling in contextual-bandit problems.
\newblock {\em Journal of Machine Learning Research}, 13:2069--2106, 2012.

\bibitem{tang2014ensemble}
L.~Tang, Y.~Jiang, L.~Li, and T.~Li.
\newblock Ensemble contextual bandits for personalized recommendation.
\newblock In {\em Proceedings of the 8th ACM Conference on Recommender Systems
  (RecSys 2014)}, pages 73--80, 2014.

\bibitem{rendle2012factorization}
S.~Rendle.
\newblock Factorization machines with lib{F}{M}.
\newblock {\em ACM Trans. Intell. Syst. Technol.}, 3(3):57:1--57:22, 2012.

\bibitem{juan2016field}
Y.~Juan, Y.~Zhuang, W.-S. Chin, and C.-J. Lin.
\newblock Field-aware factorization machines for {C}{T}{R} prediction.
\newblock In {\em Proceedings of the 10th ACM Conference on Recommender Systems
  (RecSys 2016)}, pages 43--50, 2016.

\bibitem{covington2016deep}
P.~Covington, J.~Adams, and E.~Sargin.
\newblock Deep neural networks for {Y}ou{T}ube recommendations.
\newblock In {\em Proceedings of the 10th ACM Conference on Recommender Systems
  (RecSys 2016)}, pages 191--198, 2016.

\bibitem{contal2014gaussian}
E.~Contal, V.~Perchet, and N.~Vayatis.
\newblock Gaussian process optimization with mutual information.
\newblock In {\em Proceedings of 31st International Conference on Machine
  Learning (ICML 2014)}, 2014.

\bibitem{tak2017data}
H.~Tak and C.~N. Morris.
\newblock Data-dependent posterior propriety of a bayesian beta-binomial-logit
  model.
\newblock {\em Bayesian Analysis}, 12(2):533--555, 2017.

\bibitem{goldberg1998regression}
P.~W. Goldberg, C.~K.~I. Williams, and C.~M. Bishop.
\newblock Regression with input-dependent noise: A {G}aussian process
  treatment.
\newblock In {\em Proceedings of the 1997 Conference on Advances in Neural
  Information Processing Systems 10}, pages 493--499, Cambridge, MA, USA, 1998.
  MIT Press.

\bibitem{rasmussen2005gaussian}
C.~E. Rasmussen and C.~K.~I. Williams.
\newblock {\em Gaussian Processes for Machine Learning (Adaptive Computation
  and Machine Learning)}.
\newblock The MIT Press, 2005.

\bibitem{minka2001expectation}
T.~P. Minka.
\newblock Expectation propagation for approximate {B}ayesian inference.
\newblock In {\em Proceedings of the 17th Conference on Uncertainty in
  Artificial Intelligence (UAI 2001)}, pages 362--369, 2001.

\bibitem{shahbaba2009nonlinear}
B.~Shahbaba and R.~Neal.
\newblock Nonlinear models using {D}irichlet process mixtures.
\newblock {\em Journal of Machine Learning Research}, 10:1829--1850, 2009.

\bibitem{sugiyama2010conditional}
M.~Sugiyama, I.~Takeuchi, T.~Suzuki, T.~Kanamori, H.~Hachiya, and D.~Okanohara.
\newblock Conditional density estimation via least-squares density ratio
  estimation.
\newblock In {\em Proceedings of the 13th International Conference on
  Artificial Intelligence and Statistics (AISTATS 2010)}, volume~9, pages
  781--788, 2010.

\bibitem{ding2010logistic}
N.~Ding and S.v.n. Vishwanathan.
\newblock $t$-logistic regression.
\newblock In {\em Advances in Neural Information Processing Systems 23}, pages
  514--522. 2010.

\bibitem{pericchi1991robust}
L.~R. Pericchi and P.~Walley.
\newblock Robust {B}ayesian credible intervals and prior ignorance.
\newblock {\em International Statistical Review / Revue Internationale de
  Statistique}, 59(1):1--23, 1991.

\bibitem{robbins1952some}
H.~Robbins.
\newblock Some aspects of the sequential design of experiments.
\newblock {\em Bulletin of the American Mathematical Society}, 58:527--535,
  1952.

\bibitem{kahneman1979prospect}
D.~Kahneman and A.~Tversky.
\newblock Prospect theory: An analysis of decision under risk.
\newblock {\em Econometrica}, 47(2):263--291, 1979.

\bibitem{usher2004loss}
M.~Usher and J.~L. McClelland.
\newblock Loss aversion and inhibition in dynamical models of multialternative
  choice.
\newblock {\em Psychological Review}, 111:757--769, 2004.

\bibitem{higgins2000making}
E.~Higgins.
\newblock Making a good decision: Value from fit.
\newblock {\em American Psychologist}, 55(11):1217--1230, 2000.

\bibitem{avnet2006how}
T.~Avnet and E.T. Higgins.
\newblock How regulatory fit impacts value in consumer choices and opinions.
\newblock {\em Journal of Marketing Research}, 43(1):1--10, 2006.

\bibitem{frazier2014incentivizing}
P.~Frazier, D.~Kempe, J.~Kleinberg, and R.~Kleinberg.
\newblock Incentivizing exploration.
\newblock In {\em Proceedings of the 15th ACM Conference on Economics and
  Computation (EC 2014)}, pages 5--22. ACM, 2014.

\bibitem{smith1999diagnosing}
G.E. Smith, M.P. Venkatraman, and R.R. Dholakia.
\newblock Diagnosing the search cost effect: Waiting time and the moderating
  impact of prior category knowledge.
\newblock {\em Journal of Economic Psychology}, 20:285--314, 1999.

\bibitem{kumar2006estimation}
A.~Kumar and T.~Minakshi.
\newblock Estimation of variety seeking for segmentation and targeting: An
  empirical analysis.
\newblock {\em Journal of Targeting, Measurement and Analysis for Marketing},
  15(1):21--29, 2006.

\bibitem{ho2014towards}
E.~Ho and A.~Ilic.
\newblock Towards high resolution identification of variety seeking behavior.
\newblock In {\em Proceedings of the 22nd European Conference on Information
  Systems (ECIS 2014)}, 2014.

\bibitem{raju1980optimum}
P.S. Raju.
\newblock Optimum stimulation level: Its relationship to personality,
  demographics, and exploratory behavior.
\newblock {\em Journal of Consumer Research}, 7(3):272--282, 1980.

\bibitem{kulesza2011learning}
A.~Kulesza and B.~Taskar.
\newblock Learning determinantal point processes.
\newblock In {\em Proceedings of the 27th Conference on Uncertainty in
  Artificial Intelligence (UAI 2011)}, 2011.

\bibitem{chapelle2014simple}
O.~Chapelle, E.~Manavoglu, and R.~Rosales.
\newblock Simple and scalable response prediction for display advertising.
\newblock {\em ACM Trans. Intell. Syst. Technol.}, 5(4):61:1--61:34, 2014.

\end{thebibliography}
