# From bilinear regression to inductive matrix completion:
a quasi-Bayesian analysis

The Tien Mai

(  Department of Mathematical Sciences,
Norwegian University of Science and Technology, Trondheim, Norway.
Email<EMAIL_ADDRESS>)

###### Abstract

In this paper we study the problem of bilinear regression and we further
address the case when the response matrix contains missing data that referred
as the problem of inductive matrix completion. We propose a quasi-Bayesian
approach first to the problem of bilinear regression where a quasi-likelihood
is employed. Then, we adapt this approach to the context of inductive matrix
completion. Under a low-rankness assumption and leveraging PAC-Bayes bound
technique, we provide statistical properties for our proposed estimators and
for the quasi-posteriors. We propose a Langevin Monte Carlo method to
approximately compute the proposed estimators. Some numerical studies are
conducted to demonstrated our methods.

#### Keywords:

bilinear regression, matrix completion, low-rank model, PAC-Bayesian bound,
Langevin Monte Carlo.

##  1 Introduction

In this paper we study the bilinear regression model where a linear
relationship is assumed between a set multiple response variables and two sets
of covariates. Often this model is referred as the growth curve model or
generalized multivariate analysis model that is often used for analyzing
longitudinal data, see for example [45, 47, 40, 48, 31, 27]. However, these
mentioned works cover only the case in which the response matrix is fully
observed. The bilinear regression model with incomplete response is recently
introduced and studied as the so-called inductive matrix completion, a
generalization of the matrix completion problem, [38, 51]. Inductive matrix
completion drawn significant interest in different fields such as drug
repositioning [50], collaborative filtering [23], genomics [38].

Under a low-rank constraint on the coefficient matrix, most approaches for
bilinear regression and inductive matrix completion are frequentist methods
such as maximum likelihood estimation [47] or penalized optimization [38]. One
the other hand, Bayesian approaches have been recently considered discretely
for these problems, for example, the paper [26] proposed a Bayesian approach
for bilinear regression where as a Bayesian method was proposed for inductive
matrix completion in the work [50]. However, unlike frequentist approaches,
statistical properties of the Bayesian approach for these models have not been
carried out yet.

In this work, we aim at filling these gaps by deriving some theoretical
properties of a Bayesian approach for these two problems. More specifically,
we propose a quasi-Bayesian approach to the problem of bilinear regression
where a quasi-likelihood is used. Then, we further generalize this approach to
the problem of inductive matrix completion. Using a spectral scaled Student
prior, we show that the posterior mean satisfy a tight oracle inequality
proving that our method is adaptive to the rank of the coefficient matrix.
Contraction properties of the posteriors are also proved.

Quasi-Bayesian method is an expansion of the Bayesian approach and it is
increasingly popular in statistics and machine learning, as the so-called
generalized Bayesian inference, [29, 4, 19]. This approach allows to relax the
generating model assumption where the likelihood can be replaced by a more
general notion of risk or a quasi-likelihood. Using PAC-Bayesian technique
[37, 46, 8], we derive theoretical results for our proposed quasi-posteriors.
We refer the reader to [20, 1] for recent reviews and introductions to this
technique. PAC-Bayes bounds were successfully used before in the context
matrix estimation problems, see for example [33, 9, 34, 35]. Interestingly,
using this technique for inductive matrix completion, we do not need to put
any assumption on the distribution of the missing entries in the response
matrix. This is in contrast with previous works on matrix completion such as
[25, 5, 30, 17, 28, 39].

The use of a spectral scaled Student prior is motivated from recent works in
which it had been shown to lead to optimal rates in different problems such as
high-dimensional regression [15], image denoising [12] and reduced rank
regression [35]. Although it is not a conjugate prior in our problems, it
allows to implement conveniently gradient-based sampling methods efficiently.
For the computation of the proposed estimators and sampling from the quasi
posterior, we employ a Langevin Monte Carlo (LMC) method.

The rest of the paper is structured as follows. In Section 2, we introduce and
study the problem of bilinear regression, the low-rank promoting prior
distribution is also presented in this section. Section 3 is devoted for the
analysis of the problem of inductive matrix completion. In Section 4, we
present shortly LMC methods for these two problems, after that some
simulations are demonstrated. We conclude our works in Section 5. All of the
technical proofs are given in the appendix A.

###  1.1 Notations

Let $\mathbb{R}^{n_{1}\times n_{2}}$ denote the set of $n_{1}\times n_{2}$
matrices with real coefficients. Let $A^{\intercal}\in\mathbb{R}^{n_{2}\times
n_{1}}$ denote the transpose of $A$. For any $A\in\mathbb{R}^{n_{1}\times
n_{2}}$ and $I=(i,j)\in\\{1,\dots,n_{1}\\}\times\\{1,\dots,n_{2}\\}$, we
denote by $A_{I}=A_{(i,j)}=A_{i,j}$ the coefficient on the $i$-th row and
$j$-th column of $A$. The matrix in $\mathbb{R}^{n_{1}\times n_{2}}$ with all
entries equal to $0$ is denoted by $\mathbf{0}_{n_{1}\times n_{2}}$. For a
matrix $B\in\mathbb{R}^{n_{1}\times n_{1}}$ we let ${\rm Tr}(B)$ denote its
trace. The identity matrix in $\mathbb{R}^{n_{1}\times n_{1}}$ is denoted by
$\mathbf{I}_{n_{1}}$. For $A\in\mathbb{R}^{n_{1}\times n_{2}}$, we define its
sup-norm $\|A\|_{\infty}=\max_{i,j}|A_{i,j}|$; its Frobenius norm $\|A\|_{F}$
is defined by $\|A\|_{F}^{2}={\rm Tr}(A^{\intercal}A)=\sum_{i,j}A_{i,j}^{2}$
and ${\rm rank}(A)$ its rank.

##  2 Bilinear linear regression

###  2.1 Model

Let $Y\in\mathbb{R}^{n\times q}$ consist of $n$ independent response vectors,
$X\in\mathbb{R}^{n\times p}$ be a given between-individuals design matrix and
$Z\in\mathbb{R}^{k\times q}$ is a known within-individuals design matrix.
Consider the bilinear regression model as follows

$Y=XM^{*}Z+E,$ (1)

where $M^{*}\in\mathbb{R}^{p\times k}$ is the unknown parameter matrix. The
random noise matrix $E$ is assumed to have zero mean, $\mathbb{E}(E)=0.$ The
main assumption here is the reduced-rank restriction on the model parameter
that ${\rm rank}(M^{*})<\min(p,k)$.

In model (1), as soon as $k=q$ and $Z=\mathbf{I}_{q}$ we recover, as a special
case, the reduced rank regression problem [3, 24]. However, we consider in
this paper the case where $Z$ contains explanatory variables. The low-rank
assumption can be interpreted as there exists a latent process so that the
response data is not only affected by the “between-individuals” structure of
the model but can also be interpretable in terms of the “within-individuals”
structure. This model is often called the growth curve model or generalized
multivariate analysis model (GMANOVA), see [47] for more details.

######  Assumption 1.

There is a known constant $C<+\infty$ such that $\|XM^{*}Z\|_{\infty}\leq C.$

From Assumption 1, it is not reliable to return predictions $XMZ$ with entries
that are outside of interval $[-C,C]$. However, for computational reason, it
is extremely convenient to employ an unbounded prior for $M$. Therefore, we
propose to use unbounded distributions for $M$, but to use, as a predictor, a
truncated version of $XMZ$ rather than $M$ itself. For a matrix $A$, let
$\Pi_{C}(A)=\arg\min_{\|B\|_{\infty}\leq C}\|A-B\|_{F}$ be the orthogonal
projection of $A$ on matrices with entries bounded by $C$. Note that $B$ is
simply obtained by replacing entries of $A$ larger than $C$ by $C$, and
entries smaller than $-C$ by $-C$.

For a matrix $M\in\mathbb{R}^{p\times k}$, we denote by $r(M)$ the empirical
risk of $M$ as

$r(M)=\dfrac{1}{nq}\|Y-\Pi_{C}(XMZ)\|_{F}^{2}$

and its expectation is denoted by

$R(M)=\mathbb{E}\left[r(M)\right]=\mathbb{E}\left[\left(Y_{11}-(\Pi_{C}(XMZ))_{{}_{11}}\right)^{2}\right].$

The focus of our work in this paper is on the predictive aspects of the model
that is a matrix $M$ predicts almost as well as $M^{*}$ if $R(M)-R(M^{*})$ is
small. Under the assumption that $E_{ij}$ has a finite variance, using
Pythagorean theorem, we have

$R(M)-R(M^{*})=\dfrac{1}{nq}\|\Pi_{C}(XMZ)-XM^{*}Z\|^{2}_{F}$ (2)

for any $M$, which means that our results can also be interpreted in terms of
Frobenius norm.

Let $\pi$ be a prior distribution on $\mathbb{R}^{p\times k}$ (see Subsection
2.2). For any $\lambda>0$, we define the quasi-posterior

$\widehat{\rho}_{\lambda}(dM)\propto\exp(-\lambda r(M))\pi(dM).$

Note that, for $\lambda=nq/(2\sigma^{2})$, this is exactly the posterior that
we would be obtain for a Gaussian noise $E_{ij}\sim\mathcal{N}(0,\sigma^{2})$.
However, our theoretical results will hold under a more general class of
noise. It is known that a small enough $\lambda$ serves the purpose when the
model is misspecified [19]. Moreover, even in the case of a Gaussian noise, in
high-dimensional settings, taking $\lambda$ smaller than
$\lambda=n/(2\sigma^{2})$ leads to better adaptation properties [10, 15]. The
choice of $\lambda$ in our method will be specified below.

We consider the following posterior mean of $XMZ$, given by

$\widehat{XMZ}_{\lambda}=\int\Pi_{C}(XMZ)\widehat{\rho}_{\lambda}(dM).$ (3)

Let us briefly discuss the above estimator. It is noted from our simulations
experiments that using reasonable values for $C$, the Monte Carlo algorithm
never samples matrices $M$ such that $\Pi_{C}(XMZ)\neq XMZ$. In other words,
it has very little impact in practice and it is necessary for technical
proofs. If one would like to obtain an estimator of $M^{*}$ instead of an
estimator of $XM^{*}Z$, when $X^{\intercal}X$ and $ZZ^{\intercal}$ are
invertible, one can consider
$\widehat{M}_{\lambda}=(X^{\intercal}X)^{-1}X^{\intercal}\widehat{XMZ}_{\lambda}Z^{\intercal}(ZZ^{\intercal})^{-1}$
and note that $X\widehat{M}_{\lambda}Z=\widehat{XMZ}_{\lambda}$.

The quasi-posterior is often referred to as the so-called “Gibbs posterior” in
the PAC-Bayes approach [7, 8, 2, 20, 1]. Moreover, $\widehat{M}_{\lambda}$ is
sometimes referred to as the Gibbs estimator, or the exponentially weighted
aggregate (EWA) [10, 42, 13].

###  2.2 Prior specification

We consider, in this paper, the following spectral scaled Student prior, with
parameter $\tau>0$,

$\displaystyle\pi(M)\propto\det(\tau^{2}\mathbf{I}_{m}+MM^{\intercal})^{-(p+m+2)/2}.$
(4)

One can verify that
$\pi(M)\propto\prod_{j=1}^{m}(\tau^{2}+s_{j}(M)^{2})^{-(p+m+2)/2},$ where
$s_{j}(M)$ denotes the $j^{th}$ largest singular value of $M$. Intuitively,
one can recognize a scaled Student distribution evaluated at $s_{j}(M)$ in the
last display above which induces approximate sparsity on the $s_{j}(M)$ [15].
Also, the log-sum function $\sum_{j=1}^{m}\log(\tau^{2}+s_{j}(M)^{2})$ used by
[6, 49] is well known to enforce approximate sparsity on the singular values
$s_{j}(M)$. In other words, under this prior, most of the $s_{j}(M)$ are close
to $0$, which means that $M$ is well approximated by a low-rank matrix. Thus
it has the ability to promote the low-rankness of $M$.

The spectral scaled Student prior has been used before in the context of image
denoising [12] and in reduced rank regression [35]. It can also be seen as the
marginal distribution of the Gaussian-inverse Wishart prior where the
precision matrix is integrated out that is explored in [49] for matrix
completion. As mentioned in the introduction, this prior is not conjugate in
our problems, however, it is particularly convenient to implement gradient-
based sampling algorithms, the Langevin Monte Carlo, see Section 4.

###  2.3 Theoretical results

We assume the sub-exponential distribution assumption on the noise.

######  Assumption 2.

The entries $E_{i,j}$ of $E$ are independent. There exist two known constants
$\sigma>0$ and $\xi>0$ such that $\forall k\geq
2,\,\mathbb{E}(|E_{i,j}|^{k})\leq\sigma^{2}k!\xi^{k-2}/2.$

Let us put
$C_{1}=8(\sigma^{2}+C^{2});\,C_{2}=64C\max(\xi,C);\,\tau^{*}=\sqrt{C_{1}(k+p)/(nkq\|X\|_{F}^{2}\|Z\|_{F}^{2})}.$
The statistical properties of mean estimator are given in the following
theorem where we propose a non-asymptotic analysis for our mean estimator.

######  Theorem 1.

