From Bilinear Regression to Inductive Matrix Completion: A Quasi-Bayesian Analysis

Abstract

In this paper, we study the problem of bilinear regression, a type of statistical modeling that deals with multiple variables and multiple responses. One of the main difficulties that arise in this problem is the presence of missing data in the response matrix, a problem known as inductive matrix completion. To address these issues, we propose a novel approach that combines elements of Bayesian statistics with a quasi-likelihood method. Our proposed method starts by addressing the problem of bilinear regression using a quasi-Bayesian approach. The quasi-likelihood method that we employ in this step allows us to handle the complex relationships between the variables in a more robust way. Next, we adapt our approach to the context of inductive matrix completion. We make use of a low-rankness assumption and leverage the powerful PAC-Bayes bound technique to provide statistical properties for our proposed estimators and for the quasi-posteriors. To compute the estimators, we propose a Langevin Monte Carlo method to obtain approximate solutions to the problem of inductive matrix completion in a computationally efficient manner. To demonstrate the effectiveness of our proposed methods, we conduct a series of numerical studies. These studies allow us to evaluate the performance of our estimators under different conditions and provide a clear illustration of the strengths and limitations of our approach.

1. Introduction

In this paper, we investigate the bilinear regression model, a statistical method that assumes a linear relationship between a set of multiple response variables and two sets of covariates. This model, also known as the growth curve model or generalized multivariate analysis model, is commonly used for analyzing longitudinal data, as shown in previous studies such as [1,2,3,4,5,6]. However, these studies only cover the scenario in which the response matrix is fully observed.

Recently, the bilinear regression model with incomplete response has been introduced and studied as the so-called inductive matrix completion, which is a generalization of the matrix completion problem [7,8]. This problem has attracted significant attention in various fields, such as drug repositioning [9], collaborative filtering [10], and genomics [7]. Inductive matrix completion is a challenging problem that arises when some of the entries in the response matrix are missing, which makes it difficult to infer the underlying relationship between the variables.

In this work, we explore the problem of bilinear regression and inductive matrix completion under a low-rank constraint on the coefficient matrix. Most existing approaches for these problems are frequentist methods, such as maximum likelihood estimation [2] or penalized optimization [7]. These methods are effective in providing point estimates for the parameters of the model but lack the ability to provide a full probabilistic characterization of the uncertainty. Recently, Bayesian approaches have been considered for these problems. For example, the paper [11] proposed a Bayesian approach for bilinear regression, and a Bayesian method was proposed for inductive matrix completion in the work [9]. However, unlike frequentist approaches, the statistical properties of the Bayesian approach for these models have not been fully explored yet.

The aim of this paper is to address an existing gap in the understanding of the bilinear regression and inductive matrix completion problems. To achieve this goal, we propose a novel approach that combines elements of Bayesian statistics with a quasi-likelihood method. Specifically, we start by addressing the problem of bilinear regression using a quasi-Bayesian approach, where a quasi-likelihood is employed. We then generalize this approach to the problem of inductive matrix completion. To ensure that our method is adaptive to the rank of the coefficient matrix, we use a spectral scaled Student prior distribution, which allows us to prove that the posterior mean satisfies a tight oracle inequality. This result demonstrates that our method is able to accurately estimate the parameters of the model, even when the rank of the coefficient matrix is unknown. Additionally, we also prove the contraction properties of the posteriors, which further enhances the performance of our method.

The proposed method in this paper, the quasi-Bayesian approach, is an extension of the traditional Bayesian approach and is becoming increasingly popular in statistics and machine learning as a technique for generalized Bayesian inference, as noted in studies such as [12,13,14]. This approach allows for more flexibility in the modeling assumptions by replacing the likelihood function with a more general notion of risk or quasi-likelihood.

To provide theoretical guarantees for our proposed quasi-posteriors, we make use of the PAC-Bayesian technique [15,16,17]. This technique provides bounds on the generalization error of a learned estimator, and has been widely used in the literature as described in recent reviews and introductions, such as [18,19]. The PAC-Bayes bounds have been successfully applied in the context of matrix estimation problems as shown in studies such as [20,21,22,23]. Interestingly, by using the PAC-Bayesian technique for inductive matrix completion, we do not need to make any assumptions about the distribution of the missing entries in the response matrix. This is in contrast with previous works on matrix completion, such as [24,25,26,27,28,29], which typically require assumptions about the missing data. This makes our method more versatile and applicable to a wider range of problems.

The proposed method in this paper makes use of a spectral scaled Student prior, which is a specific choice of prior distribution. This choice is motivated by recent works in which it has been shown to lead to optimal rates in a variety of problems, including high-dimensional regression [30], image denoising [31] and reduced rank regression [23]. Although this prior is not conjugate to our problems, it allows for the convenient implementation of gradient-based sampling methods, which makes it computationally efficient.

To compute the proposed estimators and sample from the quasi-posterior, we employ a Langevin Monte Carlo (LMC) method. This method is a widely used algorithm for approximating complex distributions and allows for efficient computation of the proposed estimators. The LMC method allows us to obtain approximate solutions to the problem of bilinear regression and inductive matrix completion in a computationally efficient manner.