Let Assumptions 1 and 2 be satisfied. Fix the parameter $\tau=\tau^{*}$ in the
prior. Fix $\delta>0$ and define
$\lambda^{*}:=nq\min(1/(2C_{2}),\delta/[C_{1}(1+\delta)])$. Then, for any
$\varepsilon\in(0,1)$, we have, with probability at least $1-\varepsilon$ on
the sample,

$\left\|\widehat{XMZ}_{\lambda^{*}}-XM^{*}Z\right\|_{F}^{2}\leq\inf_{0\leq
r\leq pk}\inf_{\tiny\begin{array}[]{c}\bar{M}\in\mathbb{R}^{p\times k}\\\ {\rm
rank}(\bar{M})\leq
r\end{array}}\Biggl{\\{}(1+\delta)\left\|X\bar{M}Z-XM^{*}Z\right\|_{F}^{2}\,+\\\
\frac{C_{1}(1+\delta)^{2}}{\delta}\left[4r(k+p+2)\log\left(1+\frac{\|X\|_{F}\|Z\|_{F}\|\bar{M}\|_{F}}{\sqrt{C_{1}}}\sqrt{\frac{nkq}{r(k+p)}}\right)+k+p+2\log\frac{2}{\varepsilon}\right]\Biggr{\\}}.$

The choice of $\lambda=\lambda^{*}$ arises from the optimization of a upper
bound on the risk $R$ (in the proof of this theorem). However this choice may
not be necessarily the best choice in practice even though it gives the good
order of magnitude for $\lambda$. The user could use cross-validation to
properly tune the temperature parameter. It is also noted that ${\rm
rank}(\bar{M})\neq 0$ is not required in the above formula, if ${\rm
rank}(\bar{M})=0$ then $\bar{M}=0$ and we interpret $0\log(1+0/0)$ as $0$. The
proof of this theorem relies on the PAC-Bayes theory and it is provided in the
appendix A. In particular, we can upper bound the infimum on $\bar{M}$ by
taking $\bar{M}=M^{*}$, which leads to the following result.

######  Corollary 1.

Under the same assumptions and the same $\tau,\lambda^{*}$ as in Theorem 1,
let $r^{*}={\rm rank}(M^{*})$. Then, for any $\varepsilon\in(0,1)$, we have,
with probability at least $1-\varepsilon$ on the sample,

$\left\|\widehat{XMZ}_{\lambda^{*}}-XM^{*}Z\right\|_{F}^{2}\leq\frac{C_{1}(1+\delta)^{2}}{\delta}\left[4r^{*}(k+p+2)\log\left(1+\frac{\|X\|_{F}\|Z\|_{F}\|M^{*}\|_{F}}{\sqrt{C_{1}}}\sqrt{\frac{nkq}{r(k+p)}}\right)+\right.\\\
\left.k+p+2\log\frac{2}{\varepsilon}\right].$

While Theorem 1 is on the statistical properties of the posterior mean, we
further show the contraction of the posterior in the following theorem.

######  Theorem 2.

Under the assumptions for Theorem 1, let $\varepsilon_{n}$ be any sequence in
$(0,1)$ such that $\varepsilon_{n}\rightarrow 0$ when $n\rightarrow\infty$.
Define