Furthermore, we use numerical studies to demonstrate the effectiveness of our proposed methods. These studies enable us to evaluate the performance of our estimators in various scenarios and provide insight into the capabilities and limitations of our approach. By conducting a thorough evaluation, we can gain a better understanding of how our method performs in different settings and make any necessary adjustments to improve its performance. We also compare our method with the ordinary least squared method. This comparison allows us to demonstrate the superiority of our method over the traditional approach in certain scenarios, such as when the response matrix contains missing data. Overall, the numerical studies serve as a valuable tool for assessing the effectiveness of our proposed method and provide a clear illustration of its strengths and weaknesses, as well as its improvement over the traditional method.

The remainder of this paper is organized as follows: In Section 2, we present the problem of bilinear regression, and introduce the low-rank promoting prior distribution that we use to address this problem. In Section 3, we extend our approach to the problem of inductive matrix completion. In Section 4, we discuss the Langevin Monte Carlo method used for the computation of the estimators, and present numerical studies to demonstrate the effectiveness of our proposed methods. Finally, we conclude our work and provide a summary of our findings in Section 5. The technical proofs are provided in Appendix A for the interested readers.

Notation 1.

Let ${\mathbb{R}}^{n_1\times n_2}$ denote the set of $n_1\times n_2$ matrices with real elements. Let $A^{⊺}\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 i-th row and j-th column elements 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 $\mathrm{Tr}\left(B\right)$ 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 ${\parallel A\parallel }_{\infty }={max}_{i,j}\left|A_{i,j}\right|$; its Frobenius norm ${\parallel A\parallel }_F$ is defined by ${\parallel A\parallel }_F^2=\mathrm{Tr}\left(A^{⊺}A\right)={\sum }_{i,j}{A}_{i,j}^2$ and $\mathrm{rank}\left(A\right)$ 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}$ be a known within-individuals design matrix. Consider the bilinear regression model as follows:

$$Y=X{M}^{*}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}\left(E\right)=0.$ The main assumption here is the low-rank restriction on the model parameter that $\mathrm{rank}\left(M^{*}\right)<min(p,k)$.

The model presented in Equation (1) is a bilinear regression model, which can be seen as a generalization of the reduced rank regression problem. In the case where $k=q$ and $Z={\mathbf{I}}_q$, the model simplifies to the traditional reduced rank regression problem, which has been well-studied in the literature, such as [32,33]. However, in this paper, we consider the more general case where the matrix Z contains additional explanatory variables.

The low-rank assumption in this model can be interpreted as indicating the presence of a latent process that affects the response data, not only through the “between-individuals” structure of the model but also through the “within-individuals” structure. This model is often referred to as the growth curve model or the generalized multivariate analysis model (GMANOVA) and has been studied in depth in the literature, such as [2].

Assumption 1.

There is a known constant $C<+\infty$ such that $\parallel X{M}^{*}{Z\parallel }_{\infty }\le C.$

From Assumption 1, it is not reliable to return predictions $XMZ$ with entries that are outside of interval $[-C,C]$. However, for computational reasons, 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\left(A\right)=arg\underset{{\parallel B\parallel }_{\infty }\le C}{min}{\parallel A-B\parallel }_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\left(M\right)$ the empirical risk of M as

$$r\left(M\right)={\displaystyle \frac{1}{nq}}{\parallel Y-{\Pi }_C\left(XMZ\right)\parallel }_F^2$$

and its expectation is denoted by