$\mathcal{M}_{n}=\Biggl{\\{}M\in\mathbb{R}^{p\times
k}:\left\|\widehat{XMZ}_{\lambda^{*}}-XM^{*}Z\right\|_{F}^{2}\leq\inf_{0\leq
r\leq pk}\inf_{\tiny\begin{array}[]{c}\bar{M}\in\mathbb{R}^{p\times k}\\\ {\rm
rank}(\bar{M})\leq
r\end{array}}\Biggl{\\{}(1+\delta)\left\|X\bar{M}Z-XM^{*}Z\right\|_{F}^{2}\,+\\\
\frac{C_{1}(1+\delta)^{2}}{\delta}\left[4r(k+p+2)\log\left(1+\frac{\|X\|_{F}\|Z\|_{F}\|\bar{M}\|_{F}}{\sqrt{C_{1}}}\sqrt{\frac{nkq}{r(k+p)}}\right)+k+p+2\log\frac{2}{\varepsilon_{n}}\right]\Biggr{\\}}.$

Then

$\mathbb{E}\Bigl{[}\mathbb{P}_{M\sim\widehat{\rho}_{\lambda}}(M\in\mathcal{M}_{n})\Bigr{]}\geq
1-\varepsilon_{n}\xrightarrow[n\rightarrow\infty]{}1.$

The proof of this theorem is provided in the appendix A.

##  3 Inductive matrix completion

###  3.1 Model and method

In the context of inductive matrix completion, given two side information
matrices $X$ and $Z$, we assume that only a random subset $\Omega$ of the
entries of $Y$ in model (1) is observed. More precisely, we assume that we
observe $m$ i.i.d random pairs
$(\mathcal{I}_{1},Y_{1}),\ldots,(\mathcal{I}_{m},Y_{m})$ given by

$Y_{i}=(XM^{*}Z)_{\mathcal{I}_{i}}+\mathcal{E}_{i},\quad i=1,\ldots,m$ (5)

where $M^{*}\in\mathbb{R}^{p\times k}$ is the unknown parameter matrix
expected to be low-rank. The noise variables $\mathcal{E}_{i}$ are assumed to
be independent with $\mathbb{E}(\mathcal{E}_{i})=0.$ The variables
$\mathcal{I}_{i}$ are i.i.d copies of a random variable $\mathcal{I}$ having
distribution $\Pi$ on the set $\\{1,\ldots,n\\}\times\\{1,\ldots,q\\}$, we
denote $\Pi_{x,y}:=\Pi(\mathcal{I}=(x,y))$.

When $p=n,k=q$ and that $X=\mathbf{I}_{n},Z=\mathbf{I}_{q}$ are the identity
matrices, we recover the matrix completion as a special case [30]; whereas
when $k=q$ and that $Z=\mathbf{I}_{q}$ is the identity matrix, we recover the
reduced rank regression problem with incomplete response as a special case
[32, 35]. However, in the context of inductive matrix completion, we focus in
the cases where $X$ and $Z$ contains explanatory variables. In order words,
the side information from $n$ users and the side information from the $q$
items will be taken into account in this model [51].

It is also noted that there are two approaches can be used to model the
observed entries of $Y$: with, or without replacement. For matrix completion,
both had been studied, see for example [5] for the case without replacement
and [30] with replacement. Both settings have their own practical
applications, and the same estimation methods are used in both cases. In this
paper, we study the case for i.i.d variables $\mathcal{I}_{i}$, which means
that it is possible to observe the same entry multiple times. Note that, from
the results in Section 6 of [22], our results can be extended directly to the
case of sampling $\mathcal{I}_{i}$ without replacement on the condition that
they are sampled uniformly, and that there is no observation noise:
$\mathcal{E}_{i}=0$.

We are now adapting the quasi-Bayesian approach for bilinear regression in
Section 2 to the context of inductive matrix completion. For a probability
distribution $P$ on $\\{1,\dots,n_{1}\\}\times\\{1,\dots,n_{2}\\}$, we
generalize the Frobenius norm by
$\|A\|_{F,P}^{2}=\sum_{i,j}P[(i,j)]A_{i,j}^{2}$; note that when $P$ is the
uniform distribution, then $\|A\|_{F,P}^{2}=\|A\|_{F}^{2}/(n_{1}n_{2})$.

For a matrix $M\in\mathbb{R}^{p\times k}$, we denote the empirical risk of
$M$, $r^{\prime}(M)$, and its expected risk $R^{\prime}(M)$ respectively as

$r^{\prime}(M)=\frac{1}{m}\sum_{i=1}^{m}\left(Y_{i}-(\Pi_{C}(XMZ))_{\mathcal{I}_{i}}\right)^{2},\quad
R^{\prime}(M)=\mathbb{E}[r^{\prime}(M)]=\mathbb{E}[\left(Y_{1}-(\Pi_{C}(XMZ))_{\mathcal{I}_{1}}\right)^{2}].$

As in Section 2, we will focus on the predictive aspects of the model that is
a matrix $M$ predicts almost as well as $M^{*}$ if
$R^{\prime}(M)-R^{\prime}(M^{*})$ is small. Under the assumption that
$\mathcal{E}_{i}$ has a finite variance, thanks to the Pythagorean theorem, we
have

$R^{\prime}(M)-R^{\prime}(M^{*})=\|\Pi_{C}(XMZ)-XM^{*}Z\|^{2}_{F,\Pi}$ (6)

for any $M$, which means that our results can also be interpreted in terms of
estimation of $M^{*}$ with respect to a generalized Frobenius norm.

Here, the prior $\pi$ is the low-rank inducing prior specified in the
Subsection 2.2 above. For any $\lambda>0$, we define the quasi-posterior

$\hat{\rho}^{\prime}_{\lambda}(dM)\propto\exp(-\lambda r^{\prime}(M))\pi(dM).$

We will actually specify our choice of $\lambda$ below.

The truncated posterior mean of $XMZ$ is given by

$\widehat{XMZ}_{\lambda}=\int\Pi_{C}(XM)\hat{\rho}^{\prime}_{\lambda}(dM).$
(7)

Here, for the same technical reasons as in the context of bilinear regression
problem, this truncation has a very little impact in practice for reasonable
values of $C$.

###  3.2 Theoretical results

In this section, we derive the statistical properties of the posterior
$\hat{\rho}^{\prime}_{\lambda}$ and the mean estimator
$\widehat{XMZ}_{\lambda}$ for the context of inductive matrix completion. Let
us first state our assumptions on this model.

######  Assumption 3.

The noise variables $\mathcal{E}_{1},\ldots,\mathcal{E}_{m}$ are independent
and independent of $\mathcal{I}_{1},\ldots,\mathcal{I}_{m}$. There exist two
known constants $\sigma^{\prime}>0$ and $\xi^{\prime}>0$ such that $\forall
k\geq 2,\,\mathbb{E}(|\mathcal{E}_{i}|^{k})\leq\sigma^{\prime 2}k!\xi^{\prime
k-2}/2.$

Assumption 1 and 3 are both standard, they have been used in [32] for
theoretical analysis of reduced rank regression and in [30] for trace
regression and matrix completion.

Let us put $C^{\prime}_{1}=8(\sigma^{\prime
2}+C^{2});\,C^{\prime}_{2}=64C\max(\xi^{\prime},C);\,\tau^{*}=\sqrt{C^{\prime}_{1}(k+p)/(mkp\|X\|_{F}^{2}\|Z\|_{F}^{2})}.$

######  Theorem 3.

Let Assumptions 1 and 3 be satisfied. Fix the parameter $\tau=\tau^{*}$ in the
prior. Fix $\delta>0$ and define
$\lambda^{\prime*}:=m\min(1/(2C^{\prime}_{2}),\delta/[C^{\prime}_{1}(1+\delta)])$.
Then, for any $\varepsilon\in(0,1)$, we have, with probability at least
$1-\varepsilon$ on the sample,

$\left\|\widehat{XMZ}_{\lambda^{\prime*}}-XM^{*}Z\right\|_{F,\Pi}^{2}\leq\inf_{0\leq
r\leq pk}\inf_{\tiny\begin{array}[]{c}\bar{M}\in\mathbb{R}^{p\times k}\\\ {\rm
rank}(\bar{M})\leq
r\end{array}}\Biggl{\\{}(1+\delta)\|X\bar{M}Z-XM^{*}Z\|_{F,\Pi}^{2}+\\\
\frac{C^{\prime}_{1}(1+\delta)^{2}}{\delta}\frac{\left(4r(k+p+2)\log\left(1+\frac{\|X\|_{F}\|Z\|_{F}\|\bar{M}\|_{F}}{\sqrt{C_{1}}}\sqrt{\frac{mkp}{r(k+p)}}\right)+k+p+2\log\frac{2}{\varepsilon}\right)}{m}\Biggr{\\}}.$

Similar to the context bilinear regression, the choices of
$\lambda=\lambda^{*},\tau=\tau^{*}$ come from the optimization of a upper
bound on the risk $R$ (in the proof of this theorem). Therefore these choices
may not be necessarily the best choice in practice even though it gives the
good order of magnitude for tuning these parameters. The user could use cross-
validation to properly tune them in practice. Note again that ${\rm
rank}(\bar{M})\neq 0$ is not required in the above formula, if ${\rm
rank}(\bar{M})=0$ then $\bar{M}=0$ and we interpret $0\log(1+0/0)$ as $0$. The
proof of this theorem is provided in the appendix A. In particular, we can
upper bound the infimum on $\bar{M}$ by taking $\bar{M}=M^{*}$, which leads to
the following result.

######  Corollary 2.

Under the assumptions that Theorem 3 holds, let $r^{*}={\rm rank}(M^{*})$. Put

$R_{\delta,m,p,k,r^{*},\varepsilon}:=\frac{C^{\prime}_{1}(1+\delta)^{2}}{\delta}\frac{\left(4r(k+p+2)\log\left(1+\frac{\|X\|_{F}\|Z\|_{F}\|\bar{M}\|_{F}}{\sqrt{C_{1}}}\sqrt{\frac{mkp}{r(k+p)}}\right)+k+p+2\log\frac{2}{\varepsilon}\right)}{m},$

then

$\displaystyle\left\|\widehat{XMZ}_{\lambda^{\prime*}}-XM^{*}Z\right\|_{F,\Pi}^{2}\leq
R_{\delta,m,p,k,r^{*},\varepsilon}$

and in particular, if the sampling distribution $\Pi$ is uniform,

$\displaystyle\dfrac{\|\widehat{XMZ}_{\lambda^{\prime*}}-XM^{*}Z\|^{2}_{F}}{nq}\leq
R_{\delta,m,p,k,r^{*},\varepsilon}.$

######  Remark 1.

Up to a log-term, our error rate $r(k+p)/m$ is similar to the best known up-
to-date rate derived in [51].

While Theorem 3 is about the finite sample convergence rate of the posterior
mean, it is actually possible to prove that the quasi-posterior
$\hat{\rho}^{\prime}_{\lambda}$ contracts around $M^{*}$ at the same rate.

######  Theorem 4.

Under the same assumptions for Theorem 3, and the same definition for $\tau$
and $\lambda^{*}$, let $\varepsilon_{m}$ be any sequence in $(0,1)$ such that
$\varepsilon_{m}\rightarrow 0$ when $m\rightarrow\infty$. Define

$\mathcal{E}_{m}=\Biggl{\\{}M\in\mathbb{R}^{p\times
k}:\|\Pi_{C}(XMZ)-XM^{*}Z\|^{2}_{F,\Pi}\leq\inf_{1\leq r\leq
pk}\inf_{\tiny\begin{array}[]{c}\bar{M}\in\mathbb{R}^{p\times k}\\\ {\rm
rank}(\bar{M})\leq
r\end{array}}\Biggl{[}(1+\delta)\|X\bar{M}Z-XM^{*}Z\|_{F,\Pi}^{2}+\\\
\frac{C^{\prime}_{1}(1+\delta)^{2}}{\delta}\frac{\left(4r(k+p+2)\log\left(1+\frac{\|X\|_{F}\|Z\|_{F}\|\bar{M}\|_{F}}{\sqrt{C_{1}}}\sqrt{\frac{mkp}{r(k+p)}}\right)+k+p+2\log\frac{2}{\varepsilon_{m}}\right)}{m}\Biggr{]}\Biggr{\\}}.$

Then

$\mathbb{E}\Bigl{[}\mathbb{P}_{M\sim\hat{\rho}^{\prime}_{\lambda}}(M\in\mathcal{E}_{m})\Bigr{]}\geq
1-\varepsilon_{m}\xrightarrow[m\rightarrow\infty]{}1.$

The proof of this theorem is provided in the appendix A.

##  4 Numerical studies

###  4.1 Langevin Monte Carlo implementation

In this section, we propose to sample from the (quasi) posterior, in Section 2
and 3, by a suitable version of the Langevin Monte Carlo (LMC) algorithm, a
gradient-based sampling method. We propose use a constant step-size unadjusted
LMC algorithm, see [16] for more details. The algorithm is given by an initial
matrix $M_{0}$ and the recursion

$M_{k+1}=M_{k}-h\nabla\log\widehat{\rho}_{\lambda}(M_{k})+\sqrt{2h}\,N_{k}\qquad
k=0,1,\ldots$ (8)

where $h>0$ is the step-size and $N_{0},N_{1},\ldots$ are independent random
matrices with i.i.d standard Gaussian entries. We provide a pseudo-code for
LMC in Algorithm 1. For small values of the step-size $h$, the output of
Algorithm 1, $\widehat{M}$, is very close to the integral (3) of interest.
However, for some $h$ that may not small enough, the sum can explode [44]. In
such case, we consider to include a Metropolis–Hastings correction in the
algorithm. Another possible choice is to take a smaller $h$ and restart the
algorithm, although it slows down the algorithm but we keep some control on
its time of execution. On the other hand, the Metropolis–Hastings approach
ensures the convergence to the desired distribution, however, the algorithm is
greatly slowed down because of an additional acception/rejection step at each
iteration.

Next, we propose a Metropolis-Hasting correction to the LMC algorithm. It
guarantees the convergence to the (quasi) posterior and it also provides a
useful way for choosing $h$. More precisely, we consider the update rule in
(8) as a proposal for a new candidate,

$\displaystyle\tilde{M}_{k+1}=M_{k}-h\nabla\log\widehat{\rho}_{\lambda}(M_{k})+\sqrt{2h}\,N_{k},\qquad
k=0,1,\ldots,$ (9)

Note that the matrix $\tilde{M}_{k+1}$ is normally distributed with mean
$M_{k}-h\nabla\log\widehat{\rho}_{\lambda}(M_{k})$ and the covariance matrices
equal to $2h$ times the identity matrices. This proposal is then accepted or
rejected according to the Metropolis-Hastings algorithm that the proposal is
accepted with probability:

$A_{MALA}:=\min\left\\{1,\frac{\widehat{\rho}_{\lambda}(\tilde{M}_{k+1})q(M_{k}|\tilde{M}_{k+1})}{\widehat{\rho}_{\lambda}(M_{k})q(\tilde{M}_{k+1}|M_{k})}\right\\},$
(10)

where

$q(x^{\prime}|x)\propto\exp\left(-\frac{1}{4h}\|x^{\prime}-x+h\nabla\log\widehat{\rho}_{\lambda}(x)\|^{2}_{F}\right)$

is the transition probability density from $x$ to $x^{\prime}$. The details of
the Metropolis-adjusted Langevin algorithm (denoted by MALA) are presented
Algorithm 2. Compared to random-walk Metropolis–Hastings, MALA usually
proposes moves into regions of higher probability, which are then more likely
to be accepted.

We note that the step-size $h$ for MALA is chosen such that the acceptance
rate is approximate $0.5$ following [43], while the step-size for LMC in the
same setting should be smaller than the one for MALA [11].

###  4.2 Simulation studies for biliear regression

We perform some numerical studies on simulated data to access the performance
of our proposed algorithms. All simulations were conducted using the R
statistical software [41].

For fixed dimensions $q=10,k=20$ of the data, we vary $n=100$ and $n=1000$ to
check the effect of the samples whereas the dimensions of the coefficient
matrix are varied by $p=10$ and $p=100$. The entries of the design matrices
$X$ and $Z$ are independently simulated from the standard Gaussian
$\mathcal{N}(0,1)$. Then, given a matrix $M^{*}$, we simulate the response
matrix $Y$ from model (1) whose entries of the noise matrix $E$ are i.i.d
sampled from $\mathcal{N}(0,1)$. We consider the following setups for the true
coefficient matrix:

  * •

Model I: The true coefficient matrix $M^{*}$ is a rank-2 matrix that is
generated as $M^{*}=B_{1}B_{2}^{\top}$ where $B_{1}\in\mathbb{R}^{p\times
2},B_{2}\in\mathbb{R}^{k\times 2}$ and all entries in $B_{1}$ and $B_{2}$ are
i.i.d sampled from $\mathcal{N}(0,1)$.

  * •

Model II: An approximate low-rank set up is studied. This series of simulation
is similar to the Model I, except that the true coefficient is no longer rank
2, but it can be well approximated by a rank 2 matrix:

$M^{*}=2\cdot B_{1}B_{2}^{\top}+U,$

where $U$ is a matrix whose entries are i.i.d sampled from
$\mathcal{N}(0,0.1)$.

We compare our approaches denoted by LMC and MALA against the (generalized)
ordinary least square [47], denoted by OLS. The OLS is defined as follow

$\hat{M}_{\rm
OLS}=(X^{\top}X)^{\dagger}X^{\top}YZ^{\top}(ZZ^{\top})^{\dagger}$

where $A^{\dagger}$ denotes the Moore-Penrose matrix pseudo-inversion. We
fixed $\lambda=nq,\tau=1$ and the LMC and MALA methods are initiated at the
OLS estimator and are run with 10000 iterations where the first 1000 steps are
removed as burn-in periods.

The evaluations are done by using the mean squared estimation error (Est) and
the normalized (relative) mean square error (Nmse)

${\rm Est}:=\|\hat{M}-M^{*}\|^{2}_{F}/(pk),\quad{\rm
Nmse}:=\|\hat{M}-M^{*}\|^{2}_{F}/\|M^{*}\|^{2}_{F},$

and the prediction error (Pred) as

${\rm Pred}:=\|X(\hat{M}-M^{*})Z\|^{2}_{F}/(nq),$

where $\hat{M}$ here is one of the estimator from LMC, MALA or OLS. We report
the averages and the standard deviation of these errors over 100 data
replications.

The results are given in Table 2 and Table 3. We can observe that our proposed
methods behave quite similar to OLS method. However, the estimation method
obtained from MALA algorithm often yields smaller prediction errors,
especially in high-dimensional settings. This feature is even more visible in
the context of inductive matrix completion in the next subsection.

###  4.3 Simulation studies for inductive matrix completion

The simulation settings for inductive matrix completion are similar to the
settings for bilinear regression, subsection 4.2. However, after obtaining the
response matrix $Y$, we remove uniformly at random $\kappa=10\%$ and
$\kappa=30\%$ of the entries of $Y$. Here $\kappa$ denotes the missing rate.
We denote the response matrix with missing entries by $Y_{\rm miss}$.

As in the context of inductive matrix completion, we only observe the response
matrix with missing entries, $Y_{\rm miss}$, and thus we can not construct the
OLS estimator as in the case of bilinear regression. For this purpose, we
first impute the missing entries in $Y_{\rm miss}$ by using the R package
softImpute [21] where the rank of $M^{*}$ is specified as the true rank for
matrix $Y_{\rm miss}$. We denote the resulting imputed matrix by $Y_{\rm
imp}$.

We compare our approaches denoted by LMC and MALA against the (imputed and
generalized) ordinary least square, denoted by OLS_imp. The OLS_imp is defined
as follow

$\hat{M}_{\rm OLS\\_imp}=(X^{\top}X)^{\dagger}X^{\top}Y_{\rm
imp}Z^{\top}(ZZ^{\top})^{\dagger}$

where $A^{\dagger}$ denotes the Moore-Penrose matrix pseudo-inversion. The LMC
and MALA methods are initiated at the OLS_imp estimator and are run with 10000
iterations where the first 1000 steps are removed as burn-in periods.

As in subsection 4.2, we report the averages and the standard deviation of the
the mean squared estimation error (Est), the normalized (relative) mean square
error (Nmse) and of the prediction error (Pred) over 100 data replications.

The results are given in Table 1 and Table 4. Here, it is more clear that the
results from our MALA method outperforms the other methods in term of
prediction error in most considered settings. It is more significant when the
missing rate in the response matrix increases. It is also noted that our MALA
method robustly works well in the case of approximate low-rank setting (model
II) while the OLS and LMC are not.

##  5 Discussion and conclusion

In this work, we presented a unified quasi-Bayesian analysis for the problem
of bilinear regression and for the problem of inductive matrix completion. We
focused on filling the gaps in theoretical understanding of Bayesian-like
methods in these problem. By leveraging PAC-Bayesian bound techniques, we
proved some non-asymptotic results for our proposed estimator, we also derived
the contraction properties for the corresponding (quasi) posteriors.

We further propose an efficient gradient-based sampling algorithm to sample
from the (quasi) posterior and to approximately compute the mean estimators.
The resulting methods, denoted by LMC and MALA, were shown to work well in
some simulation studies and in compared with a ordinary least squared method.

Some questions remain open: such as missing entries in the covariate matrices
$X,Z$ or relaxing the i.i.d assumption by considering models with dispersion
matrix. These might be the objects of future works

### Acknowledgments

TTM is supported by the Norwegian Research Council, grant number 309960
through the Centre for Geophysical Forecasting at NTNU. The R codes and data
used in the numerical experiments are available at:
https://github.com/tienmt/blr_imc .

### Conflicts of interest/Competing interests

The authors declare no potential conflict of interests.

Table 1: Simulation results on simulated data in Model I in inductive matrix completion, with fixed $q=10,k=20$, for different methods, with their standard error in parentheses over 100 replications. ($\kappa$ is the missing rate; Est: average of estimation errors; Pred: average of prediction errors; Nmse: average of normalized estimation errors). |  |  Errors | LMC | MALA | OLS_imp |
---|---|---|---|---|---|---
|  |  Est | 1.0559 (.5060) | 1.0803 (.5122) | 1.0559 (.5060) |
$n=100$ |  $p=10$ |  Pred | 0.1028 (.0193) | 0.1082 (.0143) |  0.1020 (.0197) |
$\kappa=10\%$ |  |  Nmse | 0.4986 (.1116) | 0.5139 (.1197) | 0.4986 (.1116) |
|  |  Est | 1.4008 (.8555) |  1.3987 (.8542) | 1.4009 (.8555) |
|  $p=100$ |  Pred |  1.2250 (.4568) | 1.4468 (.4137) | 1.2252 (.4570) |
|  |  Nmse | 0.7148 (.3591) |  0.7136 (.3581) | 0.7148 (.3591) |
|  |  Est | 1.0432 (.4963) | 1.0917 (.5085) | 1.0432 (.4963) |
$n=100$ |  $p=10$ |  Pred | 0.2402 (.2705) |  0.1447 (.0204) | 0.2446 (.2780) |
$\kappa=30\%$ |  |  Nmse | 0.5242 (.1257) | 0.5538 (.1335) | 0.5242 (.1257) |
|  |  Est | 1.6242 (.8179) |  1.6224 (.8169) | 1.6242 (.8179) |
|  $p=100$ |  Pred |  9.8879 (14.11) | 10.807 (13.84) | 9.8901 (14.11) |
|  |  Nmse | 0.7993 (.3340) |  0.7985 (.3334) | 0.7993 (.3340) |
|  |  Est | 0.9810 (.4532) | 0.9882 (.4478) | 0.9810 (.4532) |
$n=1000$ |  $p=10$ |  Pred | 0.0114 (.0033) |  0.0112 (.0015) | 0.0114 (.0033) |
$\kappa=10\%$ |  |  Nmse | 0.4933 (.1076) | 0.4984 (.1075) | 0.4933 (.1076) |
|  |  Est | 1.0063 (.3465) | 1.0088 (.3471) | 1.0063 (.3465) |
|  $p=100$ |  Pred | 0.1902 (.1758) |  0.1116 (.0049) | 0.1902 (.1759) |
|  |  Nmse | 0.5069 (.1049) | 0.5082 (.1050) | 0.5069 (.1049) |
|  |  Est | 1.0110 (.4886) | 1.0223 (.4872) | 1.0110 (.4886) |
$n=1000$ |  $p=10$ |  Pred | 0.0539 (.0599) |  0.0141 (.0019) | 0.0540 (.0599) |
$\kappa=30\%$ |  |  Nmse | 0.5129 (.1030) | 0.5206 (.1043) | 0.5129 (.1030) |
|  |  Est | 1.0291 (.3567) | 1.0312 (.3555) | 1.0291 (.3567) |
|  $p=100$ |  Pred | 1.7529 (1.914) |  0.1475 (.0078) | 1.7530 (1.913) |
|  |  Nmse | 0.5054 (.1055) | 0.5067 (.1053) | 0.5054 (.1055) |
Table 2: Simulation results on simulated data in Model I in bilinear regression, with fixed $q=10,k=20$, for different methods, with their standard error in parentheses over 100 replications. (Est: average of estimation errors; Pred: average of prediction errors; Nmse: average of normalized estimation errors). |  |  Errors | LMC | MALA | OLS |
---|---|---|---|---|---|---
|  |  Est | 1.0053 (.5480) | 1.0342 (.5559) |  1.0052 (.5478) |
$n=100$ |  $p=10$ |  Pred | 0.1138 (.0171) |  0.0985 (.0151) | 0.1014 (.0154) |
|  |  Nmse | 0.4931 (.1178) | 0.5100 (.1207) |  0.4930(.1178) |
|  |  Est | 1.3544 (.5867) |  1.3384 (.5836) | 1.3544 (.5867) |
|  $p=100$ |  Pred | 1.0066 (.0430) |  0.8761 (.0756) | 1.0030 (.0424) |
|  |  Nmse | 0.7049 (.2944) |  0.6963 (.2927) | 0.7049 (.2944) |
|  |  Est | 1.0776 (.5671) | 1.0900 (.5670) | 1.0776 (.5671) |
$n=1000$ |  $p=10$ |  Pred | 0.0099 (.0013) | 0.0099 (.0013) | 0.0099 (.0013) |
|  |  Nmse | 0.5185 (.1198) | 0.5264 (.1219) | 0.5185 (.1198) |
|  |  Est | 0.9662 (.3240) | 0.9688 (.3244) | 0.9662 (.3240) |
|  $p=100$ |  Pred | 0.0999 (.0051) |  0.0989 (.0049) | 0.0998 (.0051) |
|  |  Nmse | 0.4961 (.1183) | 0.4976 (.1191) | 0.4961 (.1183) |
Table 3: Simulation results on simulated data in Model II (approximate low-rank) in bilinear regression, with fixed $q=10,k=20$, for different methods, with their standard error in parentheses over 100 replications. (Est: average of estimation errors; Pred: average of prediction errors; Nmse: average of normalized estimation errors). |  |  Errors | LMC | MALA | OLS |
---|---|---|---|---|---|---
|  |  Est | 4.0731 (1.828) | 4.0989 (1.821) | 4.0731 (1.828) |
$n=100$ |  $p=10$ |  Pred | 0.1090 (.0160) |  0.0969 (.0140) | 0.0987 (.0145) |
|  |  Nmse | 0.5119 (.1226) | 0.5162 (.1241) |  0.5118 (.1226) |
|  |  Est | 4.6047 (1.812) |  4.6038 (1.813) | 4.6047 (1.812) |
|  $p=100$ |  Pred | 1.0062 (.0462) | 1.0597 (.0495) |  1.0006 (.0469) |
|  |  Nmse | 0.5801 (.1942) |  0.5800 (.1941) | 0.5801 (.1942) |
|  |  Est | 3.6733 (1.606) | 3.6884 (1.606) | 3.6733 (1.606) |
$n=1000$ |  $p=10$ |  Pred | 0.0098 (.0015) | 0.0098 (.0015) | 0.0098 (.0015) |
|  |  Nmse |  0.4812 (.1271) | 0.4835 (.1260) | 0.4813 (.1271) |
|  |  Est | 3.9972 (1.375) | 3.9986 (1.376) | 3.9972 (1.375) |
|  $p=100$ |  Pred | 0.1000 (.0043) | 0.1032 (.0057) |  0.0999 (.0043) |
|  |  Nmse | 0.5013 (.1061) | 0.5014 (.1063) | 0.5013 (.1062) |
Table 4: Simulation results on simulated data in Model II (approximate low-rank) in inductive matrix completion, with fixed $q=10,k=20$, for different methods, with their standard error in parentheses over 100 replications. ($\kappa$ is the missing rate; Est: average of estimation errors; Pred: average of prediction errors; Nmse: average of normalized estimation errors). |  |  Errors | LMC | MALA | OLS_imp |
---|---|---|---|---|---|---
|  |  Est | 3.8319 (1.691) | 3.8749 (1.719) | 3.8319 (1.690) |
$n=100$ |  $p=10$ |  Pred | 0.1604 (.1271) |  0.1092 (.0153) | 0.1598 (.1322) |
imis 10% |  |  Nmse | 0.5116 (.1154) | 0.5169 (.1147) | 0.5116 (.1155) |
|  |  Est | 5.9500 (2.834) | 5.9452 (2.835) | 5.9500 (2.834) |
|  $p=100$ |  Pred | 4.7640 (5.272) |  4.6964 (5.515) | 4.7658 (5.275) |
|  |  Nmse | 0.7313 (.3454) | 0.7307 (.3455) | 0.7313 (.3454) |
|  |  Est |  4.1838 (1.850) | 4.2535 (1.859) | 4.1839 (1.850) |
$n=100$ |  $p=10$ |  Pred | 0.7221 (.7562) |  0.1498 (.0183) | 0.7371 (.7741) |
imis 30% |  |  Nmse | 0.5182 (.1128) | 0.5283 (.1147) | 0.5182 (.1128) |
|  |  Est | 7.1589 (4.084) | 7.1558 (4.083) | 7.1589 (4.084) |
|  $p=100$ |  Pred |  39.899 (52.40) | 40.233 (51.76) | 39.908 (52.41) |
|  |  Nmse | 0.8998 (.3821) |  0.8994 (.3820) | 0.8998 (.3821) |
|  |  Est | 3.9618 (1.678) | 3.9788 (1.677) | 3.9618 (1.678) |
$n=1000$ |  $p=10$ |  Pred | 0.0409 (.0269) |  0.0110 (.0015) | 0.0409 (.0269) |
imis 10% |  |  Nmse | 0.4968 (.1196) | 0.4989 (.1195) | 0.4968 (.1196) |
|  |  Est | 4.1153 (1.295) | 4.1163 (1.294) | 4.1153 (1.295) |
|  $p=100$ |  Pred | 1.0250 (.9988) |  0.1135 (.0051) | 1.0250 (.9988) |
|  |  Nmse | 0.5060 (.1096) | 0.5062 (.1096) | 0.5060 (.1096) |
|  |  Est | 4.1647 (1.990) | 4.1836 (1.995) | 4.1647 (1.990) |
$n=1000$ |  $p=10$ |  Pred | 0.4615 (.3497) |  0.0141 (.0017) | 0.4616 (.3498) |
imis 30% |  |  Nmse | 0.4905 (.1157) | 0.4933 (.1171) | 0.4905 (.1157) |
|  |  Est | 4.0578 (1.400) |  4.0565 (1.397) | 4.0578 (1.400) |
|  $p=100$ |  Pred | 8.5608 (6.419) |  0.1538 (.0069) | 8.5609 (6.419) |
|  |  Nmse | 0.4944 (.1184) |  0.4943 (.1180) | 0.4944 (.1184) |

##  Appendix A Appendix: proofs

The main technique for our proofs is the oracle type PAC-Bayes bounds, in the
spirit of [7]. We start with a few preliminary lemmas.

###  A.1 Preliminary lemmas

First, we state a version of Bernstein’s inequality from Proposition 2.9 page
24 in [36].

######  Lemma 1 (Bernstein’s inequality).

Let $U_{1}$, …, $U_{n}$ be independent real valued random variables. Let us
assume that there are two constants $v$ and $w$ such that
$\sum_{i=1}^{n}\mathbb{E}[U_{i}^{2}]\leq v$ and for all integers $k\geq 3$,
$\sum_{i=1}^{n}\mathbb{E}\left[(U_{i})^{k}\right]\leq v\frac{k!w^{k-2}}{2}.$
Then, for any $\zeta\in(0,1/w)$,

$\mathbb{E}\exp\left[\zeta\sum_{i=1}^{n}\left[U_{i}-\mathbb{E}(U_{i})\right]\right]\leq\exp\left(\frac{v\zeta^{2}}{2(1-w\zeta)}\right).$

Another basic tool to derive PAC-Bayes bounds is Donsker and Varadhan’s
variational inequality, that we state here, see Lemma 1.1.3 in Catoni [8] for
a proof (among others). From now, for any $\Theta\subset\mathbb{R}^{n_{1}}$ or
$\Theta\subset\mathbb{R}^{n_{1}\times n_{2}}$, we let $\mathcal{P}(\Theta)$
denote the set of all probability distributions on $\Theta$ equipped with the
Borel $\sigma$-algebra. We remind that when
$(\mu,\nu)\in\mathcal{P}(\Theta)^{2}$, the Kullback-Leibler divergence is
defined by $\mathcal{K}(\nu,\mu)=\int\log\left(\frac{{\rm d}\nu}{{\rm
d}\mu}(\theta)\right)\nu({\rm d}\theta)$ if $\nu$ admits a density $\frac{{\rm
d}\nu}{{\rm d}\mu}$ with respect to $\mu$, and $\mathcal{K}(\nu,\mu)=+\infty$
otherwise.

######  Lemma 2 (Donsker and Varadhan’s variational formula).

Let $\mu\in\mathcal{P}(\Theta)$. For any measurable, bounded function
$h:\Theta\rightarrow\mathbb{R}$ we have:

$\log\int{\rm e}^{h(\theta)}\mu({\rm
d}\theta)=\sup_{\rho\in\mathcal{P}(\Theta)}\left[\int h(\theta)\rho({\rm
d}\theta)-KL(\rho\|\mu)\right].$

Moreover, the supremum with respect to $\rho$ in the right-hand side is
reached for the Gibbs measure $\mu_{h}$ defined by its density with respect to
$\pi$

$\frac{{\rm d}\mu_{h}}{{\rm d}\mu}(\theta)=\frac{{\rm e}^{h(\theta)}}{\int{\rm
e}^{h(\vartheta)}\mu({\rm d}\vartheta)}.$ (11)

These two lemmas are the only tools we need to prove Theorem 1 and Theorem 2.
Their proofs are quite similar, with a few differences. For sake of
simplicity, we will state the common parts of the proofs as a separate result
in Lemma 3. Note that the proof of this lemma will use Lemmas 1 and 2.

######  Lemma 3.

Under Assumptions 1 and 2, put

$\alpha=\left(\lambda-\frac{\lambda^{2}C_{1}}{2nq(1-\frac{C_{2}\lambda}{nq})}\right)\quad\text{
and
}\quad\beta=\left(\lambda+\frac{\lambda^{2}C_{1}}{2nq(1-\frac{C_{2}\lambda}{nq})}\right).$
(12)

Then for any $\varepsilon\in(0,1)$, and $\lambda\in(0,nq/C_{2})$,

$\mathbb{E}\int\exp\Biggl{\\{}\alpha\Bigl{(}R(M)-R(M^{*})\Bigr{)}+\lambda\Bigl{(}-r(M)+r(M^{*})\Bigr{)}-\log\left[\frac{d\widehat{\rho}_{\lambda}}{d\pi}(M)\right]-\log\frac{2}{\varepsilon}\Biggr{\\}}\widehat{\rho}_{\lambda}(dM)\leq\frac{\varepsilon}{2}$
(13)

and

$\displaystyle\mathbb{E}\sup_{\rho\in\mathcal{P}(\mathbb{R}^{p\times
k})}\exp\Biggl{[}\beta\left(-\int Rd\rho+R(M^{*})\right)+\lambda\left(\int
rd\rho-r(M^{*})\right)-\mathcal{K}(\rho,\pi)-\log\frac{2}{\varepsilon}\Biggr{]}\leq\frac{\varepsilon}{2}.$
(14)

###### Proof of Lemma 3.

We prove the first inequality (20) as follows. Fix any $M$ with
$\|XMZ\|_{\infty}\leq C$ and put

$T_{ij}=\left(Y_{ij}-(XM^{*}Z)_{ij}\right)^{2}-\left(Y_{ij}-(\Pi_{C}(XMZ))_{ij}\right)^{2}.$

Note that the random variables $T_{ij}$ with $i=1,\ldots,n;j=1,\ldots,q$ are
independent by construction. We have

$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{q}\mathbb{E}[T_{ij}^{2}]$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{q}\mathbb{E}\left[\left(2Y_{ij}-(XM^{*}Z)_{ij}-(XMZ)_{ij}\right)^{2}\left((XM^{*}Z)_{ij}-(XMZ)_{ij}\right)^{2}\right]$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{q}\mathbb{E}\left[\left(2E_{ij}+(XM^{*}Z)_{ij}-(XMZ)_{ij}\right)^{2}\left((XM^{*}Z)_{ij}-(XMZ)_{ij}\right)^{2}\right]$
$\displaystyle\leq\sum_{i=1}^{n}\sum_{j=1}^{q}\mathbb{E}\left[8\left[E_{ij}^{2}+C^{2}\right]\left[(XM^{*}Z)_{ij}-(XMZ)_{ij}\right]^{2}\right]$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{q}8\mathbb{E}\left[E_{ij}^{2}+C^{2}\right]\left[(XM^{*}Z)_{ij}-(XMZ)_{ij}\right]^{2}$
$\displaystyle\leq 8nq(\sigma^{2}+C^{2})\left[R(M)-R(M^{*})\right]$
$\displaystyle=nqC_{1}\left[R(M)-R(M^{*})\right]=:v(M,M^{*}).$

Next we have, for any integer $k\geq 3$, that

$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{q}\mathbb{E}\left[(T_{ij})^{k}\right]\leq$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{q}\mathbb{E}\left[\left|2Y_{ij}-(XM^{*}Z)_{ij}-(XMZ)_{ij}\right|^{k}\left|(XM^{*}Z)_{ij}-(XMZ)_{ij}\right|^{k}\right]$
$\displaystyle\leq$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{q}\mathbb{E}\left[2^{k-1}\left[|2E_{ij}|^{k}+(2C)^{k}\right]\left|(XM^{*}Z)_{ij}-(XMZ)_{ij}\right|^{k}\right]$
$\displaystyle\leq$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{q}\mathbb{E}\left[2^{2k-1}\left(|E_{ij}|^{k}+C^{k}\right)(2C)^{k-2}\left|(XM^{*}Z)_{ij}-(XMZ)_{ij}\right|^{2}\right]$
$\displaystyle\leq$ $\displaystyle
2^{2k-1}\left[\sigma^{2}k!\xi^{k-2}+C^{k}\right](2C)^{k-2}\sum_{i=1}^{n}\sum_{j=1}^{q}\mathbb{E}\left|(XM^{*}Z)_{ij}-(XMZ)_{ij}\right|^{2}$
$\displaystyle\leq$
$\displaystyle\frac{2^{3k-3}\left[\sigma^{2}k!\xi^{k-2}+C^{k}\right]C^{k-2}}{8(\sigma^{2}+C^{2})}v(M,M^{*})$
$\displaystyle\leq$
$\displaystyle\frac{2^{3k-6}\left[\sigma^{2}\xi^{k-2}+C^{k}\right]C^{k-2}}{(\sigma^{2}+C^{2})}k!v(M,M^{*})$
$\displaystyle\leq$ $\displaystyle
2^{3k-5}\left[\xi^{k-2}+C^{k-2}\right]C^{k-2}k!v(M,M^{*})$ $\displaystyle\leq$
$\displaystyle
2^{3k-4}\max(\xi,C)^{k-2}C^{k-2}k!v(M,M^{*})=[2^{3}\max(\xi,C)C]^{k-2}2^{2}k!v(M,M^{*})$

and use the fact that, for any $k\geq 3$, $2^{2}\leq 2^{3(k-2)}/2$ to obtain

$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{q}\mathbb{E}\left[(T_{ij})^{k}\right]\leq\frac{[2^{6}\max(\xi,C)C]^{k-2}k!v(M,M^{*})}{2}=v(M,M^{*})\frac{k!C_{2}^{k-2}}{2}.$

Thus, we can apply Lemma 1 with $U_{i}:=T_{i}$, $v:=v(M,M^{*})$, $w:=C_{2}$
and $\zeta:=\lambda/nq$. We obtain, for any $\lambda\in(0,nq/w)=(0,nq/C_{2})$,

$\displaystyle\mathbb{E}\exp\left[\lambda\Bigl{(}R(M)-R(M^{*})-r(M)+r(M^{*})\Bigr{)}\right]\leq\exp\left[\frac{v\lambda^{2}}{2(nq)^{2}(1-\frac{w\lambda}{nq})}\right]=\exp\left[\frac{C_{1}\left[R(M)-R(M^{*})\right]\lambda^{2}}{2nq(1-\frac{C_{2}\lambda}{nq})}\right].$

Rearranging terms, and using the definition of $\alpha$ in (12),

$\mathbb{E}\exp\left[\alpha\Bigl{(}R(M)-R(M^{*})\Bigr{)}+\lambda\Bigl{(}-r(M)+r(M^{*})\Bigr{)}\right]\leq
1.$

Multiplying both sides by $\varepsilon/2$ and then integrating w.r.t. the
probability distribution $\pi(.)$, we get

$\int\mathbb{E}\exp\Biggl{[}\alpha\Bigl{(}R(M)-R(M^{*})\Bigr{)}+\lambda\Bigl{(}-r(M)+r(M^{*})\Bigr{)}-\log\frac{2}{\varepsilon}\Biggr{]}\pi(dM)\leq\frac{\varepsilon}{2}.$

Next, Fubini’s theorem gives

$\mathbb{E}\int\exp\Biggl{[}\alpha\Bigl{(}R(M)-R(M^{*})\Bigr{)}+\lambda\Bigl{(}-r(M)+r(M^{*})\Bigr{)}-\log\frac{2}{\varepsilon}\Biggr{]}\pi(dM)\leq\frac{\varepsilon}{2}.$

and note that for any measurable function $h$,

$\displaystyle\int\exp[h(M)]\pi(dM)=\int\exp\left[h(M)-\log\frac{{\rm
d}\widehat{\rho}_{\lambda}}{{\rm d}\pi}(M)\right]\widehat{\rho}_{\lambda}(dM)$

to get (13).

Let us now prove (14). Here again, we start with an application of Lemma 1,
but this time with $U_{i}:=-T_{i}$ (we keep $v:=v(M,M^{*})$, $w:=C_{2}$ and
$\zeta:=\lambda/nq$). We obtain, for any $\lambda\in(0,nq/C_{2})$,

$\mathbb{E}\exp\left[\lambda\Bigl{(}r(M)+r(M^{*})-R(M)+R(M^{*})\Bigr{)}\right]\leq\exp\left[\frac{C_{1}\left[R(M)-R(M^{*})\right]\lambda^{2}}{2nq(1-\frac{C_{2}\lambda}{nq})}\right].$

Rearranging terms, using the definition of $\beta$ in (12) and multiplying
both sides by $\varepsilon/2$, we obtain

$\mathbb{E}\exp\Biggl{[}\beta\left(-R(M)+R(M^{*})\right)+\lambda\left(r(M)-r(M^{*})\right)-\log\frac{2}{\varepsilon}\Biggr{]}\leq\frac{\varepsilon}{2}.$

We integrate with respect to $\pi$ and use Fubini to get:

$\mathbb{E}\int\exp\Biggl{[}\beta\left(-R(M)+R(M^{*})\right)+\lambda\left(r(M)-r(M^{*})\right)-\log\frac{2}{\varepsilon}\Biggr{]}\pi({\rm
d}M)\leq\frac{\varepsilon}{2}.$

Here, we use a different argument from the proof of the first inequality: we
use Lemma 2 on the integral, this gives directly (14). ∎

Finally, in both proofs, we will use quite often distributions
$\rho\in\mathcal{P}(\mathbb{R}^{p\times k})$ that will be defined as
translations of the prior $\pi$. We introduce the following notation.

######  Definition 1.

For any matrix $\bar{M}\in\mathbb{R}^{p\times k}$, we define
$\rho_{\bar{M}}\in\mathcal{P}(\mathbb{R}^{p\times k})$ by

$\rho_{\bar{M}}(M)=\pi(\bar{M}-M).$

The following technical lemmas from [12] will be useful in the proofs.

######  Lemma 4 (Lemma 1 in [12]).

We have $\int\|M\|_{F}^{2}\pi({\rm d}M)\leq pk\tau^{2}.$

######  Lemma 5 (Lemma 2 in [12]).

For any $\bar{M}\in\mathbb{R}^{p\times k}$, we have

$\mathcal{K}(\rho_{\bar{M}},\pi)\leq 2{\rm
rank}(\bar{M})(k+p+2)\log\left(1+\frac{\|\bar{M}\|_{F}}{\tau\sqrt{2{\rm
rank}(\bar{M})}}\right)$

with the convention $0\log(1+0/0)=0$.

###  A.2 Proof of Theorem 1

###### Proof of Theorem 1.

An application of Jensen’s inequality on inequality (13) yields

$\mathbb{E}\exp\Biggl{[}\alpha\left(\int
Rd\widehat{\rho}_{\lambda}-R(M^{*})\right)+\lambda\left(-\int
rd\widehat{\rho}_{\lambda}+r(M^{*})\right)-\mathcal{K}(\widehat{\rho}_{\lambda},\pi)-\log\frac{2}{\varepsilon}\Biggr{]}\leq\frac{\varepsilon}{2}.$

Using the standard Chernoff’s trick to transform an exponential moment
inequality into a deviation inequality, that is:
$\exp(x)\geq\mathbf{1}_{\mathbb{R}_{+}}(x)$, we obtain

$\displaystyle\mathbb{P}\Biggl{\\{}\Biggr{[}\alpha\left(\int
Rd\widehat{\rho}_{\lambda}-R(M^{*})\right)+\lambda\left(-\int
rd\widehat{\rho}_{\lambda}+r(M^{*})\right)-\mathcal{K}(\widehat{\rho}_{\lambda},\pi)-\log\frac{2}{\varepsilon}\Biggr{]}\geq
0\Biggr{\\}}\leq\frac{\varepsilon}{2}$ (15)

Using (2) we have

$\displaystyle\int
Rd\widehat{\rho}_{\lambda}-R(M^{*})=\dfrac{1}{nq}\int\left\|\Pi_{C}(XMZ)-XM^{*}Z\right\|_{F}^{2}\widehat{\rho}_{\lambda}({\rm
d}M)$
$\displaystyle\geq\dfrac{1}{nq}\left\|\int\Pi_{C}(XMZ)\widehat{\rho}_{\lambda}({\rm
d}M)-XM^{*}Z\right\|_{F}^{2}$
$\displaystyle\geq\dfrac{1}{nq}\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F}^{2}$