$$R\left(M\right)=\mathbb{E}\left[r\left(M\right)\right]=\mathbb{E}\left[{\left(Y_{11}-{\left({\Pi }_C\left(XMZ\right)\right)}_{{}_{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\left(M\right)-R\left(M^{*}\right)$ is small. Under the assumption that $E_{ij}$ has a finite variance, using the Pythagorean theorem, we have

$$R\left(M\right)-R\left(M^{*}\right)={\displaystyle \frac{1}{nq}}{\parallel {\Pi }_C\left(XMZ\right)-X{M}^{*}Z\parallel }_F^2$$

(2)

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

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

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

It is worth noting that for a specific choice of $\lambda =nq/\left(2{\sigma }^2\right)$, the posterior distribution obtained corresponds to the case where the noise term $E_{ij}$ is assumed to be Gaussian distributed with a mean of 0 and a variance of ${\sigma }^2$. However, our theoretical results hold under a more general class of noise distributions. It is known that a small enough $\lambda$ is sufficient when the model is misspecified [14]. Additionally, even in the case of Gaussian noise, in high-dimensional settings, a smaller value of $\lambda$ than $n/\left(2{\sigma }^2\right)$ leads to better adaptation properties [30,34]. The precise choice of $\lambda$ in our method will be further discussed below.

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

$${\widehat{XMZ}}_{\lambda }=\int {\Pi }_C\left(XMZ\right){\widehat{\rho }}_{\lambda }\left(dM\right).$$

(3)

It is worth noting that from the simulation experiments, it is observed that using reasonable values for C, the Monte Carlo algorithm never samples matrices M such that ${\Pi }_C\left(XMZ\right)\ne XMZ$. In other words, the boundedness constraint has very little impact on practice, and it is mainly necessary for technical proofs. If one is interested in obtaining an estimator of $M^{*}$ instead of an estimator of $X{M}^{*}Z$, when $X^{⊺}X$ and $Z{Z}^{⊺}$ are invertible, one can consider the estimator ${\widehat{M}}_{\lambda }={\left(X^{⊺}X\right)}^{-1}{X}^{⊺}{\widehat{XMZ}}_{\lambda }{Z}^{⊺}{\left(Z{Z}^{⊺}\right)}^{-1}$ and note that $X{\widehat{M}}_{\lambda }Z={\widehat{XMZ}}_{\lambda }$. This estimator can be used to obtain the estimator of $M^{*}$ in case the inverses of $X^{⊺}X$ and $Z{Z}^{⊺}$ are computationally feasible.

The quasi-posterior distribution investigated in this paper is often referred to as the “Gibbs posterior” in the PAC-Bayes approach. This terminology is used in the literature such as [17,18,19,35,36]. Additionally, the estimator ${\widehat{M}}_{\lambda }$ is sometimes referred to as the Gibbs estimator or the exponentially weighted aggregate (EWA) in the literature, such as [34,37,38].

2.2. Prior Specification

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

$$\begin{array}{c}\hfill \pi \left(M\right)\propto det{({\tau }^2{\mathbf{I}}_m+M{M}^{⊺})}^{-(p+m+2)/2}.\end{array}$$

(4)

This prior distribution is designed to promote low-rankness by placing more probability mass on matrices with smaller singular values, which leads to sparse solutions. This prior has a similar form to the one used in related works such as [39,40], and it is known to lead to good performance in different problems, such as high-dimensional regression and image denoising. It can be verified that

$$\pi \left(M\right)\propto \prod _{j=1}^m{({\tau }^2+s_j{\left(M\right)}^2)}^{-(p+m+2)/2},$$

where $s_j\left(M\right)$ denotes the jth largest singular value of M. The above expression is a scaled Student distribution evaluated at $s_j\left(M\right)$, which can be seen as a way to approximate sparsity on $s_j\left(M\right)$ [30]. The log-sum function ${\sum }_{j=1}^{m}log({\tau }^2+s_j{\left(M\right)}^2)$ used by [39,40] is also known to enforce approximate sparsity on the singular values of $s_j\left(M\right)$. This means that under this prior, most of the $s_j\left(M\right)$ are close to 0, which implies that M is well approximated by a low-rank matrix. Therefore, it has the ability to promote the low-rankness of M.

As previously stated, this prior is not conjugate to the problem at hand. However, it is particularly convenient to implement gradient-based sampling algorithms, such as the Langevin Monte Carlo method, which will be discussed in more detail in Section 4. This is because the gradient of the log-posterior can be computed efficiently, and it allows for an efficient implementation of the LMC algorithm.

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\ge 2,\mathbb{E}\left(\right|E_{i,j}{|}^k)\le {\sigma }^{2}k!{\xi }^{k-2}/2.$$

Let us put

$$C_1=8({\sigma }^2+C^2);C_2=64Cmax(\xi ,C);{\tau }^{*}=\sqrt{C_1(k+p)/{(nkq\parallel X\parallel }_F^2{\parallel Z\parallel }_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 }^{*}:=nqmin(1/\left(2C_2\right),\delta /\left[C_1(1+\delta )\right])$. Then, for any $\epsilon \in (0,1)$, we have, with probability at least $1-\epsilon$ on the sample,

$$\begin{array}{c}{∥{\widehat{XMZ}}_{{\lambda }^{*}}-X{M}^{*}Z∥}_F^2\le \underset{0\le r\le pk}{inf}\underset{\begin{array}{c}\hfill \begin{array}{c}\overline{M}\in {\mathbb{R}}^{p\times k}\\ \mathrm{rank}\left(\overline{M}\right)\le r\end{array}\end{array}}{inf}\{(1+\delta ){∥X\overline{M}Z-X{M}^{*}Z∥}_F^2+\hfill \\ \hfill \frac{C_1{(1+\delta )}^2}{\delta }\left[4r(k+p+2)log\left(1+\frac{{\parallel X\parallel }_F{\parallel Z\parallel }_F{\parallel \overline{M}\parallel }_F}{\sqrt{C_1}}\sqrt{\frac{nkq}{r(k+p)}}\right)+k+p+2log\frac{2}{\epsilon }\right]\}.\end{array}$$

The choice of $\lambda ={\lambda }^{*}$ is determined by optimizing an upper bound on the risk R (as shown in the proof of this theorem). However, it is important to note that this choice may not necessarily be the best choice in practice, even though it gives a good estimate of the order of magnitude for $\lambda$. To ensure optimal performance, the user can use cross-validation to properly adjust the temperature parameter. Additionally, it is worth noting that $\mathrm{rank}\left(\overline{M}\right)\ne 0$ is not a requirement in the above formula. If $\mathrm{rank}\left(\overline{M}\right)=0$, then $\overline{M}=0$ and we interpret $0log(1+0/0)$ as 0.

The proof of this theorem is based on the PAC-Bayes theory, and it is provided in the Appendix A. By taking $\overline{M}=M^{*}$, we can obtain an upper bound on the infimum, leading to the following result.

Corollary 1.

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

$$\begin{array}{c}{∥{\widehat{XMZ}}_{{\lambda }^{*}}-X{M}^{*}Z∥}_F^2\le \hfill \\ \hfill \frac{C_1{(1+\delta )}^2}{\delta }\left[4r^{*}(k+p+2)log\left(1+\frac{{\parallel X\parallel }_F{\parallel Z\parallel }_F{\parallel M^{*}\parallel }_F}{\sqrt{C_1}}\sqrt{\frac{nkq}{r(k+p)}}\right)+k+p+2log\frac{2}{\epsilon }\right].\end{array}$$

Theorem 1 provides an understanding of the statistical properties of the posterior mean. However, it is also important to understand the contraction properties of the quasi-posterior distribution. In the following theorem, we aim to provide a result that demonstrates this aspect of the proposed method.

Theorem 2.

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

$$\begin{array}{c}{\mathcal{M}}_n=\{M\in {\mathbb{R}}^{p\times k}:{∥{\widehat{XMZ}}_{{\lambda }^{*}}-X{M}^{*}Z∥}_F^2\le \underset{0\le r\le pk}{inf}\underset{\begin{array}{c}\hfill \begin{array}{c}\overline{M}\in {\mathbb{R}}^{p\times k}\\ \mathrm{rank}\left(\overline{M}\right)\le r\end{array}\end{array}}{inf}\{(1+\delta ){∥X\overline{M}Z-X{M}^{*}Z∥}_F^2+\hfill \\ \hfill \frac{C_1{(1+\delta )}^2}{\delta }\left[4r(k+p+2)log\left(1+\frac{{\parallel X\parallel }_F{\parallel Z\parallel }_F{\parallel \overline{M}\parallel }_F}{\sqrt{C_1}}\sqrt{\frac{nkq}{r(k+p)}}\right)+k+p+2log\frac{2}{{\epsilon }_n}\right]\}.\end{array}$$

Then,

$$\mathbb{E}\left[{\mathbb{P}}_{M\sim {\widehat{\rho }}_{\lambda }}(M\in {\mathcal{M}}_n)\right]\ge 1-{\epsilon }_n\underset{n\to \infty }{\to }1.$$

The proof of this theorem is provided in 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 independent and identically random pairs $({\mathcal{I}}_1,Y_1),\dots ,({\mathcal{I}}_m,Y_m)$ given by

$$Y_i={\left(X{M}^{*}Z\right)}_{{\mathcal{I}}_i}+{\mathcal{E}}_i,i=1,\dots ,m$$

(5)

where $M^{*}\in {\mathbb{R}}^{p\times k}$ is the unknown parameter matrix expected to be low-rank and observation sample size is assumed that $m<nq$. The noise variables ${\mathcal{E}}_i$ are assumed to be independent with $\mathbb{E}\left({\mathcal{E}}_i\right)=0.$ The variables ${\mathcal{I}}_i$ are independent and identical copies of a random variable $\mathcal{I}$ having distribution $\Pi$ on the set $\{1,\dots ,n\}\times \{1,\dots ,q\}$, we denote ${\Pi }_{x,y}:=\Pi (\mathcal{I}=(x,y))$.

Our goal in this paper is to investigate the problem of bilinear regression and also address the case where the response matrix contains missing data, a problem known as inductive matrix completion. In particular, when $p=n,k=q$ and $X={\mathbf{I}}_n,Z={\mathbf{I}}_q$ are the identity matrices, the problem reduces to the traditional matrix completion problem, which has been well studied in the literature [26]. Similarly, when $k=q$ and $Z={\mathbf{I}}_q$ is the identity matrix, the problem becomes the reduced rank regression problem with incomplete response, which has also been studied in recent works, such as [23,41]. However, in the context of inductive matrix completion, we focus on the more general case where X and Z contain additional explanatory variables; in other words, we consider the side information from the n users and the q items in our model [8].

It has been acknowledged that there are two different ways to model the observed values of Y, either by including or excluding the possibility of observing the same entry multiple times. Previous studies have examined both of these methods, such as the examination of matrix completion without replacement in [25] and with replacement in [26]. Both methods have practical uses and use similar techniques for estimation. This particular study focuses on the scenario where the variables ${\mathcal{I}}_i$ are independently and identically distributed, meaning that it is possible to observe the same entry multiple times. Additionally, it is important to note that, according to the findings presented in Section 6 of [42], the results of this study can also be applied to the scenario of sampling without replacement, as long as the sampling is performed uniformly and there is no observation noise.

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 ${\parallel A\parallel }_{F,P}^2={\sum }_{i,j}P\left[(i,j)\right]A_{i,j}^2$; note that when P is the uniform distribution, then ${\parallel A\parallel }_{F,P}^2={\parallel A\parallel }_F^2/\left(n_1{n}_2\right)$.

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

$$r^{\prime }\left(M\right)=\frac{1}{m}\sum _{i=1}^m{\left(Y_i-{\left({\Pi }_C\left(XMZ\right)\right)}_{{\mathcal{I}}_i}\right)}^2,$$

$$R^{\prime }\left(M\right)=\mathbb{E}\left[r^{\prime }\left(M\right)\right]=\mathbb{E}\left[{\left(Y_1-{\left({\Pi }_C\left(XMZ\right)\right)}_{{\mathcal{I}}_1}\right)}^2\right].$$

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 }\left(M\right)-R^{\prime }\left(M^{*}\right)$ is small. Under the assumption that ${\mathcal{E}}_i$ has a finite variance, based on the Pythagorean theorem, we have

$$R^{\prime }\left(M\right)-R^{\prime }\left(M^{*}\right)={\parallel {\Pi }_C\left(XMZ\right)-X{M}^{*}Z\parallel }_{F,\Pi }^2$$

(6)

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

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

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

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

The truncated posterior mean of $XMZ$ is given by

$${\widehat{XMZ}}_{\lambda }=\int {\Pi }_C\left(XM\right){\widehat{\rho }}_{\lambda }^{\prime }\left(dM\right).$$

(7)

Here, for the same technical reasons as in the context of bilinear regression, 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 ${\widehat{\rho }}_{\lambda }^{\prime }$ 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,\dots ,{\mathcal{E}}_m$ are independent of ${\mathcal{I}}_1,\dots ,{\mathcal{I}}_m$. There exist two known constants ${\sigma }^{\prime }>0$ and ${\xi }^{\prime }>0$ such that

$$\forall k\ge 2,\mathbb{E}\left(\right|{\mathcal{E}}_i{|}^k)\le {\sigma }^{\prime 2}k!{\xi }^{\prime k-2}/2.$$

Assumptions 1 and 3 are both standard; they have been used in [41] for theoretical analysis of reduced rank regression and in [26] for trace regression and matrix completion.

Let us put

$$C_1^{\prime }=8({\sigma }^{\prime 2}+C^2);C_2^{\prime }=64Cmax({\xi }^{\prime },C);{\tau }^{*}=\sqrt{C_1^{\prime }(k+p)/{(mkp\parallel X\parallel }_F^2{\parallel Z\parallel }_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 *}:=mmin(1/\left(2C_2^{\prime }\right),\delta /\left[C_1^{\prime }(1+\delta )\right])$. Then, for any $\epsilon \in (0,1)$, we have, with probability at least $1-\epsilon$ on the sample,

$$\begin{array}{c}{∥{\widehat{XMZ}}_{{\lambda }^{\prime *}}-X{M}^{*}Z∥}_{F,\Pi }^2\le \underset{0\le r\le pk}{inf}\underset{\begin{array}{c}\hfill \begin{array}{c}\overline{M}\in {\mathbb{R}}^{p\times k}\\ \mathrm{rank}\left(\overline{M}\right)\le r\end{array}\end{array}}{inf}\{(1+\delta ){\parallel X\overline{M}Z-X{M}^{*}Z\parallel }_{F,\Pi }^2+\hfill \\ \hfill \frac{C_1^{\prime }{(1+\delta )}^2}{\delta }\frac{\left(4r(k+p+2)log\left(1+\frac{{\parallel X\parallel }_F{\parallel Z\parallel }_F{\parallel \overline{M}\parallel }_F}{\sqrt{C_1}}\sqrt{\frac{mkp}{r(k+p)}}\right)+k+p+2log\frac{2}{\epsilon }\right)}{m}\}.\end{array}$$

Similar to the context of bilinear regression, the choices of $\lambda ={\lambda }^{*},\tau ={\tau }^{*}$ come from the optimization of an 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 a good order of magnitude for tuning these parameters. The user could use cross-validation to properly tune them in practice. Note again that $\mathrm{rank}\left(\overline{M}\right)\ne 0$ is not required in the above formula, if $\mathrm{rank}\left(\overline{M}\right)=0$ then $\overline{M}=0$ and we interpret $0log(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 $\overline{M}$ by taking $\overline{M}=M^{*}$, which leads to the following result.

Corollary 2.

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

$$R_{\delta ,m,p,k,r^{*},\epsilon }:=\frac{C_1^{\prime }{(1+\delta )}^2}{\delta }\frac{\left(4r(k+p+2)log\left(1+\frac{{\parallel X\parallel }_F{\parallel Z\parallel }_F{\parallel \overline{M}\parallel }_F}{\sqrt{C_1}}\sqrt{\frac{mkp}{r(k+p)}}\right)+k+p+2log\frac{2}{\epsilon }\right)}{m},$$

then

$$\begin{array}{c}\hfill {∥{\widehat{XMZ}}_{{\lambda }^{\prime *}}-X{M}^{*}Z∥}_{F,\Pi }^2\le R_{\delta ,m,p,k,r^{*},\epsilon }\end{array}$$

and in particular, if the sampling distribution Π is uniform,

$$\begin{array}{c}\hfill {\displaystyle \frac{\parallel {\widehat{XMZ}}_{{\lambda }^{\prime *}}-X{M}^{*}{Z\parallel }_F^2}{nq}}\le R_{\delta ,m,p,k,r^{*},\epsilon }.\end{array}$$

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 [8].

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

Theorem 4.

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

$$\begin{array}{c}{\Omega }_m=\{M\in {\mathbb{R}}^{p\times k}:{\parallel {\Pi }_C\left(XMZ\right)-X{M}^{*}Z\parallel }_{F,\Pi }^2\le \hfill \\ \underset{1\le r\le pk}{inf}\underset{\begin{array}{c}\hfill \begin{array}{c}\overline{M}\in {\mathbb{R}}^{p\times k}\\ \mathrm{rank}\left(\overline{M}\right)\le r\end{array}\end{array}}{inf}[(1+\delta ){\parallel X\overline{M}Z-X{M}^{*}Z\parallel }_{F,\Pi }^2+\\ \hfill \frac{C_1^{\prime }{(1+\delta )}^2}{\delta }\frac{\left(4r(k+p+2)log\left(1+\frac{{\parallel X\parallel }_F{\parallel Z\parallel }_F{\parallel \overline{M}\parallel }_F}{\sqrt{C_1}}\sqrt{\frac{mkp}{r(k+p)}}\right)+k+p+2log\frac{2}{{\epsilon }_m}\right)}{m}\left]\right\}.\end{array}$$

Then

$$\mathbb{E}\left[{\mathbb{P}}_{M\sim {\widehat{\rho }}_{\lambda }^{\prime }}(M\in {\Omega }_m)\right]\ge 1-{\epsilon }_m\underset{m\to \infty }{\to }1.$$

The proof of this theorem is provided in 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 Section 3, by a suitable version of the Langevin Monte Carlo (LMC) algorithm, a gradient-based sampling method. We propose to use a constant step-size unadjusted LMC algorithm; see [43] 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 }\left(M_k\right)+\sqrt{2h}{N}_{k}k=0,1,\dots$$

(8)

where $h>0$ is the step-size, ${\widehat{\rho }}_{\lambda }$ is the (quasi) posterior and $N_0,N_1,\dots$ are independent random matrices with independent and identical standard Gaussian entries. We provide a pseudo-code for LMC in Algorithm A1. For small values of the step-size h, the output of Algorithm A1, $\widehat{M}$, is very close to the integral (3) of interest. However, for some h that may not be small enough, the sum can explode [44]. In such cases, 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, we keep some control over 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 acceptance/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:

$$\begin{array}{c}\hfill {\tilde{M}}_{k+1}=M_k-h\nabla log{\widehat{\rho }}_{\lambda }\left(M_k\right)+\sqrt{2h}{N}_k,k=0,1,\dots ,\end{array}$$

(9)

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

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

(10)

where

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

is the transition probability density from x to $x^{\prime }$. The details of the Metropolis-adjusted Langevin algorithm (denoted by MALA) are presented in Algorithm A2. Compared to the 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 approximately $0.5$ following [45], while the step-size for LMC in the same setting should be smaller than the one for MALA [46].

4.2. Simulation Studies for Biliear Regression

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

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 independent and identically 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 independent and identically sampled from $\mathcal{N}(0,1)$.
• Model II: An approximate low-rank set up is studied. This series of simulations 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·B_1{B}_2^{\top }+U,$$

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

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

$${\widehat{M}}_{\mathrm{OLS}}={\left(X^{\top }X\right)}^{†}{X}^{\top }Y{Z}^{\top }{\left(Z{Z}^{\top }\right)}^{†}$$

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

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

$$\mathrm{Est}:=\parallel \widehat{M}-M^{*}{\parallel }_F^2/\left(pk\right),\mathrm{Nmse}:=\parallel \widehat{M}-M^{*}{\parallel }_F^2/{\parallel M^{*}\parallel }_F^2,$$

and the prediction error (Pred) as

$$\mathrm{Pred}:=\parallel X(\widehat{M}-M^{*}){Z\parallel }_F^2/\left(nq\right),$$

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

The results of our study are presented in

Table 1. 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
n = 100		Est	1.0053 (0.5480)	1.0342 (0.5559)	1.0052 (0.5478)
	$p=10$	Pred	0.1138 (0.0171)	0.0985 (0.0151)	0.1014 (0.0154)
		Nmse	0.4931 (0.1178)	0.5100 (0.1207)	0.4930 (0.1178)
		Est	1.3544 (0.5867)	1.3384 (0.5836)	1.3544 (0.5867)
	$p=100$	Pred	1.0066 (0.0430)	0.8761 (0.0756)	1.0030 (0.0424)
		Nmse	0.7049 (0.2944)	0.6963 (0.2927)	0.7049 (0.2944)
n = 1000		Est	1.0776 (0.5671)	1.0900 (0.5670)	1.0776 (0.5671)
	$p=10$	Pred	0.0099 (0.0013)	0.0099 (0.0013)	0.0099 (0.0013)
		Nmse	0.5185 (0.1198)	0.5264 (0.1219)	0.5185 (0.1198)
		Est	0.9662 (0.3240)	0.9688 (0.3244)	0.9662 (0.3240)
	$p=100$	Pred	0.0999 (0.0051)	0.0989 (0.0049)	0.0998 (0.0051)
		Nmse	0.4961 (0.1183)	0.4976 (0.1191)	0.4961 (0.1183)

and

Table 2. 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
n = 100		Est	4.0731 (1.828)	4.0989 (1.821)	4.0731 (1.828)
	$p=10$	Pred	0.1090 (0.0160)	0.0969 (0.0140)	0.0987 (0.0145)
		Nmse	0.5119 (0.1226)	0.5162 (0.1241)	0.5118 (0.1226)
		Est	4.6047 (1.812)	4.6038 (1.813)	4.6047 (1.812)
	$p=100$	Pred	1.0062 (0.0462)	1.0597 (0.0495)	1.0006 (0.0469)
		Nmse	0.5801 (0.1942)	0.5800 (0.1941)	0.5801 (0.1942)
n = 1000		Est	3.6733 (1.606)	3.6884 (1.606)	3.6733 (1.606)
	$p=10$	Pred	0.0098 (0.0015)	0.0098 (0.0015)	0.0098 (0.0015)
		Nmse	0.4812 (0.1271)	0.4835 (0.1260)	0.4813 (0.1271)
		Est	3.9972 (1.375)	3.9986 (1.376)	3.9972 (1.375)
	$p=100$	Pred	0.1000 (0.0043)	0.1032 (0.0057)	0.0999 (0.0043)
		Nmse	0.5013 (0.1061)	0.5014 (0.1063)	0.5013 (0.1062)

. From the tables, it can be observed that our proposed methods perform similarly to the OLS method. However, the estimation method obtained from the MALA algorithm often results in smaller prediction errors, particularly in high-dimensional settings. This advantage is even more pronounced when the method is applied in the context of inductive matrix completion, as discussed 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, Section 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_{\mathrm{miss}}$.

As in the context of inductive matrix completion, we only observe the response matrix with missing entries, $Y_{\mathrm{miss}}$, and thus we cannot construct the OLS estimator as in the case of bilinear regression. For this purpose, we first impute the missing entries in $Y_{\mathrm{miss}}$ by using the R package softImpute [48], where the rank of $M^{*}$ is specified as the true rank for matrix $Y_{\mathrm{miss}}$. We denote the resulting imputed matrix by $Y_{\mathrm{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 follows:

$${\widehat{M}}_{\mathrm{OLS}\_\mathrm{imp}}={\left(X^{\top }X\right)}^{†}{X}^{\top }{Y}_{\mathrm{imp}}{Z}^{\top }{\left(Z{Z}^{\top }\right)}^{†}$$

where $A^{†}$ denotes the Moore–Penrose inverse of matrix A. The LMC and MALA methods are initiated at the OLS_imp estimator and are run with 10,000 iterations, where the first 1000 steps are removed as burn-in periods.

As previously discussed in Section 4.2, we present the averages and the standard deviation of the mean squared estimation error (Est), the normalized (relative) mean square error (Nmse), and the prediction error (Pred) over 100 data replications in our results.

The results are detailed in

Table 3. 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
n = 100 $\kappa$ = 10%		Est	1.0559 (0.5060)	1.0803 (0.5122)	1.0559 (0.5060)
	$p=10$	Pred	0.1028 (0.0193)	0.1082 (0.0143)	0.1020 (0.0197)
		Nmse	0.4986 (0.1116)	0.5139 (0.1197)	0.4986 (0.1116)
		Est	1.4008 (0.8555)	1.3987 (0.8542)	1.4009 (0.8555)
	$p=100$	Pred	1.2250 (0.4568)	1.4468 (0.4137)	1.2252 (0.4570)
		Nmse	0.7148 (0.3591)	0.7136 (0.3581)	0.7148 (0.3591)
n = 100 $\kappa$ = 30%		Est	1.0432 (0.4963)	1.0917 (0.5085)	1.0432 (0.4963)
	$p=10$	Pred	0.2402 (0.2705)	0.1447 (0.0204)	0.2446 (0.2780)
		Nmse	0.5242 (0.1257)	0.5538 (0.1335)	0.5242 (0.1257)
		Est	1.6242 (0.8179)	1.6224 (0.8169)	1.6242 (0.8179)
	$p=100$	Pred	9.8879 (14.11)	10.807 (13.84)	9.8901 (14.11)
		Nmse	0.7993 (0.3340)	0.7985 (0.3334)	0.7993 (0.3340)
n = 1000 $\kappa$ = 10%		Est	0.9810 (0.4532)	0.9882 (0.4478)	0.9810 (0.4532)
	$p=10$	Pred	0.0114 (0.0033)	0.0112 (0.0015)	0.0114 (0.0033)
		Nmse	0.4933 (0.1076)	0.4984 (0.1075)	0.4933 (0.1076)
		Est	1.0063 (0.3465)	1.0088 (0.3471)	1.0063 (0.3465)
	$p=100$	Pred	0.1902 (0.1758)	0.1116 (0.0049)	0.1902 (0.1759)
		Nmse	0.5069 (0.1049)	0.5082 (0.1050)	0.5069 (0.1049)
n = 1000 $\kappa$ = 30%		Est	1.0110 (0.4886)	1.0223 (0.4872)	1.0110 (0.4886)
	$p=10$	Pred	0.0539 (0.0599)	0.0141 (0.0019)	0.0540 (0.0599)
		Nmse	0.5129 (0.1030)	0.5206 (0.1043)	0.5129 (0.1030)
		Est	1.0291 (0.3567)	1.0312 (0.3555)	1.0291 (0.3567)
	$p=100$	Pred	1.7529 (1.914)	0.1475 (0.0078)	1.7530 (1.913)
		Nmse	0.5054 (0.1055)	0.5067 (0.1053)	0.5054 (0.1055)

and

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
n = 100 imis 10%		Est	3.8319 (1.691)	3.8749 (1.719)	3.8319 (1.690)
	$p=10$	Pred	0.1604 (0.1271)	0.1092 (0.0153)	0.1598 (0.1322)
		Nmse	0.5116 (0.1154)	0.5169 (0.1147)	0.5116 (0.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 (0.3454)	0.7307 (0.3455)	0.7313 (0.3454)
n = 100 imis 30%		Est	4.1838 (1.850)	4.2535 (1.859)	4.1839 (1.850)
	$p=10$	Pred	0.7221 (0.7562)	0.1498 (0.0183)	0.7371 (0.7741)
		Nmse	0.5182 (0.1128)	0.5283 (0.1147)	0.5182 (0.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 (0.3821)	0.8994 (0.3820)	0.8998 (0.3821)
n = 1000 imis 10%		Est	3.9618 (1.678)	3.9788 (1.677)	3.9618 (1.678)
	$p=10$	Pred	0.0409 (0.0269)	0.0110 (0.0015)	0.0409 (0.0269)
		Nmse	0.4968 (0.1196)	0.4989 (0.1195)	0.4968 (0.1196)
		Est	4.1153 (1.295)	4.1163 (1.294)	4.1153 (1.295)
	$p=100$	Pred	1.0250 (0.9988)	0.1135 (0.0051)	1.0250 (0.9988)
		Nmse	0.5060 (0.1096)	0.5062 (0.1096)	0.5060 (0.1096)
n = 1000 imis 30%		Est	4.1647 (1.990)	4.1836 (1.995)	4.1647 (1.990)
	$p=10$	Pred	0.4615 (0.3497)	0.0141 (0.0017)	0.4616 (0.3498)
		Nmse	0.4905 (0.1157)	0.4933 (0.1171)	0.4905 (0.1157)
		Est	4.0578 (1.400)	4.0565 (1.397)	4.0578 (1.400)
	$p=100$	Pred	8.5608 (6.419)	0.1538 (0.0069)	8.5609 (6.419)
		Nmse	0.4944 (0.1184)	0.4943 (0.1180)	0.4944 (0.1184)

. It is evident from these tables that the results obtained from our MALA method surpass those of the other methods in terms of prediction error in most of the settings considered. This advantage becomes more pronounced as the missing rate in the response matrix increases. Additionally, it is worth noting that our MALA method is robust and performs well in the approximate low-rank setting (model II), while the OLS and LMC methods do not.

5. Discussion and Conclusions

In this paper, we focus on the problem of bilinear regression and its extension, the problem of inductive matrix completion, where the response matrix contains missing data. We propose a novel approach that combines elements of Bayesian statistics with a quasi-likelihood method. Our proposed method first addresses the problem of bilinear regression using a quasi-Bayesian approach and then adapts this approach to the problem of inductive matrix completion. By making use of a low-rankness assumption and leveraging the powerful PAC-Bayes bound technique, we provide statistical properties for our proposed estimators and for the quasi-posteriors.

Our proposed method includes an efficient gradient-based sampling algorithm that is designed to sample from the (quasi) posterior distribution. This algorithm allows for the approximate computation of mean estimators. These methods, referred to as LMC and MALA, were tested in various simulation studies and were found to perform well when compared to the ordinary least squared method. The ability to accurately sample from the (quasi) posterior distribution and compute mean estimators makes these methods a valuable tool for data analysis and modeling.

There are still some unresolved issues that require further investigation. One of these is the presence of missing data in the covariate matrices X and Z. This can have a significant impact on the analysis and may lead to biased results. Another area that needs further exploration is the assumption of independence and identically distributed data. In some cases, this assumption may not hold, and alternative models that allow for a dispersion matrix may be needed. These are potential topics for future research to address and further our understanding of these issues.