where we used Jensen’s inequality in the first line, and the definition of
$\widehat{XMZ}_{\lambda}$ from the first to the second line. Plugging this
into our probability bound (15), and dividing both sides by $\alpha$, we
obtain

$\mathbb{P}\Biggl{\\{}\dfrac{1}{nq}\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F}^{2}\leq\frac{\int
rd\widehat{\rho}_{\lambda}-r(M^{*})+\frac{1}{\lambda}\left[\mathcal{K}(\widehat{\rho}_{\lambda},\pi)+\log\frac{2}{\varepsilon}\right]}{\frac{\alpha}{\lambda}}\Biggr{\\}}\geq
1-\frac{\varepsilon}{2}$

under the additional condition that $\lambda$ is such that $\alpha>0$, that we
will assume from now (note that this is satisfied by $\lambda^{*}$). Using
Lemma 2 we can rewrite this as

$\mathbb{P}\Biggl{\\{}\dfrac{1}{nq}\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F}^{2}\leq\inf_{\rho\in\mathcal{P}(\mathbb{R}^{p\times
k})}\frac{\int
rd\rho-r(M^{*})+\frac{1}{\lambda}\left[\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon}\right]}{\frac{\alpha}{\lambda}}\Biggr{\\}}\geq
1-\frac{\varepsilon}{2}.$ (16)

We consider now the consequences of the second inequality in Lemma 3, that is
(14). Chernoff’s trick and rearranging terms a little, we get

$\displaystyle\mathbb{P}\Biggl{\\{}\forall\rho\in\mathcal{P}(\mathbb{R}^{p\times
k}),\int rd\rho-r(M^{*})\leq\frac{\beta}{\lambda}\left[\int
Rd\rho-R(M^{*})\right]+\frac{1}{\lambda}\left[\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon}\right]\Biggr{\\}}\geq
1-\frac{\varepsilon}{2}.$

which we can rewrite as, $\forall\rho\in\mathcal{P}(\mathbb{R}^{p\times k}),$
with probability at least $1-\frac{\varepsilon}{2}$,

$\displaystyle\int
rd\rho-r(M^{*})\leq\frac{\beta}{\lambda}\int\dfrac{1}{nq}\|\Pi_{C}(XMZ)-XM^{*}Z\|_{F}^{2}\rho({\rm
d}M)+\frac{1}{\lambda}\left[\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon}\right].$
(17)

Combining (17) and (16) with a union bound argument gives the following bound,
with probability at least $1-\varepsilon$,

$\displaystyle\frac{1}{nq}\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F}^{2}\leq\inf_{\rho\in\mathcal{P}(\mathbb{R}^{p\times
k})}\frac{\beta\int\frac{1}{nq}\|\Pi_{C}(XMZ)-XM^{*}Z\|_{F}^{2}\rho({\rm
d}M)+2\left[\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon}\right]}{\alpha}.$

Noting that, for any $(i,j)$, $(XM^{*}Z)_{i,j}\in[-C,C]$ implies that
$|(\Pi_{C}(XMZ))_{i,j}-(XM^{*}Z)_{i,j}|\leq|(XMZ)_{i,j}-(XM^{*}Z)_{i,j}|$ and
thus

$\displaystyle\frac{1}{nq}\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F}^{2}\leq\inf_{\rho\in\mathcal{P}(\mathbb{R}^{p\times
k})}\frac{\beta\int\frac{1}{nq}\|XMZ-XM^{*}Z\|_{F}^{2}\rho({\rm
d}M)+2\left[\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon}\right]}{\alpha}.$

The end of the proof consists in making the right-hand side in the inequality
more explicit. In order to do so, we restrict the infimum bound above to the
distributions given by Definition 1.

$\mathbb{P}\Biggl{\\{}\frac{1}{nq}\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F}^{2}\leq\inf_{\bar{M}\in\mathbb{R}^{p\times
k}}\frac{\beta\int\frac{1}{nq}\|XMZ-XM^{*}Z\|_{F}^{2}\rho_{\bar{M}}({\rm
d}M)+2\left[\mathcal{K}(\rho_{\bar{M}},\pi)+\log\frac{2}{\varepsilon}\right]}{\alpha}\Biggr{\\}}\geq
1-\varepsilon.$ (18)

We see immediately that Dalalyan’s lemma will be extremely useful for that.
First, Lemma 5 provides an upper bound on $\mathcal{K}(\rho_{\bar{M}},\pi)$.
Moreover,

$\displaystyle\int\|XMZ-XM^{*}Z\|_{F}^{2}\rho_{\bar{M}}({\rm d}M)$
$\displaystyle\leq\int\|X\bar{M}Z-XM^{*}Z-XMZ\|_{F}^{2}\pi({\rm d}M)$
$\displaystyle=\|X\bar{M}Z-XM^{*}Z\|_{F}^{2}-2\int\sum_{i,j}(X\bar{M}Z-XM^{*}Z)_{j,i}(XMZ)_{i,j}\pi({\rm
d}M)+\int\|XMZ\|_{F}^{2}\pi({\rm d}M).$

The second term in the right-hand side is null because $\pi$ is centered, and
thus

$\displaystyle\int\|XMZ-XM^{*}Z\|_{F}^{2}\rho_{\bar{M}}({\rm d}M)$
$\displaystyle\leq\|X\bar{M}Z-XM^{*}Z\|_{F}^{2}+\int\|XMZ\|_{F}^{2}\,\pi({\rm
d}M)$
$\displaystyle\leq\|X\bar{M}Z-XM^{*}Z\|_{F}^{2}+\|X\|_{F}^{2}\|Z\|_{F}^{2}\int\|M\|_{F}^{2}\,\pi({\rm
d}M)$
$\displaystyle\leq\|X\bar{M}Z-XM^{*}Z\|_{F}^{2}+\|X\|_{F}^{2}\|Z\|_{F}^{2}pk\tau^{2}$

where we used elementary properties of the Frobenius norm, and Lemma 4 in the
last line. We can now plug this (and Lemma 5) back into (18) to get:

$\mathbb{P}\Biggl{\\{}\frac{1}{nq}\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F}^{2}\leq\inf_{\bar{M}\in\mathbb{R}^{p\times
k}}\Biggl{[}\frac{\beta}{\alpha}\frac{1}{nq}\|X\bar{M}Z-XM^{*}Z\|_{F}^{2}+\frac{\beta}{\alpha}\frac{1}{nq}\|X\|_{F}^{2}\|Z\|_{F}^{2}pk\tau^{2}\\\
+\frac{1}{\alpha}\left(4{\rm
rank}(\bar{M})(k+p+2)\log\left(1+\frac{\|\bar{M}\|_{F}}{\tau\sqrt{2{\rm
rank}(\bar{M})}}\right)+2\log\frac{2}{\varepsilon}\right)\Biggr{]}\Biggr{\\}}\geq
1-\varepsilon.$

We are now making the constants be explicit. First, if $\lambda\leq
nq/(2C_{2})$, then $2nq(1-C_{2}\lambda/nq)\geq np$ and thus

$\frac{\beta}{\alpha}=\frac{1+\frac{\lambda
C_{1}}{2nq(1-\frac{C_{2}\lambda}{nq})}}{1-\frac{\lambda
C_{1}}{2nq(1-\frac{C_{2}\lambda}{nq})}}\leq\frac{1+\frac{\lambda
C_{1}}{nq}}{1-\frac{\lambda C_{1}}{nq}}.$

Then, $\lambda\leq\frac{nq\delta}{C_{1}(1+\delta)}$ leads to
$\frac{\beta}{\alpha}\leq(1+\delta).$

Note that $\lambda^{*}=nq\min(1/(2C_{2}),\delta/[C_{1}(1+\delta)])$ satisfies
these two conditions, so from now $\lambda=\lambda^{*}$. We also use the
following:

$\displaystyle\frac{1}{\alpha}=\frac{1}{\lambda^{*}\left(1-\frac{\lambda^{*}C_{1}}{2nq(1-C_{2}\lambda^{*}/nq)}\right)}\leq\frac{\beta}{\lambda^{*}\alpha}\leq\frac{(1+\delta)}{nq\min(1/(2C_{2}),\delta/[C_{1}(1+\delta)])}\leq\frac{C_{1}(1+\delta)^{2}}{nq\delta}.$

So far the bound is:

$\mathbb{P}\Biggl{\\{}\frac{1}{nq}\left\|\widehat{XMZ}_{\lambda^{*}}-XM^{*}Z\right\|_{F}^{2}\leq\inf_{\bar{M}\in\mathbb{R}^{p\times
k}}\Biggl{[}\frac{(1+\delta)}{nq}\|X\bar{M}Z-XM^{*}Z\|_{F}^{2}+\frac{(1+\delta)}{nq}\|X\|_{F}^{2}\|Z\|_{F}^{2}pk\tau^{2}+\\\
\frac{C_{1}(1+\delta)^{2}\left(4{\rm
rank}(\bar{M})(k+p+2)\log\left(1+\frac{\|\bar{M}\|_{F}}{\tau\sqrt{2{\rm
rank}(\bar{M})}}\right)+2\log\frac{2}{\varepsilon}\right)}{nq\delta}\Biggr{]}\Biggr{\\}}\geq
1-\varepsilon.$

In particular, with probability at least $1-\varepsilon$, the choice
$\tau^{2}=C_{1}(k+p)/(nkq\|X\|_{F}^{2}\|Z\|_{F}^{2})$ gives

$\frac{1}{nq}\left\|\widehat{XMZ}_{\lambda^{*}}-XM^{*}Z\right\|_{F}^{2}\leq\inf_{\bar{M}\in\mathbb{R}^{p\times
k}}\Biggl{[}\frac{(1+\delta)}{nq}\|X\bar{M}Z-XM^{*}Z\|_{F}^{2}+\frac{C_{1}(1+\delta)(k+p)}{nq}+\\\
\frac{C_{1}(1+\delta)^{2}\left(4{\rm
rank}(\bar{M})(k+p+2)\log\left(1+\frac{\|X\|_{F}\|Z\|_{F}\|\bar{M}\|_{F}}{\sqrt{C_{1}}}\sqrt{\frac{nkq}{(k+p){\rm
rank}(\bar{M})}}\right)+2\log\frac{2}{\varepsilon}\right)}{nq\delta}\Biggr{]}.$

∎

###  A.3 Proof of theorem 2

###### Proof of Theorem 2.

We also start with an application of Lemma 3, and focus on (13), applied to
$\varepsilon:=\varepsilon_{n}$, that is:

$\displaystyle\mathbb{E}\int\exp\Biggl{\\{}\alpha\Bigl{(}R(M)-R(M^{*})\Bigr{)}+\lambda\Bigl{(}-r(M)+r(M^{*})\Bigr{)}-\log\left[\frac{d\widehat{\rho}_{\lambda}}{d\pi}(M)\right]-\log\frac{2}{\varepsilon_{n}}\Biggr{\\}}\widehat{\rho}_{\lambda}(dM)\leq\frac{\varepsilon_{n}}{2}.$

Using Chernoff’s trick, this gives:

$\mathbb{E}\Bigl{[}\mathbb{P}_{M\sim\widehat{\rho}_{\lambda}}(M\in\mathcal{A}_{n})\Bigr{]}\geq
1-\frac{\varepsilon_{n}}{2}$

where

$\mathcal{A}_{n}=\left\\{M:\alpha\Bigl{(}R(M)-R(M^{*})\Bigr{)}+\lambda\Bigl{(}-r(M)+r(M^{*})\Bigr{)}\leq\log\left[\frac{d\widehat{\rho}_{\lambda}}{d\pi}(M)\right]+\log\frac{2}{\varepsilon_{n}}\right\\}.$

Using the definition of $\widehat{\rho}_{\lambda}$, for $M\in\mathcal{A}_{n}$
we have

$\displaystyle\alpha\Bigl{(}R(M)-R(M^{*})\Bigr{)}$
$\displaystyle\leq\lambda\Bigl{(}r(M)-r(M^{*})\Bigr{)}+\log\left[\frac{d\widehat{\rho}_{\lambda}}{d\pi}(M)\right]+\log\frac{2}{\varepsilon_{n}}$
$\displaystyle\leq-\log\int\exp\left[-\lambda r(M)\right]\pi({\rm d}M)-\lambda
r(M^{*})+\log\frac{2}{\varepsilon_{n}}$ $\displaystyle=\lambda\Bigl{(}\int
r(M)\widehat{\rho}_{\lambda}({\rm
d}M)-r(M^{*})\Bigr{)}+\mathcal{K}(\widehat{\rho}_{\lambda},\pi)+\log\frac{2}{\varepsilon_{n}}$
$\displaystyle=\inf_{\rho}\left\\{\lambda\Bigl{(}\int r(M)\rho({\rm
d}M)-r(M^{*})\Bigr{)}+\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon_{n}}\right\\}.$

Now, let us define

$\mathcal{B}_{n}=\left\\{\forall\rho\text{, }\beta\left(-\int
Rd\rho+R(M^{*})\right)+\lambda\left(\int
rd\rho-r(M^{*})\right)\leq\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon_{n}}\right\\}.$

Using (14), we have that

$\mathbb{E}\Bigl{[}\mathbf{1}_{\mathcal{B}_{n}}\Bigr{]}\geq
1-\frac{\varepsilon_{n}}{2}.$

We will now prove that, if $\lambda$ is such that $\alpha>0$,

$\mathbb{E}\Bigl{[}\mathbb{P}_{M\sim\widehat{\rho}_{\lambda}}(M\in\mathcal{M}_{n})\Bigr{]}\geq\mathbb{E}\Bigl{[}\mathbb{P}_{M\sim\widehat{\rho}_{\lambda}}(M\in\mathcal{A}_{n})\mathbf{1}_{\mathcal{B}_{n}}\Bigr{]}$

which, together with

$\displaystyle\mathbb{E}\Bigl{[}\mathbb{P}_{M\sim\widehat{\rho}_{\lambda}}(M\in\mathcal{A}_{n})\mathbf{1}_{\mathcal{B}_{n}}\Bigr{]}$
$\displaystyle=\mathbb{E}\Bigl{[}(1-\mathbb{P}_{M\sim\widehat{\rho}_{\lambda}}(M\notin\mathcal{A}_{n}))(1-\mathbf{1}_{\mathcal{B}^{c}_{n}})\Bigr{]}$
$\displaystyle\geq\mathbb{E}\Bigl{[}1-\mathbb{P}_{M\sim\widehat{\rho}_{\lambda}}(M\notin\mathcal{A}_{n})-\mathbf{1}_{\mathcal{B}^{c}_{n}}\Bigr{]}$
$\displaystyle\geq 1-\varepsilon_{n}$

will bring

$\mathbb{E}\Bigl{[}\mathbb{P}_{M\sim\widehat{\rho}_{\lambda}}(M\in\mathcal{M}_{n})\Bigr{]}\geq
1-\varepsilon_{n}.$

In order to do so, assume that we are on the set $\mathcal{B}_{n}$, and let
$M\in\mathcal{A}_{n}$. Then,

$\displaystyle\alpha\Bigl{(}R(M)-R(M^{*})\Bigr{)}$
$\displaystyle\leq\inf_{\rho}\left\\{\lambda\Bigl{(}\int r(M)\rho({\rm
d}M)-r(M^{*})\Bigr{)}+\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon_{n}}\right\\}$
$\displaystyle\leq\inf_{\rho}\left\\{\beta\Bigl{(}\int R(M)\rho({\rm
d}M)-R(M^{*})\Bigr{)}+2\mathcal{K}(\rho,\pi)+2\log\frac{2}{\varepsilon_{n}}\right\\}$

that is,

$R(M)-R(M^{*})\leq\inf_{\rho\in\mathcal{P}(\mathbb{R}^{p\times
k})}\frac{\beta\left[\int
Rd\rho-R(M^{*})\right]+2\left[\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon}\right]}{\alpha}$

or, rewriting it in terms of norms,

$\left\|\Pi_{C}(XMZ)-XM^{*}Z\right\|_{F}^{2}\\\
\leq\inf_{\bar{M}\in\mathbb{R}^{p\times k}}\frac{\beta\int\|XMZ-
XM^{*}Z\|_{F}^{2}\rho_{\bar{M}}({\rm
d}M)+2\left[\mathcal{K}(\rho_{\bar{M}},\pi)+\log\frac{2}{\varepsilon}\right]}{\alpha}.$

We upper-bound the right-hand side exactly as in the proof of Theorem 1, this
gives $M\in\mathcal{M}_{n}$.

∎

###  A.4 Proof of Theorem 3

######  Lemma 6.

Under Assumptions 1 and 3, put

$\alpha^{\prime}=\left(\lambda-\frac{\lambda^{2}C^{\prime}_{1}}{2m(1-\frac{C^{\prime}_{2}\lambda}{m})}\right)\quad\text{
and
}\quad\beta^{\prime}=\left(\lambda+\frac{\lambda^{2}C^{\prime}_{1}}{2m(1-\frac{C^{\prime}_{2}\lambda}{m})}\right).$
(19)

Then for any $\varepsilon\in(0,1)$, and $\lambda\in(0,m/C^{\prime}_{2})$,

$\mathbb{E}\int\exp\Biggl{\\{}\alpha^{\prime}\Bigl{(}R^{\prime}(M)-R^{\prime}(M^{*})\Bigr{)}+\lambda\Bigl{(}-r^{\prime}(M)+r^{\prime}(M^{*})\Bigr{)}-\log\left[\frac{d\widehat{\rho}^{\prime}_{\lambda}}{d\pi}(M)\right]-\log\frac{2}{\varepsilon}\Biggr{\\}}\widehat{\rho}^{\prime}_{\lambda}(dM)\leq\frac{\varepsilon}{2}$
(20)

and

$\displaystyle\mathbb{E}\sup_{\rho\in\mathcal{P}(\mathbb{R}^{p\times
k})}\exp\Biggl{[}\beta^{\prime}\left(-\int
R^{\prime}d\rho+R^{\prime}(M^{*})\right)+\lambda\left(\int
r^{\prime}d\rho-r^{\prime}(M^{*})\right)-\mathcal{K}(\rho,\pi)-\log\frac{2}{\varepsilon}\Biggr{]}\leq\frac{\varepsilon}{2}.$
(21)

###### Proof of Lemma 6.

The inequality (20) is proved in a similar way to the proof of Lemma 3. That
is we apply Lemma 1 to the following independent random variables

$V_{i}=\left(Y_{i}-(XM^{*}Z)_{i}\right)^{2}-\left(Y_{i}-(\Pi_{C}(XMZ))_{i}\right)^{2},i=1,\ldots,m.$

The proof of the inequality (21) is processed similar in the proof of Lemma 3
in which we apply Lemma 1 to the independent random variables
$-V_{i},i=1,\ldots,m.$ ∎

###### Proof of Theorem 3.

Similar to the proof of Theorem 1, until the (16), and noting that using (6)
we have

$\displaystyle\int
R^{\prime}d\widehat{\rho}_{\lambda}-R^{\prime}(M^{*})=\int\left\|\Pi_{C}(XMZ)-XM^{*}Z\right\|_{F,\Pi}^{2}\widehat{\rho}^{\prime}_{\lambda}({\rm
d}M)$
$\displaystyle\geq\left\|\int\Pi_{C}(XMZ)\widehat{\rho}^{\prime}_{\lambda}({\rm
d}M)-XM^{*}Z\right\|_{F,\Pi}^{2}$
$\displaystyle\geq\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F,\Pi}^{2},$

thus we obtain

$\mathbb{P}\Biggl{\\{}\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F,\Pi}^{2}\leq\inf_{\rho\in\mathcal{P}(\mathbb{R}^{p\times
k})}\frac{\int
r^{\prime}d\rho-r^{\prime}(M^{*})+\frac{1}{\lambda}\left[\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon}\right]}{\frac{\alpha^{\prime}}{\lambda}}\Biggr{\\}}\geq
1-\frac{\varepsilon}{2}.$ (22)

We consider now the consequences of inequality (21) in Lemma 6. Chernoff’s
trick and rearranging terms a little, we get,
$\forall\rho\in\mathcal{P}(\mathbb{R}^{p\times k})$, with probability at least
$1-\frac{\varepsilon}{2}$,

$\displaystyle\int
r^{\prime}d\rho-r^{\prime}(M^{*})\leq\frac{\beta^{\prime}}{\lambda}\int\|\Pi_{C}(XMZ)-XM^{*}Z\|_{F,\Pi}^{2}\rho({\rm
d}M)+\frac{1}{\lambda}\left[\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon}\right].$
(23)

Combining (23) and (22) with a union bound argument gives the bound and noting
that for any $(i,j)$, $(XM^{*}Z)_{i,j}\in[-C,C]$ implies that
$|(\Pi_{C}(XMZ))_{i,j}-(XM^{*}Z)_{i,j}|\leq|(XMZ)_{i,j}-(XM^{*}Z)_{i,j}|$ and
thus

$\displaystyle\mathbb{P}\Biggl{\\{}\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F,\Pi}^{2}\leq\inf_{\rho\in\mathcal{P}(\mathbb{R}^{p\times
k})}\frac{\beta^{\prime}\int\|XMZ-XM^{*}Z\|_{F,\Pi}^{2}\rho({\rm
d}M)+2\left[\mathcal{K}(\rho,\pi)+\log\frac{2}{\varepsilon}\right]}{\alpha^{\prime}}\Biggr{\\}}\geq
1-\varepsilon.$

We are now making the right-hand side in the inequality more explicit. In
order to do so, we restrict the infimum bound above to the distributions given
by Definition 1.

$\mathbb{P}\Biggl{\\{}\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F,\Pi}^{2}\leq\inf_{\bar{M}\in\mathbb{R}^{p\times
k}}\frac{\beta^{\prime}\int\|XMZ-XM^{*}Z\|_{F,\Pi}^{2}\rho_{\bar{M}}({\rm
d}M)+2\left[\mathcal{K}(\rho_{\bar{M}},\pi)+\log\frac{2}{\varepsilon}\right]}{\alpha^{\prime}}\Biggr{\\}}\geq
1-\varepsilon.$ (24)

We see immediately that Dalalyan’s lemma will be extremely useful for that.
First, Lemma 5 provides an upper bound on $\mathcal{K}(\rho_{\bar{M}},\pi)$.

Moreover,

$\displaystyle\int\|XMZ-XM^{*}Z\|_{F,\Pi}^{2}\rho_{\bar{M}}({\rm d}M)$
$\displaystyle\leq\int\|X\bar{M}Z-XM^{*}Z-XMZ\|_{F,\Pi}^{2}\pi({\rm d}M)$
$\displaystyle=\|X\bar{M}Z-XM^{*}Z\|_{F,\Pi}^{2}-2\int\sum_{i,j}\Pi_{i,j}(X\bar{M}Z-XM^{*}Z)_{j,i}(XMZ)_{i,j}\pi({\rm
d}M)+\int\|XMZ\|_{F,\Pi}^{2}\pi({\rm d}M).$

The second term in the above right-hand side is null because $\pi$ is
centered, and thus

$\displaystyle\int\|XMZ-XM^{*}Z\|_{F,\Pi}^{2}\rho_{\bar{M}}({\rm d}M)$
$\displaystyle\leq\|X\bar{M}Z-XM^{*}Z\|_{F,\Pi}^{2}+\int\|XMZ\|_{F,\Pi}^{2}\pi({\rm
d}M)$
$\displaystyle\leq\|X\bar{M}Z-XM^{*}Z\|_{F,\Pi}^{2}+\int\|XMZ\|_{F}^{2}\,\pi({\rm
d}M)$
$\displaystyle\leq\|X\bar{M}Z-XM^{*}Z\|_{F,\Pi}^{2}+\|X\|_{F}^{2}\|Z\|_{F}^{2}\int\|M\|_{F}^{2}\,\pi({\rm
d}M)$
$\displaystyle\leq\|X\bar{M}Z-XM^{*}Z\|_{F,\Pi}^{2}+\|X\|_{F}^{2}\|Z\|_{F}^{2}pk\tau^{2}$

where we used elementary properties of the Frobenius norm, and Lemma 4 in the
last line. We can now plug this (and Lemma 5) back into (24) to get:

$\mathbb{P}\Biggl{\\{}\left\|\widehat{XMZ}_{\lambda}-XM^{*}Z\right\|_{F,\Pi}^{2}\leq\inf_{\bar{M}\in\mathbb{R}^{p\times
k}}\Biggl{[}\frac{\beta^{\prime}}{\alpha^{\prime}}\|X\bar{M}Z-XM^{*}Z\|_{F,\Pi}^{2}+\frac{\beta^{\prime}}{\alpha^{\prime}}\|X\|_{F}^{2}\|Z\|_{F}^{2}pk\tau^{2}\\\
+\frac{1}{\alpha^{\prime}}\left(4{\rm
rank}(\bar{M})(k+p+2)\log\left(1+\frac{\|\bar{M}\|_{F}}{\tau\sqrt{2{\rm
rank}(\bar{M})}}\right)+2\log\frac{2}{\varepsilon}\right)\Biggr{]}\Biggr{\\}}\geq
1-\varepsilon.$

We are now making the constants be explicit. First, if $\lambda\leq
m/(2C^{\prime}_{2})$, then $2m(1-C^{\prime}_{2}\lambda/m)\geq m$ and thus

$\frac{\beta^{\prime}}{\alpha^{\prime}}=\frac{1+\frac{\lambda
C^{\prime}_{1}}{2m(1-\frac{C^{\prime}_{2}\lambda}{m})}}{1-\frac{\lambda
C^{\prime}_{1}}{2m(1-\frac{C^{\prime}_{2}\lambda}{m})}}\leq\frac{1+\frac{\lambda
C^{\prime}_{1}}{m}}{1-\frac{\lambda C^{\prime}_{1}}{m}}.$

Then, $\lambda\leq\frac{m\delta}{C^{\prime}_{1}(1+\delta)}$ leads to
$\frac{\beta^{\prime}}{\alpha^{\prime}}\leq(1+\delta).$

Note that
$\lambda^{\prime*}=m\min(1/(2C^{\prime}_{2}),\delta/[C^{\prime}_{1}(1+\delta)])$
satisfies these two conditions, so from now $\lambda=\lambda^{\prime*}$. We
also use the following:

$\displaystyle\frac{1}{\alpha^{\prime}}=\frac{1}{\lambda^{\prime*}\left(1-\frac{\lambda^{\prime*}C^{\prime}_{1}}{2m(1-C_{2}\lambda^{\prime*}/m)}\right)}\leq\frac{\beta^{\prime}}{\lambda^{\prime*}\alpha^{\prime}}\leq\frac{(1+\delta)}{m\min(1/(2C^{\prime}_{2}),\delta/[C^{\prime}_{1}(1+\delta)])}\leq\frac{C^{\prime}_{1}(1+\delta)^{2}}{m\delta}.$

So far the bound is:

$\mathbb{P}\Biggl{\\{}\left\|\widehat{XMZ}_{\lambda^{\prime*}}-XM^{*}Z\right\|_{F,\Pi}^{2}\leq\inf_{\bar{M}\in\mathbb{R}^{p\times
k}}\Biggl{[}(1+\delta)\|X\bar{M}Z-XM^{*}Z\|_{F,\Pi}^{2}+(1+\delta)\|X\|_{F}^{2}\|Z\|_{F}^{2}pk\tau^{2}\\\
+\frac{C^{\prime}_{1}(1+\delta)^{2}\left(4{\rm
rank}(\bar{M})(k+p+2)\log\left(1+\frac{\|\bar{M}\|_{F}}{\tau\sqrt{2{\rm
rank}(\bar{M})}}\right)+2\log\frac{2}{\varepsilon}\right)}{m\delta}\Biggr{]}\Biggr{\\}}\geq
1-\varepsilon.$

In particular, with probability at least $1-\varepsilon$, the choice
$\tau^{2}=C^{\prime}_{1}(k+p)/(mkp\|X\|_{F}^{2}\|Z\|_{F}^{2})$ gives

$\left\|\widehat{XMZ}_{\lambda^{\prime*}}-XM^{*}Z\right\|_{F,\Pi}^{2}\leq\inf_{\bar{M}\in\mathbb{R}^{p\times
k}}\Biggl{[}(1+\delta)\|X\bar{M}Z-XM^{*}Z\|_{F,\Pi}^{2}+\frac{C^{\prime}_{1}(1+\delta)(k+p)}{m}+\\\
\frac{C^{\prime}_{1}(1+\delta)^{2}\left(4{\rm
rank}(\bar{M})(k+p+2)\log\left(1+\frac{\|X\|_{F}\|Z\|_{F}\|\bar{M}\|_{F}}{\sqrt{C^{\prime}_{1}}}\sqrt{\frac{mkp}{(k+p){\rm
rank}(\bar{M})}}\right)+2\log\frac{2}{\varepsilon}\right)}{m\delta}\Biggr{]}.$

∎

###  A.5 Proof of theorem 4

###### Proof.

The proof is proceeded completely similar to the proof of Theorem 2, in
Section A.3.

∎

##  Appendix B Comments on algorithm implementation

For the case of inductive matrix completion, we write the logarithm of the
density of the posterior

$\log\widehat{\rho}_{\lambda}(M)=-\frac{\lambda}{n}\sum_{i=1}^{n}(Y_{i}-(\Pi_{C}(XMZ))_{i})^{2}-\frac{p+m+2}{2}\log\det(\tau^{2}\mathbf{I}_{m}+MM^{\intercal}).$

Let us now differentiate this expression in $M$. Note that the term
$(Y_{i}-(\Pi_{C}(XMZ))_{i})^{2}$ does actually not depend on $M$ locally if
$(XMZ)_{i}\notin[-C,C]$, in this case its differential with respect to $M$ is
$\mathbf{0}_{p\times m}$. Otherwise,
$(Y_{i}-(\Pi_{C}(XMZ))_{i})^{2}=(Y_{i}-(XMZ)_{i})^{2}$. In order to be able to
differentiate the term $(XMZ)_{i}$, let us introduce a notation for the
entries of $\mathcal{I}_{i}$: $\mathcal{I}_{i}=(a_{i},b_{i})$. Then
$\nabla(XMZ)_{i}=D$ where the matrix $D\in\mathbb{R}^{p\times m}$ satisfies
$D_{x,y}=\mathbf{1}_{\\{x=b_{j}\\}}X_{a_{j},y}$. Then

$\nabla\log\widehat{\rho}_{\lambda}(M)=\frac{2\lambda}{n}\sum_{i=1}^{n}(\nabla(XMZ)_{i})(Y_{i}-(XMZ)_{i})\mathbf{1}_{\\{|(XMZ)_{i}|<C\\}}-(p+m+2)(\tau^{2}\mathbf{I}_{m}+MM^{\intercal})^{-1}M.$

In the above calculation still requires to calculate a $p\times p$ matrix
inversion at each iteration, for very large $p$, this might be expensive and
can slow down the algorithm. Therefore, we could replace this matrix inversion
by its accurately approximation through a convex optimization. It is noted
that the matrix $\mathbf{B}:=(\tau^{2}\mathbf{I}_{m}+MM^{\intercal})^{-1}M$ is
the solution to the following convex optimization problem:
$\min_{\mathbf{B}}\big{\\{}\|\mathbf{I}_{p}-M^{\top}\mathbf{B}\|_{F}^{2}+\tau^{2}\|\mathbf{B}\|_{F}^{2}\big{\\}}.$
The solution of this optimization problem can be conveniently obtained by
using the package ‘glmnet’ [18] (with the family option ‘mgaussian’). This
avoids to perform matrix inversion or other costly calculation. However, we
note here that the LMC algorithm is being used with approximate gradient
evaluation, theoretical assessment of this approach can be found in [14].

Algorithm 1 LMC

Input: The data.

Parameters: Positive real numbers $\lambda,\tau,h,T$.

Output: The matrix $\widehat{M}$

Initialize: $M_{0},\widehat{M}=\mathbf{0}_{m\times p}$

for $k\leftarrow 1$ to $T$ do

Sample $M_{k}$ from (8);

$\widehat{M}\leftarrow\widehat{M}+M_{k}/T$

end for

Algorithm 2 MALA

Input: The data.

Parameters: Positive real numbers $\lambda,\tau,h,T$

Output: The matrix $\widehat{M}$

Initialize: $M_{0};\widehat{M}=\mathbf{0}_{m\times p}$

for $k=1$ to $T$ do

Sample $\tilde{M}_{k}$ from (9).

Set $M_{k}=\tilde{M}_{k}$ with probability $A_{MALA}$, from (10), otherwise
$M_{k}=M_{k-1}$ .

$\widehat{M}\leftarrow\widehat{M}+M_{k}/T$ .

end for

## References

  * [1] P. Alquier.  User-friendly introduction to PAC-Bayes bounds.  preprint arXiv:2110.11216, 2021.
  * [2] P. Alquier, J. Ridgway, and N. Chopin.  On the properties of variational approximations of Gibbs posteriors.  The Journal of Machine Learning Research, 17(1):8374–8414, 2016\.
  * [3] T. W. Anderson.  Estimating linear restrictions on regression coefficients for multivariate normal distributions.  Annals of mathematical statistics, 22(3):327–351, 1951.
  * [4] P. G. Bissiri, C. C. Holmes, and S. G. Walker.  A general framework for updating belief distributions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), pages n/a–n/a, 2016.
  * [5] E. J. Candès and Y. Plan.  Matrix completion with noise.  Proceedings of the IEEE, 98(6):925–936, 2010.
  * [6] E. J. Candes, M. B. Wakin, and S. P. Boyd.  Enhancing sparsity by reweighted $\ell_{1}$ minimization.  Journal of Fourier analysis and applications, 14(5-6):877–905, 2008\.
  * [7] O. Catoni.  Statistical learning theory and stochastic optimization, volume 1851 of Saint-Flour Summer School on Probability Theory 2001 (Jean Picard ed.), Lecture Notes in Mathematics.  Springer-Verlag, Berlin, 2004.
  * [8] O. Catoni.  PAC-Bayesian supervised classification: the thermodynamics of statistical learning.  IMS Lecture Notes—Monograph Series, 56. Institute of Mathematical Statistics, Beachwood, OH, 2007.
  * [9] V. Cottet and P. Alquier.  1-Bit matrix completion: PAC-Bayesian analysis of a variational approximation.  Machine Learning, 107(3):579–603, 2018.
  * [10] A. Dalalyan and A. B. Tsybakov.  Aggregation by exponential weighting, sharp PAC-Bayesian bounds and sparsity.  Machine Learning, 72(1-2):39–61, 2008.
  * [11] A. S. Dalalyan.  Theoretical guarantees for approximate sampling from smooth and log-concave densities.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 3(79):651–676, 2017.
  * [12] A. S. Dalalyan.  Exponential weights in multivariate regression and a low-rankness favoring prior.  In Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, volume 56, pages 1465–1483. Institut Henri Poincaré, 2020\.
  * [13] A. S. Dalalyan, E. Grappin, and Q. Paris.  On the exponentially weighted aggregate with the laplace prior.  The Annals of Statistics, 46(5):2452–2478, 2018.
  * [14] A. S. Dalalyan and A. Karagulyan.  User-friendly guarantees for the langevin monte carlo with inaccurate gradient.  Stochastic Processes and their Applications, 129(12):5278–5311, 2019.
  * [15] A. S. Dalalyan and A. B. Tsybakov.  Sparse regression learning by aggregation and langevin monte-carlo.  Journal of Computer and System Sciences, 78(5):1423–1443, 2012\.
  * [16] A. Durmus and E. Moulines.  High-dimensional Bayesian inference via the unadjusted langevin algorithm.  Bernoulli, 25(4A):2854–2882, 2019.
  * [17] R. Foygel, O. Shamir, N. Srebro, and R. Salakhutdinov.  Learning with the weighted trace-norm under arbitrary sampling distributions.  In Advances in Neural Information Processing Systems, pages 2133–2141, 2011.
  * [18] J. Friedman, T. Hastie, and R. Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  Journal of Statistical Software, 33(1):1–22, 2010.
  * [19] P. Grünwald and T. Van Ommen.  Inconsistency of Bayesian inference for misspecified linear models, and a proposal for repairing it.  Bayesian Analysis, 12(4):1069–1103, 2017.
  * [20] B. Guedj.  A primer on PAC-Bayesian learning.  In Proceedings of the second congress of the French Mathematical Society, 2019.
  * [21] T. Hastie and R. Mazumder.  softImpute: Matrix Completion via Iterative Soft-Thresholded SVD, 2021.  R package version 1.4-1.
  * [22] W. Hoeffding.  Probability inequalities for sums of bounded random variables.  Journal of the American Statistical Association, 58(301):13–30, 1963.
  * [23] C.-J. Hsieh, N. Natarajan, and I. Dhillon.  Pu learning for matrix completion.  In International conference on machine learning, pages 2445–2453. PMLR, 2015.
  * [24] A. J. Izenman.  Modern multivariate statistical techniques.  Regression, classification and manifold learning, 10:978–0, 2008\.
  * [25] P. Jain and I. S. Dhillon.  Provable inductive matrix completion.  arXiv preprint arXiv:1306.0626, 2013.
  * [26] S. Jana, N. Balakrishnan, and J. S. Hamid.  Bayesian growth curve model useful for high-dimensional longitudinal data.  Journal of Applied Statistics, 46(5):814–834, 2019.
  * [27] S. Jana et al.  Inference for Generalized Multivariate Analysis of Variance (GMANOVA) Models and High-dimensional Extensions.  PhD thesis, 2017.
  * [28] O. Klopp.  Noisy low-rank matrix completion with general sampling distribution.  Bernoulli, 20(1):282–303, 2014.
  * [29] J. Knoblauch, J. Jewson, and T. Damoulas.  An optimization-centric view on bayes’ rule: Reviewing and generalizing variational inference.  Journal of Machine Learning Research, 23(132):1–109, 2022.
  * [30] V. Koltchinskii, K. Lounici, and A. B. Tsybakov.  Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion.  Ann. Statist., 39(5):2302–2329, 2011.
  * [31] A. Kshirsagar and W. Smith.  Growth curves, volume 145.  CRC Press, 1995.
  * [32] C. Luo, J. Liang, G. Li, F. Wang, C. Zhang, D. K. Dey, and K. Chen.  Leveraging mixed and incomplete outcomes via reduced-rank modeling.  Journal of Multivariate Analysis, 167:378–394, 2018.
  * [33] T. T. Mai and P. Alquier.  A Bayesian approach for noisy matrix completion: Optimal rate under general sampling distribution.  Electron. J. Statist., 9(1):823–841, 2015.
  * [34] T. T. Mai and P. Alquier.  Pseudo-Bayesian quantum tomography with rank-adaptation.  Journal of Statistical Planning and Inference, 184:62–76, 2017\.
  * [35] T. T. Mai and P. Alquier.  Optimal quasi-bayesian reduced rank regression with incomplete response.  arXiv preprint arXiv:2206.08619, 2022.
  * [36] P. Massart.  Concentration inequalities and model selection, volume 1896 of Lecture Notes in Mathematics.  Springer, Berlin, 2007.  Lectures from the 33rd Summer School on Probability Theory held in Saint-Flour, July 6–23, 2003, Edited by Jean Picard.
  * [37] D. McAllester.  Some PAC-Bayesian theorems.  In Proceedings of the Eleventh Annual Conference on Computational Learning Theory, pages 230–234, New York, 1998. ACM.
  * [38] N. Natarajan and I. S. Dhillon.  Inductive matrix completion for predicting gene–disease associations.  Bioinformatics, 30(12):i60–i68, 2014.
  * [39] S. Negahban and M. J. Wainwright.  Restricted strong convexity and weighted matrix completion: optimal bounds with noise.  J. Mach. Learn. Res., 13:1665–1697, 2012.
  * [40] R. F. Potthoff and S. Roy.  A generalized multivariate analysis of variance model useful especially for growth curve problems.  Biometrika, 51(3-4):313–326, 1964.
  * [41] R Core Team.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria, 2022.
  * [42] P. Rigollet and A. B. Tsybakov.  Sparse estimation by exponential weighting.  Statistical Science, 27(4):558–575, 2012.
  * [43] G. O. Roberts and J. S. Rosenthal.  Optimal scaling of discrete approximations to langevin diffusions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1):255–268, 1998.
  * [44] G. O. Roberts and O. Stramer.  Langevin diffusions and metropolis-hastings algorithms.  Methodology and computing in applied probability, 4(4):337–357, 2002.
  * [45] D. v. Rosen.  Bilinear regression with rank restrictions on the mean and dispersion matrix.  In Methodology and Applications of Statistics, pages 193–211. Springer, 2021.
  * [46] J. Shawe-Taylor and R. Williamson.  A PAC analysis of a Bayes estimator.  In Proceedings of the Tenth Annual Conference on Computational Learning Theory, pages 2–9, New York, 1997. ACM.
  * [47] D. Von Rosen.  Bilinear regression analysis.  Lecture notes in statistics, 220, 2018.
  * [48] R. F. Woolson and J. D. Leeper.  Growth curve analysis of complete and incomplete longitudinal data.  Communications in Statistics-theory and Methods, 9(14):1491–1513, 1980.
  * [49] L. Yang, J. Fang, H. Duan, H. Li, and B. Zeng.  Fast low-rank Bayesian matrix completion with hierarchical gaussian prior models.  IEEE Transactions on Signal Processing, 66(11):2804–2817, 2018\.
  * [50] W. Zhang, H. Xu, X. Li, Q. Gao, and L. Wang.  Drimc: an improved drug repositioning approach using bayesian inductive matrix completion.  Bioinformatics, 36(9):2839–2847, 2020.
  * [51] P. Zilber and B. Nadler.  Inductive matrix completion: No bad local minima and a fast algorithm.  In K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu, and S. Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 27671–27692. PMLR, 2022.