Low-rank model with covariates for count data analysis

TL;DR

This paper introduces LORI, a comprehensive low-rank modeling approach with covariates for count data, providing theoretical guarantees, an algorithm, and software, and demonstrating superior performance in estimation and imputation tasks.

Contribution

The paper presents LORI, a novel methodology combining a Poisson model, an algorithm, and automatic regularization selection for count data with covariates, including theoretical error bounds.

Findings

LORI outperforms existing methods in estimation accuracy.

LORI effectively imputes missing count data.

The method is validated on ecological and biological datasets.

Abstract

Count data are collected in many scientific and engineering tasks including image processing, single-cell RNA sequencing and ecological studies. Such data sets often contain missing values, for example because some ecological sites cannot be reached in a certain year. In addition, in many instances, side information is also available, for example covariates about ecological sites or species. Low-rank methods are popular to denoise and impute count data, and benefit from a substantial theoretical background. Extensions accounting for covariates have been proposed, but to the best of our knowledge their theoretical and empirical properties have not been thoroughly studied, and few softwares are available for practitioners. We propose a complete methodology called LORI (Low-Rank Interaction), including a Poisson model, an algorithm, and automatic selection of the regularization parameter, to analyze count tables with covariates. We also derive an upper bound on the estimation error. We provide a simulation study with synthetic data, revealing empirically that LORI improves on state of the art methods in terms of estimation and imputation of the missing values. We illustrate how the method can be interpreted through visual displays with the analysis of a well-know plant abundance data set, and show that the LORI outputs are consistent with known results. Finally we demonstrate the relevance of the methodology by analyzing a water-birds abundance table from the French national agency for wildlife and hunting management (ONCFS). The method is available in the R package lori on the Comprehensive Archive Network (CRAN).

Figures

Tables

Table 1. Table 1 : Estimation error (RMSE) of regression coefficients ‖ α ^ − α ∗ ‖ 2 2 + ‖ β ^ − β ∗ ‖ 2 2 superscript subscript norm ^ 𝛼 superscript 𝛼 2 2 superscript subscript norm ^ 𝛽 superscript 𝛽 2 2 \sqrt{\|\hat{\alpha}-\alpha^{*}\|_{2}^{2}+\|\hat{\beta}-\beta^{*}\|_{2}^{2}} of LORI and a Poisson GLM, for decreasing values of τ = ‖ Θ ‖ F / ‖ X 0 ‖ F 𝜏 subscript norm Θ 𝐹 subscript norm subscript 𝑋 0 𝐹 \tau=\|\Theta\|_{F}/\|X_{0}\|_{F} .

$\tau$		Mean of RMSE$\ast 100$	Standard deviation of RMSE$\ast 100$
$1$	LORI	$52$	$1.6$
	GLM	$143$	$29$
$0.5$	LORI	$17$	$0.9$
	GLM	$22$	$1.0$
$0.25$	LORI	$4.3$	$0.37$
	GLM	$5.3$	$0.33$
$0.1$	LORI	$1.7$	$0.34$
	GLM	$2.6$	$0.34$
$0$	LORI	$0.88$	$0.2$
	GLM	$0.87$	$0.19$

Table 2. Table 2 : Estimation error (Relative RMSE) of parameter matrix ‖ X ^ − X ∗ ‖ F / ‖ X ∗ ‖ F subscript norm ^ 𝑋 superscript 𝑋 𝐹 subscript norm superscript 𝑋 𝐹 \|\hat{X}-X^{*}\|_{F}/\|X^{*}\|_{F} of LORI and a Poisson GLM, for decreasing values of τ = ‖ Θ ‖ F / ‖ X 0 ‖ F 𝜏 subscript norm Θ 𝐹 subscript norm subscript 𝑋 0 𝐹 \tau=\|\Theta\|_{F}/\|X_{0}\|_{F} .

$\tau$		Mean of Relative RMSE$\ast 100$	Standard deviation of Relative RMSE$\ast 100$
$1$	LORI	$82$	$1.8$
	Poisson LRM	$95$	$4.4$
$0.5$	LORI	$40$	$0.25$
	Poisson LRM	$50$	$0.17$
$0.25$	LORI	$24$	$0.12$
	Poisson LRM	$40$	$0.12$
$0.1$	LORI	$11$	$0.081$
	Poisson LRM	$36$	$1.5$
$0$	LORI	$4.3$	$0.1$
	Poisson LRM	$34$	$1.2$

Table 3. Table 3 : Main effect of the Aravo environment characteristics estimated with LORI. The regularization parameter is tuned using QUT.

Aspect	Slope	PhysD	Snow
0.04	0.07	-0.02	-0.07

Table 4. Table 4 : Main effect of the Aravo species traits estimated with LORI. The regularization parameter is tuned using QUT.

Height	Spread	Angle	Area	Thick	SLA	Nmass	Seed
0.09	-0.24	-0.18	-0.20	-0.11	-0.17	0.18	-0.12

Table 5. Table 5 : Main effect of the sites characteristics estimated with LORI. The regularization parameter is tuned using QUT.

Agricultural surface	Latitude	Dist. to town	Dist. to coast	Surface
0.09	-0.21	0.09	-0.48	0.20

Table 6. Table 6 : Main effect of the years characteristics estimated with LORI.

Spring N/O	Spring N/E	Winter S/O	Winter S/E
-0.05	-0.03	-0.05	-0.03

Equations

$$Y_{ij}\sim \mathcal{P}({\mathrm{e}}^{X_{ij}^{\ast }}),(i,j)\in [n]\times [p],$$

$$X_{ij}^{\ast }={\mu }^{\ast }+{\alpha }_i^{\ast }+{\beta }_j^{\ast }+{\Theta }_{ij}^{\ast }\text{, }\mathrm{rank}({\Theta }^{\ast })\le \min (n-1,p-1),$$

$$X_{ij}^{\ast }={\mu }^{\ast }+\underset{\text{row effect}}{\underbrace{\sum _{k=1}^{K_1}{R}_{ik}{\alpha }_k^{\ast }}}+\underset{\text{column effect}}{\underbrace{\sum _{l=1}^{K_2}{C}_{il}{\beta }_l^{\ast }}}+{\Theta }_{ij}^{\ast },\mathrm{rank}({\Theta }^{\ast })\le \min (n-1,p-1).$$

$$X_{ij}^{\ast }={\mu }^{\ast }+{\alpha }_i^{\ast }+{\beta }_j^{\ast }+{ϵ}_{RC}{R}_i{C}_j,$$

$${\mathrm{e}}^{-\gamma }\le \mathbb{E}[Y_{ij}]\le {\mathrm{e}}^{\gamma }\text{ and }{\sigma }_{-}^2\le \mathrm{var}[Y_{ij}]\le {\sigma }_{+}^2.$$

$$X_{ij}^{\ast }={\mathrm{argmin}}_{x\in \mathbb{R}}\{-\mathbb{E}[Y_{ij}]x+\exp (x)\}.$$

$$r=\max \left(\left\{\mathrm{rank}(A):A\in {\mathcal{X}}_0\right\}\right).$$

$$X^{\ast }=X_0^{\ast }+{\Theta }^{\ast },X_0^{\ast }\in {\mathcal{X}}_0,{\Theta }^{\ast }\in \mathcal{T}.$$

$${\mathcal{X}}_0=\left\{{\left(\mu +\sum _{k=1}^{K_1}{R}_{ik}{\alpha }_k+\sum _{k=2}^{K_2}{C}_{ik}{\beta }_k\right)}_{(i,j)\in [n]\times [p]};\mu \in \mathbb{R},\alpha \in {\mathbb{R}}^{K_1},\beta \in {\mathbb{R}}^{K_2}\right\}.$$

$$\min \left\{{\pi }_{ij}:(i,j)\in [n]\times [p]\right\}=\pi >0.$$

$$\max \left(\left\{{\pi }_{i.}:i\in [n]\right\}\cup \left\{{\pi }_{.j}:j\in [p]\right\}\right)\le \beta .$$

$$\mathcal{L}(X)=\sum _{(i,j)\in [n]\times [p]}{\omega }_{ij}\left\{-Y_{ij}{X}_{ij}+\exp (X_{ij})\right\}.$$

$$({\hat{X}}_0,\hat{\Theta })$$

$$\Vert X_0+\Theta {\Vert }_{\infty }\le \gamma ,(i,j)\in [n]\times [p]$$

$$X_0\in {\mathcal{X}}_0,\Theta \in \mathcal{T}.$$

$${\Sigma }_R=\sum _{(i,j)\in [n]\times [p]}{ϵ}_{ij}{\omega }_{ij}{E}_{ij}.$$

$${\Vert X^{\ast }-\hat{X}\Vert }_F^2\le \frac{C}{{\pi }^2}\left(\left[\frac{{\lambda }^2}{{\sigma }_{-}^4}+(\mathbb{E}\Vert {\Sigma }_R\Vert {)}^2{\gamma }^2\right](\mathrm{rank}({\Theta }^{\ast })+r)+\log (n+p)\right),$$

$$\mathbb{E}\left[\exp (|Y_{ij}|/\delta )\right]<+\infty .$$

$${\Phi }_1$$

$${\Phi }_2$$

$${\Phi }_3$$

$$\lambda =\max \left\{4{\sigma }_{+}\sqrt{3\beta \log (n+p)},12\delta (e-1)\log \left(1+8{\delta }^2\frac{np}{\beta {\sigma }_{-}^2}\right)\log (n+p)\right\}.$$

$${\Vert X^{\ast }-\hat{X}\Vert }_F^2\le \frac{C}{{\pi }^2}\left\{\left(\max ({\Phi }_1,{\Phi }_2)+{\Phi }_3\right)\left(\mathrm{rank}({\Theta }^{\ast })+r\right)+\log (n+p)\right\},$$

$$X_{ij}^{\ast }={\mu }^{\ast }+\sum _{k=1}^{K_1}{R}_{ik}{\alpha }_k^{\ast }+\sum _{l=1}^{K_2}{C}_{il}{\beta }_l^{\ast }+{\Theta }_{ij}^{\ast },$$

$$\mathcal{F}(\mu ,\alpha ,\beta ,\Theta )=\mathcal{L}((\mu +R_{i,.}\alpha +C_{j,.}\beta {)}_{i,j}+\Theta ),$$

$$\tilde{\mathcal{F}}(\mu ,\alpha ,\beta ,\Theta )=\sum _{(i,j)\in [n]\times [p]}{\tilde{\omega }}_{ij}\{-Y_{ij}(\mu +R_{i,.}\alpha +C_{j,.}\beta +{\Theta }_{ij})+\exp (\mu +R_{i,.}\alpha +C_{j,.}\beta +{\Theta }_{ij})\}.$$

$$(\hat{\mu },\hat{\alpha },\hat{\beta },\hat{\Theta })$$

$$\Theta \in \mathcal{T}$$

$${\lambda }_0(Y)=\Vert \nabla \mathcal{L}({\hat{X}}_0)\Vert ,\mathrm{where}{\hat{X}}_0\in {\mathrm{argmin}}_{X\in {\mathcal{X}}_0}\mathcal{L}(X).$$

$${\mathbf{H}}_0:{\Theta }^{\ast }=0\text{ against the alternative }{\mathbf{H}}_1:{\Theta }^{\ast }\ne 0$$

$$\exp ({\hat{X}}_{ij})=\exp (\hat{\mu })\exp (R_{i,.}\hat{\alpha })\exp (C_{j,.}\hat{\beta })\exp ({\hat{\Theta }}_{ij}).$$

$$\frac{{\sigma }_{-}^2{\sum }_{(i,j)\in \Omega }({\hat{X}}_{ij}-X_{ij}^{\ast }{)}^2}{2}\le \underset{I}{\underbrace{-⟨\nabla \mathcal{L}(X^{\ast }),\hat{X}-X^{\ast }⟩}}+\underset{II}{\underbrace{\lambda (\Vert {\Theta }^{\ast }{\Vert }_{\ast }-\Vert \hat{\Theta }{\Vert }_{\ast })}}.$$

$$I\le \Vert \nabla \mathcal{L}(X^{\ast })\Vert \left(\Vert {\mathcal{P}}_{{\Theta }^{\ast }}(\hat{\Theta }-{\Theta }^{\ast }){\Vert }_{\ast }+\Vert {\mathcal{P}}_{{\Theta }^{\ast }}^{\perp }(\hat{\Theta }-{\Theta }^{\ast }){\Vert }_{\ast }+\Vert {\hat{X}}_0-X_0^{\ast }{\Vert }_{\ast }\right).$$

Full text

We propose a complete methodology called LORI (Low-Rank Interaction), including a Poisson model, an algorithm, and automatic selection of the regularization parameter, to analyze count tables with covariates. We also derive an upper bound on the estimation error. We provide a simulation study with synthetic data, revealing empirically that LORI improves on state of the art methods in terms of estimation and imputation of the missing values. We illustrate how the method can be interpreted through visual displays with the analysis of a well-know plant abundance data set, and show that the LORI outputs are consistent with known results. Finally we demonstrate the relevance of the methodology by analyzing a water-birds abundance table from the French national agency for wildlife and hunting management (ONCFS). The method is available in the R package lori on the Comprehensive Archive Network (CRAN).

Keywords: Count data; Dimensionality reduction; Ecological data; Imputation; Low-rank matrix recovery; Quantile universal threshold

1 Introduction

Let $Y$ be an $n\times p$ observation matrix of counts, $R\in\mathbb{R}^{n\times K_{1}}$ and $C\in\mathbb{R}^{p\times K_{2}}$ be two matrices containing row and column covariates, respectively. In our ecological application in Section 6, rows of the contingency table represent ecological sites, and columns represent years. For $(i,j)\in\{1,\ldots,n\}\times\{1,\ldots,p\}$, $Y_{ij}$ counts the abundance of water-birds measured in site $i$ during the year $j$. The row feature $R_{i\ell}$, $\ell\in\{1,\ldots,K_{1}\}$ embeds geographical information about the site $i$ (latitude, longitude, distance to coast, etc.) while the column feature $C_{j\ell}$, $\ell\in\{1,\ldots,K_{2}\}$ codes meteorological characteristics of the year $j$ (precipitation, etc.). In addition, some entries of $Y$ are missing. For example ecological sites are sometimes inaccessible because of meteorological or political conditions, and therefore cannot be counted.

Such count tables are often analyzed using low-rank models (Greenacre, 1984; Goodman, 1985; de Falguerolles, 1998; Christensen, 2010; Gower et al., 2011; Fithian and Josse, 2017), imposing a low-rank structure to an underlying parameter matrix. We assume a probabilistic framework with independent entries $Y_{ij}$ following a Poisson model

[TABLE]

and focus on the estimation of $X^{*}$ based on a low-rank assumption. The generalized additive main effects and multiplicative interaction model, or row-column model (see, e.g., Goodman (1985); de Falguerolles (1998)), assuming

[TABLE]

is adapted to do so. In this model, $\mu^{*}$ is an offset, the terms which only depend on the index of the row or column ($\alpha^{*}_{i}$ and $\beta^{*}_{j}$) are called main effects, and the terms which depend on both (here $\Theta^{*}_{ij}$) are called interactions (Kateri, 2014, Section 4.1.2, p.87).

A natural idea to incorporate covariates in this framework, is to express the row and column effects $\alpha^{*}_{i}$ and $\beta^{*}_{j}$ as regression terms on the covariates. In other words, for $\mu^{*}\in\mathbb{R}$, $\alpha^{*}\in\mathbb{R}^{K_{1}}$, $\beta^{*}\in\mathbb{R}^{K_{2}}$ and $\Theta^{*}\in\mathbb{R}^{n\times p}$,

[TABLE]

Such an extension is useful in practice for two main reasons. First, estimated covariates coefficients (and in particular their signs) can be used to determine whether the studied covariates have positive or negative effects on the counts; this is particularly useful in ecology, to check whether meteorological, geographical or political conditions favor or endanger species. Second, when the proportion of missing values is large, which is often the case in bird monitoring, incorporating (relevant) covariates can improve the imputation significantly.

Models related to (3) have been considered for statistical ecology applications in (Brown et al., 2014; ter Braak et al., 2017). However, to the best of our knowledge, their theoretical and empirical properties have not been thoroughly studied. On the other hand, the literature on convex low-rank matrix estimation is abundant and benefits from a substantial theoretical background, but few software with ready to use solution are available for practitioners, and applications for count data outside image analysis (Luisier et al., 2011; Salmon et al., 2014; Cao and Xie, 2016) and recommendation systems (Gopalan et al., 2014b) have not been attempted. The scope of this paper is to develop a complete methodology for the inference of model (3), bridging the gap between convex low-rank matrix completion and model-based count data analysis.

After detailing related work in Section 1.1, we introduce in Section 2 a general model which includes (3); we propose an estimation procedure through the minimization of a data fitting term penalized by the nuclear norm of the matrix $\Theta$, which acts as a convex relaxation of the rank constraint. Building up on existing results on nuclear norm regularized loss functions, we derive statistical guarantees in Section 2.1. In particular, we provide an upper bound for the Frobenius norm of the estimation error. In Section 3, we propose an optimization algorithm, and two methods to choose the regularization parameter automatically. We provide a simulation study in Section 4 revealing that LORI outperforms state-of-the-art methods when the proportion of missing values is large and the interactions are of significant order compared to the main effects. In Section 5, we show on plant abundance data with side information, how the results of our procedure can be interpreted through visual displays. In particular, the arising interpretation is consistent with known results from the original study (Choler, 2005). In Section 6, we use LORI to analyze a water-birds abundance data set from the French national agency for wildlife and hunting management (ONCFS). The proofs of the statistical guarantees are postponed to the appendix, and the method is available as an R package (R Core Team, 2016) called lori (LOw-Rank Interaction) on the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=lori.

1.1 Related work

Model (3) is closely related to other models previously suggested in the statistical ecology literature to analyze count tables with row and column covariates. For instance, Brown et al. (2014) and ter Braak et al. (2017) suggested the following model:

[TABLE]

with $R_{i}$, $1\leq i\leq n$ a row trait and $C_{j}$, $1\leq j\leq p$ a column trait. The interaction between covariates is modeled by $\epsilon_{RC}R_{i}C_{j}$, where $\epsilon_{RC}$ is an unknown parameter measuring the strength of the interaction between the two traits. The main difference with model (3) is that we incorporate the covariates in the main effects rather than the interactions, which leads to different interpretations. In terms of estimation properties, the main advantage of (3) is that, as long as $K_{1}\leq n$ and $K_{2}\leq p$, we estimate less parameters. This is an important point for us since in many applications we are interested in (see e.g. Section 6), a large proportion of entries is missing, limiting the amount of available data. Finally, model (4) was developed with the aim of testing significant associations between covariates, and its theoretical and empirical estimation properties, as far as we know, were not studied.

In the low-rank matrix completion literature, related approaches for count matrix recovery and dimensionality reduction can be embedded within the framework of low-rank exponential family estimation (Collins et al., 2001; de Leeuw, 2006; Li and Tao, 2013; Josse and Wager, 2016; Liu et al., 2016) as well as its Bayesian counterpart (Mohamed et al., 2009; Gopalan et al., 2014a). In terms of statistical guarantees, the theoretical performance of nuclear norm penalized estimators for Poisson denoising has been studied in Cao and Xie (2016), where the authors prove uniform bounds on the empirical error risk. Estimation rates are also given in Lafond (2015), where optimal bounds are proved for matrix completion in the exponential family. These two papers do not account for available covariates.

More recently, Chiquet et al. (2018) developped a probabilistic PCA framework for the exponential family, where covariates can be included in the parameter space. Fithian and Mazumder (2018) present a variety of low-rank problems including the generalized nuclear norm penalty (Angst et al., 2011), that can be used to include row and column covariates. Similar estimation problems were also considered, e.g., in Agarwal and Chen (2009); Abernethy et al. (2009). However, to the best of our knowledge, these papers did not provide statistical guarantees and the practical advantages of such extensions compared to classical low-rank methods have not been thoroughly studied.

2 General model and estimation

We now introduce a general version of the model described in the previous section. First, we relax the Poisson model and replace it with the following assumption on the distribution of $Y_{ij}$, $(i,j)\in[n]\times[p]$.

H 1**.**

The random variables $Y=\{Y_{i,j}\}_{(i,j)\in[n]\times[p]}$ are independent and there exist $\gamma>0$, $\sigma_{-}>0$ and $\sigma_{+}<\infty$ such that for all $i\in[n]$ and $j\in[p]$

[TABLE]

We define $X^{*}_{ij}$ for all $(i,j)\in[n]\times[p]$ by

[TABLE]

In other words, we do not assume that the random variable $Y_{ij}$ follows a Poisson distribution. The target parameter $X_{ij}^{*}$ minimizes the Kullback-Leibler divergence between the distribution of $Y_{ij}$ and a Poisson distribution. Let us also generalize the decomposition introduced in (3). Let $S_{1}$ and $S_{2}$ be fixed linear subspaces of $\mathbb{R}^{n}$ and $\mathbb{R}^{p}$ respectively. Let $P_{1}$ and $P_{2}$ be the orthogonal projection matrices on $S_{1}$ and $S_{2}$, $\mathcal{P}^{\perp}:X\in\mathbb{R}^{n\times p}\mapsto P_{1}XP_{2}^{\top}$, $\mathcal{P}:X\in\mathbb{R}^{n\times p}\mapsto X-\mathcal{P}^{\perp}(X)$, $\mathcal{X}_{0}\subset\{X\in\mathbb{R}^{n\times p};\mathcal{P}^{\perp}(X)=0\}$ and $\mathcal{T}=\{X\in\mathbb{R}^{n\times p};\mathcal{P}(X)=0\}$. We denote

[TABLE]

Consider the following decomposition:

[TABLE]

Denote, for $m\geq 1$, $\mathbbm{1}_{m}$ the vector of ones of length $m$. Model (3) is included in (7) by setting $S_{1}=\{u\in\mathbb{R}^{n};\mathbbm{1}_{n}^{\top}u=0\}$, $S_{2}=\{v\in\mathbb{R}^{p};\mathbbm{1}_{p}^{\top}v=0\}$, and

[TABLE]

The dimension of this subspace is at most $1+K_{1}+K_{2}$ and the rank of a matrix in $\mathcal{X}_{0}$ is less that $3$. We finally consider a setting with missing observations. Denote by $\Omega\subset[n]\times[p]$ the set of observed entries: $(i,j)\in\Omega$ if and only if $Y_{ij}$ is observed. Define also the random variables $(\omega_{ij})$ such that $\omega_{ij}=1$ if $Y_{ij}$ is observed and $\omega_{ij}=0$ otherwise. We assume that $(\omega_{ij})$ and $Y$ are independent, and a Missing Completely At Random (MCAR) scenario (Little and Rubin, 2002) where $(\omega_{ij})$ are independent Bernoulli random variables. For $(i,j)\in\{1,\ldots,n\}\times\{1,\ldots,p\}$, we denote $\pi_{ij}=\mathbb{P}(\omega_{ij}=1)$. We assume the probability of observing any entry is positive, i.e. there exists $\pi>0$ such that

[TABLE]

For $j\in[p]$, denote by $\pi_{.j}=\sum_{i=1}^{n}\pi_{ij}$ the probability of observing an element in the $j$-th column. Similarly, for $i\in[n]$, denote by $\pi_{i.}=\sum_{j=1}^{p}\pi_{ij}$ the probability of observing an element in the $i$-th row. We define the following upper bound:

[TABLE]

We can now define our data-fitting term:

[TABLE]

Denote $\|\cdot\|$ the operator norm (the largest singular value), $\|\cdot\|_{\infty}$ the infinity norm (the largest entry in absolute value) and $\|\cdot\|_{*}$ the nuclear norm (the sum of singular values). Our estimator of model (3), for a given regularization parameter $\lambda>0$, is the minimizer of the data-fitting term (10) penalized by the nuclear norm of $\Theta$:

[TABLE]

Denote by $\nabla\mathcal{L}(X)=\sum_{(i,j)\in[n]\times[p]}\omega_{ij}\left\{-Y_{ij}+\exp({X}_{ij})\right\}E_{ij}$, where $(E_{ij})$ are the matrices of the canonical basis of $\mathbb{R}^{n\times p}$, the gradient of $\mathcal{L}$ at $X$. Denote also $\partial^{2}\mathcal{L}/\partial x_{ij}^{2}$ the second derivative of $\mathcal{L}$ with respect to the $(i,j)$-th coordinate. Consider the following condition:

H 2**.**

The function $\mathcal{L}$ is strongly convex and smooth on $[-\gamma-\varepsilon,\gamma+\varepsilon]^{n\times p}$ for some $\varepsilon>0$. There exist $\sigma_{-}>0$ and $\sigma_{+}<\infty$ such that for all $X\in[-\gamma-\varepsilon,\gamma+\varepsilon]^{n\times p}$ and $(i,j)\in[n]\times[p]$, $\sigma_{-}^{2}\leq\partial^{2}\mathcal{L}(X)/\partial x_{ij}^{2}\leq\sigma_{+}^{2}$.

2.1 Statistical guarantees

We now derive an upper bound on the Frobenius estimation error of estimator (11). Let $(\epsilon_{ij})$, $(i,j)\in[n]\times[p]$ be i.i.d. Rademacher random variables independent of $Y$ and $\Omega$ and define

[TABLE]

Theorem 1**.**

Assume H 1-2, and $\lambda\geq 2\|\nabla\mathcal{L}(X^{*})\|$. Then for all $n,p\geq 1$, with probability at least $1-8(n+p)^{-1}$,

[TABLE]

where $r$, $\gamma$ are defined in (6) and (11), $C$ is a numerical constant whose value can be found in the proof and which is independent of $n$, $p$ and $X^{*}$.

Proof.

The proof is postponed to Section 8.1. ∎

We then control $\mathbb{E}\|\Sigma_{R}\|$, and compute a value of $\lambda$ such that the condition $\lambda\geq 2\|\nabla\mathcal{L}(X^{*})\|$ holds with high probability. We will need the following additional assumption on the distribution of the counts:

H 3**.**

There exists $\delta>0$ such that for all $(i,j)\in[n]\times[p]$,

[TABLE]

Define the following quantities, with $C^{*}$ a numerical constant defined in Lemma 5 and $r$, $\beta$ and $\gamma$ defined in (6), (9) and (11) respectively:

[TABLE]

Theorem 2**.**

Assume H 1-H 3 and set

[TABLE]

Then with probability at least $1-10(n+p)^{-1}$,

[TABLE]

where $C$ is a numerical constant independent of $n$, $p$ and $X^{*}$.

Proof.

The proof is postponed to Section 8.2. ∎

We recover an upper bound of order $\operatorname{rank}(\Theta^{*})\beta/\pi^{2}$, which is classical in low-rank matrix estimation and completion (Klopp, 2014; Lafond, 2015) and equal to $\operatorname{rank}(\Theta^{*})\max(n,p)/\pi$ when the sampling is almost uniform ($c_{1}\pi\leq\pi_{ij}\leq c_{2}\pi$). The additional term $r\beta/\pi^{2}$ accounts for explicit modeling of the covariates in the main effects. The constant term appearing in bound (15) grows linearly with the upper bound $\sigma_{+}^{2}$ and quadratically with the inverse of $\sigma_{-}^{2}$. This means that by relaxing Assumption 1 to allow $\operatorname{var}(Y_{ij})$ to grow as fast as $\log(n+p)$ or decrease as fast as $1/\log(n+p)$, we only lose a log-polynomial factor in bound (15).

3 Algorithm and selection of $\lambda$

3.1 Optimization algorithm

In this section, we propose an algorithm to solve the initial (11) for the initial model

[TABLE]

using alternating minimization (Csiszár and Tusnády, 1984), which consists in updating $\mu$, $\alpha$, $\beta$ and $\Theta$ alternatively, each time along a descent direction. Note that, in the algorithm and the entire numerical section, we relax the constraint $|\mu+R_{i,.}\alpha+C_{j,.}\beta+\Theta_{ij}|\leq\gamma$. Indeed, this constraint is mainly required to obtain statistical guarantees, and we observed that in practice, for $\gamma$ sufficiently large, this constraint is never reached. Denote

[TABLE]

and $\nabla_{\Theta}\mathcal{F}$ the gradient of $\mathcal{F}$ with respect to $\Theta$ defined by $(\nabla_{\Theta}\mathcal{F}(\mu,\alpha,\beta,\Theta))_{ij}=-Y_{ij}+\exp(\mu+R_{i,.}\alpha+C_{j,.}\beta+\Theta_{ij})$ if $\omega_{ij}=1$ and $(\nabla_{\Theta}\mathcal{F}(\mu,\alpha,\beta,\Theta))_{ij}=0$ otherwise. At every iteration we solve three sub-problems. The sub-problem in $\mu$ can be solved in closed form at each iteration; the updates in $\alpha$ and $\beta$ can be done simultaneously by estimating a Poisson generalized linear model (which can be done using standard algorithms implemented in available libraries); the update in $\Theta$ is along the proximal gradient direction, with a step size tuned using backtracking line search. Denote by $\mathcal{D}_{\lambda}$ the soft-thresholding operator of singular values at level $\lambda$ (Cai et al., 2010, Section 2). The procedure is sketched in Algorithm 1.

Note that if $K_{1}+K_{2}>|\Omega|$, with $|\Omega|$ denoting the cardinality of $\Omega$, the update in $\alpha$ and $\beta$ does not have a unique solution. However in our targeted applications, typically $K_{1}+K_{2}\ll|\Omega|$. In the package lori, we additionally implemented a warm-start strategy (Friedman et al., 2007), which consists in solving (11) for a large value of $\lambda$, then sequentially decreasing $\lambda$ and solving the new problem using the previous estimate as a starting point. Thus, our implementation solves (11) for the entire regularization path at once. Note that, even though our theoretical guarantees require a MCAR mechanism, the estimation method still holds when entries are missing at random (Little and Rubin (2002), Section 1.3). Its imputation properties are illustrated in an ecological application in Section 6.

3.2 Automatic selection of $\lambda$

A common way to select the regularization parameter is cross-validation, which consists in erasing a fraction of the observed cells in $Y$, estimating a complete parameter matrix $\hat{X}$ for a range of $\lambda$ values, and choosing the parameter $\lambda$ that minimizes the imputation error. This can be performed directly using LORI without modifying the code. Indeed, let $(\tilde{\omega}_{ij})$ denote the weights in $\{0,1\}$ indicating which entries are observed after removing some of them for cross-validation, and denote

[TABLE]

The optimization problem becomes

[TABLE]

which can be solved using the method described in Section 3.1 (see Algorithm 1). However, cross-validation is computationally costly. We suggest an alternative method to cross-validation, inspired by Donoho and Johnstone (1994) and the work of Giacobino et al. (2017) on quantile universal threshold. In Theorem 3 below, we define the so-called null-thresholding statistic of estimator (11), a function of the data $\lambda_{0}({Y})$ for which the estimated interaction matrix $\hat{\Theta}^{\lambda_{0}(Y)}$ is null, and the same estimate $\hat{\Theta}^{\lambda}=0$ is obtained for any $\lambda\geq\lambda_{0}({Y})$.

Theorem 3** (Null-thresholding statistic).**

The estimated interaction matrix $\hat{\Theta}^{\lambda}$ for a regularization parameter $\lambda$ is null if and only if $\lambda\geq\lambda_{0}(Y)$, where $\lambda_{0}(Y)$ is the null-thresholding statistic

[TABLE]

Proof.

The proof is postponed to Section 8.3. ∎

Here, $\left\|.\right\|$ is the operator norm (the largest singular value). We propose a heuristic selection of $\lambda$ based on this null-thresholding statistic $\lambda_{0}(Y)$. To explain further the procedure, we first need to define the following test:

[TABLE]

which tests whether the parameter matrix $X^{*}$ can be explained only in terms of linear combinations of the measured covariates. For a probability $\varepsilon\in(0,1)$, consider the upper $\varepsilon$-quantile $\lambda_{\varepsilon}$ of the null-thresholding statistics, namely that satisfies $\mathbb{P}_{\mathbf{H}_{0}}\left(\lambda_{0}(Y)>\lambda_{\varepsilon}\right)<\varepsilon$. The test which consists in comparing the statistics $\lambda_{0}(Y)$ to $\lambda_{\varepsilon}$ is of level $1-\varepsilon$ for (18). This can be seen as an alternative to the $\chi^{2}$ test for independence which handles covariates. In practice we do not have access to the distribution under the null $\mathbb{P}_{\mathbf{H}_{0}}\left(\lambda_{0}(Y)<\lambda\right)$, but perform parametric bootstrap (Efron, 1979) to compute a proxy $\check{\lambda}_{\varepsilon}$. In practice we recommend $\varepsilon=0.05$ and use $\lambda_{\textsf{QUT}}:=\check{\lambda}_{.05}$, and refer to it in what follows as quantile universal threshold (QUT). This selection of the regularization parameter is essentially the universal threshold of Donoho and Johnstone (1994) extended to our setting.

4 Simulation study

4.1 Estimation

We simulate $Y\in\mathbb{N}^{300\times 30}$ under model (1)–(3), with $R\in\mathbb{R}^{500\times 3}$ and $C\in\mathbb{R}^{300\times 4}$ drawn from multivariate Gaussian distributions with mean [math] and block-diagonal covariance matrices $\Sigma_{R}$ and $\Sigma_{C}$ respectively. We set $\mu^{*}=1$, $\alpha^{*}=(2,0,0)$, $\beta^{*}=(-2,0,0,0)$ and $\Theta$ of rank $5$. We compare the performance of LORI in terms of estimation of the regression coefficients $\alpha$ and $\beta$, and compare it to a standard Poisson Generalized Linear Model (GLM) estimated with the glm function in R. We repeat the experiment $100$ times for decreasing values of the ratio $\tau=\|\Theta\|_{F}/\|X_{0}\|_{F}$, where $X_{0}=(\mu+R_{i}\alpha+C_{j}\beta)_{ij}$ is fixed. We look at the Root Mean Square Error (RMSE) for the estimation of $X_{0}$; the results are given in Table 1, where we observe that LORI and the Poisson GLM are equivalent for $\tau=0$, and that LORI outperforms the GLM for non-zero interactions, with a gap widening as $\tau$ increases.

Second, we compare LORI to a convex low-rank matrix estimatiob procedure with a Poisson loss function and where covariates are not modeled (e.g. Lafond (2015)), in terms of the relative estimation error $\|\hat{X}-X^{*}\|_{F}/\|X^{*}\|_{F}$ (because $\|X\|_{F}$ varies with $\tau$). We refer to this competitor as ”Poisson LRM”. Again, we reproduce the experiment $100$ times for decreasing values of the ratio $\tau=\|\Theta\|_{F}/\|X_{0}\|_{F}$.

On Table 2, we observe that LORI achieves lower errors than Poisson LRM, which is expected as we simulated under the LORI model. As $\tau$ decreases – i.e. the size of the main effects increases relative to the interactions – both errors decrease as well, and the gap between LORI and the Poisson LRM widens, indicating that modeling covariates explicitly improves the estimation.

4.2 Imputation

Using the same simulation scheme, we now compare LORI in terms of missing values imputation to Correspondence Analysis (CA) and Trends & Indices for Monitoring data Pannekoek and van Strien (2001) (TRIM), a method based on a Poisson log-linear model used to impute bird abundance data. To do so we erase an increasing proportion of entries in the data and impute them using LORI, CA and TRIM, replicating the experiment 100 times. We also impute the missing values using the column means, as a baseline referred to as ”MOY”. We observe on Figure 1 that LORI performs best, which is expected as we simulate under the LORI model. Moreover, the gap widens as the percentage of missing values increases. In particular, the error of TRIM for 80% of missing entries is not represented because the method fails (we use the default parameters of the R package rtrim).

5 Analysis of the Aravo data

The Aravo data set (Choler, 2005) counts the abundance of $82$ species of alpine plants in $75$ sites in France; covariates about the environments and species are also available. We focus on $8$ species traits providing physical information about plants (height, spread, etc.), and $4$ environmental variables giving geographical and meteorological information about sites. We apply our method LORI after scaling the covariates and tuning the regularization parameter with the QUT method. This results in estimates for the main effects of the environment characteristics $\alpha$ and of the species traits $\beta$, and of the interaction matrix $\Theta$.

The main effects of environment characteristics are given in Table 3 and the main effects of the species traits in Table 4. First we observe that overall, species traits have larger effects than environment characteristics on the observed abundances. In particular, the mass-based leaf nitrogen content (Nmass) has a large positive effect, which seems to indicate that plants with a large Nmass tend to be more abundant across all environments. On the other hand, the maximum lateral spread of clonal plants (Spread), area of single leaf (Area) leaf elevation angle estimated at the middle of the lamina (Angle) and specific leaf area (SLA) have large negative effects on the abundances.

The estimated rank of the interaction matrix $\hat{\Theta}$ (number of singular values above $10^{-6}$) is $2$. The environments (rows) and species (columns) can be visualized on a biplot (de Rooij and Heiser, 2005, Section 2.5), where rows and columns are represented simultaneously in a normalized Euclidean space. In such plots, the dimensions of the Euclidean space are given by the principal directions of $\hat{\Theta}$, scaled by the square root of the singular values of $\hat{\Theta}$. Figure 2 shows such a display, which can be interpreted in terms of distance between points: a species and an environment that are close interact highly, and two species or two environments that are close have similar profiles. Justifications for such a distance interpretation can be found in (de Rooij and Heiser, 2005, Section 2.5) or (Fithian and Josse, 2017, Section 2).

We can then look at the relations between the known traits and the interaction directions of $\hat{\Theta}$. Figure 3(a) shows that the two first directions of interaction are correlated with the species covariates; the correlation is particularly high for the Nmass and SLA variables. Thus, on Figure 2, the two directions separate the plants with large SLA and Nmass (top right corner) from those with small SLA and Nmass (bottom left corner). Then, Figure 3(a) shows that the directions of interaction are also correlated with the environment covariates, and particularly with the mean snowmelt date (Snow). Thus, on Figure 2, the two directions separate the late melting environments (top right corner) from the early melting environments (bottom left corner).

Combining the interpretation of Figure 2, Figure 3(a) and Figure 3(b), we deduce that plants with large Nmass and SLA interact highly with late melting sites (large value of Snow). This was in fact the main result obtained in the original study Choler (2005) (see, e.g., the summary of findings in the abstract), which advocates the good properties of LORI in terms of interpretation.

6 Using covariates to impute ecological data

The water-birds data count the abundance of migratory water-birds in 785 wetland sites (across the 5 countries in North Africa), between 1990 and 2017 (Sayoud et al., 2017). One of the objectives is to assess the effect of time on species abundances, to monitor the populations and assess wetlands conservation policies. Ornithologists have also recorded side information concerning the sites and years, which may influence the counts. For instance, meteorological anomalies, latitude and longitude. The count table contains a large amount of missing entries (70%), but the covariate matrices which contain respectively 6 covariates about the 785 sites and 8 covariates about the 18 years, are fully observed. Our method allows to take advantage of the available covariates to provide interpretation for spatio-temporal patterns. As a by-product, it produces an imputed contingency table.

Tables 5 and 6 show the estimated main effects of some of the sites and years characteristics. Sites with large latitudes are associated to smaller counts, as well as sites which are far from the coast. Sites which are located far from towns, and sites with large water surfaces are associated to larger counts. The four year covariates given in Table 6 concern meteorological anomalies. The associated coefficients are all negative, indicating that more important anomalies are associated to smaller abundances.

The sites and years can also be displayed using the same visual tools as described in Section 5. Figure 4(a) and 4(b) show the correlations between the covariates and the directions of interaction. On the two-dimensional display on Figure 4(c), the first dimension is correlated with geographical characteristics, and the second dimension with meteorological anomalies. We observe a very clear temporal gradient along the second dimension, indicating that over time, meteorological anomalies increase (in the sense of a summary anomaly variable embodied by the second direction). One of the site (378) lays out of the point cloud, and corresponds to a site with very large surface.

LORI also returns counts estimates, which can be used to compute an estimation of the total yearly abundances (i.e. counts estimates summed across sites). To better assess the temporal trend, one can decompose the estimated counts into three factors corresponding to the site effects, year effects and interactions respectively. Indeed, for $(i,j)\in\{1,\ldots,n\}\times\{1,\ldots,p\}$, one can write

[TABLE]

Figure 5 shows the last three factors of this decomposition separately.

On Figure 5(a) we see that most sites have multiplicative effects around $1$ on count scale. One site (site 378, large red point) stands out; again, it corresponds to an extremely large site ($6000$km2, $5$ times larger than the second, $300$ times larger than the mean). In this respect, the row effects act as normalization factors accounting for surface. We also observe tenuous levels along the $x$ axis, corresponding to sites of different countries. On Figure 5(b) we observe a decreasing temporal trend. This means that, all other things being equal, later years tend to produce smaller abundances. As illustrated in Figure 4(b), this temporal trend can be associated with the effects of meteorological anomalies. Note that the temporal effects (top right) are smaller in amplitude than the spatial effects (top left). Indeed, more variability is observed between sites in a given year, than between years for a given site.

Finally, looking at the interaction matrix on Figure 5(c), we see that the interaction is mainly driven by a few sites which interact more or less highly with every year. In particular the site $582$ presents large interactions with every year, and corresponds to a national park in Morocco and is the most abundant site ($4$ million birds in total, twice as much as the second most abundant, $120$ more than the mean). Here, the large abundance is not explained by the geographical covariates available (but maybe by other unmeasured factors such as protection legislations), and thus the extreme abundance is captured in the interaction rather than the main effects. This profile was also visible on Figure 4(c) where the site $582$ lies amongst the cloud of years, indicating large interactions through small Euclidean distances.

The site 704 also presents large interactions, but they are not constant throughout the years. It corresponds to Ichkeul national park in Tunisia, which is a major site for most species; the abundances are very large in Ichkeul compared to other sites. However during several years including 2007, bad weather conditions prevented ornithologist to correctly count the birds, thus reported counts are significantly lower than expected. This explains the drop in the interaction in 2007 for Ichkeul, corresponding to an outlier behavior. Again, such a profile could not be highlighted without modeling interactions. This illustrates one of the advantages of LORI for such bird abundance data compared to state-of-the-art methods such as Pannekoek and van Strien (2001) which do not model interactions. In particular, in most cases the interaction terms absorb outlying values (small or large), and indirectly account for the over-dispersion which is known to occur in birds abundance data.

7 Discussion

We conclude by discussing some opportunities for further research. To select covariates, we could penalize the main effects with an $\ell_{1}$ penalty on $\alpha$ and $\beta$. It may be also of interest to consider other sparsity inducing penalties. In particular, penalizing the Poisson log-likelihood by the absolute values of the coefficients of the interaction matrix $\Theta$ could possibly lead to solutions where some interactions are driven to [math] and a small number of large interactions are selected. Secondly, it would be useful to develop a multiple imputation procedure based on LORI, to provide confidence regions for the estimated parameters. The properties of the thresholding test, which can be seen as an alternative to a chi-squared test for independence with covariates, also merit further investigation. In particular, the power could be assessed. Finally, we could also explore whether our model could be extended to more complex models such as the zero-inflated negative binomial models.

Acknowledgments

The authors thank Trevor Hastie, Edgar Dobriban, Olga Klopp, Kevin Bleakey and Stéphane Dray for their very helpful comments on this manuscript. We would like to thank all the organizations involved in the Mediterranean Waterbirds Network (MNW) who provided the waterbird data set used in this article: the GREPOM/BirdLife Morocco, the Direction Générale des Forêts (Algeria), the AAO/BirdLife Tunisia, the Libyan Society for Birds, the Egyptian Environment Affairs Agency, the national agency for wildlife and hunting management (ONCFS, France) and the Research Institute of the Tour du Valat (France). We are also grateful to all the field observers who participated in the North African IWC thus making this data set so rich, and to Pierre Defos du Rau and Laura Dami for their help.

8 Proofs

8.1 Proof of Theorem 1

We will first derive an upper bound for $\sum_{(i,j)\in\Omega}(\hat{X}_{ij}-X_{ij}^{*})^{2}$, then control $\|\hat{X}-X^{*}\|_{F}^{2}$ by $\|\hat{X}-X^{*}\|_{F}^{2}\leq\sum_{(i,j)\in\Omega}(\hat{X}_{ij}-X_{ij}^{*})^{2}+\mathsf{D}$, with $\mathsf{D}$ a residual term defined later on. By definition of $\hat{X}=\hat{X_{0}}+\hat{\Theta}$, $\mathcal{L}(\hat{X})+\lambda\|\hat{\Theta}\|_{*}\leq\mathcal{L}(X^{*})+\lambda\|\Theta^{*}\|_{*}.$ Using the strong convexity of $\mathcal{L}$ and substracting $\langle\nabla\mathcal{L}(X^{*}),\hat{X}-X^{*}\rangle$ on both sides of this inequality, we obtain

[TABLE]

We will bound separately the two terms on the right hand side of (19).

Given a matrix $X\in\mathbb{R}^{n\times p}$, we denote $\mathcal{S}_{1}(X)$ (resp. $\mathcal{S}_{2}(X)$) the span of left (resp. right) singular vectors of $X$. Let $P_{\mathcal{S}_{1}(X)}^{\perp}$ (resp. $P_{\mathcal{S}_{2}(X)}^{\perp}$) be the orthogonal projector in $\mathbb{R}^{n}$ on $\mathcal{S}_{1}(X)^{\perp}$ (resp. in $\mathbb{R}^{p}$ on $\mathcal{S}_{2}(X)^{\perp}$). We define the projection operator in $\mathbb{R}^{n\times p}$ $\mathcal{P}_{X}^{\perp}:\tilde{X}\mapsto P_{\mathcal{S}_{1}(X)}^{\perp}\tilde{X}P_{\mathcal{S}_{2}(X)}^{\perp}$, and $\mathcal{P}_{X}:\tilde{X}\mapsto\tilde{X}-P_{\mathcal{S}_{1}(X)}^{\perp}\tilde{X}P_{\mathcal{S}_{2}(X)}^{\perp}$. We use the following Lemma, proved in (Lafond, 2015, Lemma 16).

Lemma 1**.**

For all $M$ and $M^{\prime}$ in $\mathbb{R}^{n\times p}$,

1. (i)

$\|M+\mathcal{P}_{M}^{\perp}(M)\|_{*}=\|M\|_{*}+\|\mathcal{P}_{M}^{\perp}(M)\|_{*},$** 2. (ii)

$\|M\|_{*}-\|M^{\prime}\|_{*}\leq\|\mathcal{P}_{M}(M-M^{\prime})\|_{*}-\|\mathcal{P}_{M}^{\perp}(M-M^{\prime})\|_{*},$** 3. (iii)

$\|\mathcal{P}_{M}(M-M^{\prime})\|_{*}\leq\sqrt{2\mathrm{rk}(M)}\|M-M^{\prime}\|_{F}.$**

Using $|\langle\nabla\mathcal{L}(X^{*}),\hat{X}-X^{*}\rangle|\leq\|\hat{X}-X^{*}\|_{*}\|\nabla\mathcal{L}(X^{*})\|$ and the triangular inequality gives that

[TABLE]

Then, Lemma 1 ii applied to $\hat{\Theta}$ and $\Theta^{*}$, results in

[TABLE]

Plugging inequalities (20) and (21) in (19) we obtain

[TABLE]

We now use the condition $\lambda\geq 2\|\nabla\mathcal{L}(X^{*})\|$ in (22):

[TABLE]

Then, $\text{rk}(\hat{X}_{0}-X^{*}_{0})\leq r$ and $\|\hat{X}_{0}-X^{*}_{0}\|_{F}\leq\|\hat{X}-X^{*}\|_{F}$ imply that $\|\hat{X}_{0}-X^{*}_{0}\|_{*}\leq\sqrt{r}\|X^{*}-\hat{X}\|_{F}$, which together with Lemma 1 iii and $\|\hat{\Theta}-\Theta^{*}_{0}\|_{F}\leq\|\hat{X}-X^{*}\|_{F}$ yields

[TABLE]

We now derive the upper bound $\|\hat{X}-X^{*}\|_{F}^{2}\leq\sum_{(i,j)\in[n]\times[p]}\omega_{ij}(\hat{X}_{ij}-X_{ij}^{*})^{2}+\mathsf{D}$. Define $\eta=72\log(n+p)/(\pi\log(6/5))$,

[TABLE]

and the set

[TABLE]

We start by showing in the following Lemma that whenever $\hat{X}-X^{*}$ belongs to $\mathcal{C}(\eta,\rho)$ (for $\rho$ and $\mathsf{D}$ defined later on), a restricted strong convexity property of the form $\|\hat{X}-X^{*}\|_{F}^{2}\leq\sum_{(i,j)\in[n]\times[p]}\omega_{ij}(\hat{X}_{ij}-X_{ij}^{*})^{2}+\mathsf{D}$ holds. Define

[TABLE]

Lemma 2**.**

Let $\eta=72\log(n+p)/(\pi\log(6/5))$ and $\rho>0$. With probability at least $1-8(n+p)^{-1}$, for all $X\in\mathcal{C}(\eta,\rho)$ we get

[TABLE]

with $\Sigma_{R}$ defined in (12).

Proof.

Consider the event

[TABLE]

Define also for $l\in\mathbb{N}_{*}$

[TABLE]

for $\kappa=6/5$ and $\eta=72\log(n+p)/(\pi\log(6/5))$. On $\mathcal{B}$, there exist $l\geq 1$ and $X\in\mathcal{C}(\eta,\rho)$ such that $X\in\mathcal{C}(\eta,\rho)\bigcap\mathcal{S}_{l}$, and

[TABLE]

For $T>0$, define the set

[TABLE]

and the event

[TABLE]

It follows from (28) that $\mathcal{B}\subset\bigcup_{l=1}^{+\infty}\mathcal{B}_{l}$; thus, it is enough to estimate the probability of the events $\mathcal{B}_{l}$, $l\in\mathbb{N}$, and then apply the union bound. Such an estimation is given in the following Lemma, adapted from Klopp (2015) (see Lemma 10). Define

[TABLE]

Lemma 3**.**

Under the assumptions of Theorem 1,

[TABLE]

where $\varsigma$ is defined in (27).

Proof.

We use the following Talagrand’s concentration inequality and a symmetrization argument. Recall the statement of Talagrand’s concentration inequality. Let $f:[-1,1]^{m}\mapsto\mathbb{R}$ a convex Lipschitz function with Lipschitz constant L, $\Xi_{1},\ldots,\Xi_{m}$ be independent random variables taking values in $[-1,1]$, and $Z:=f(\Xi_{1},\ldots,\Xi_{m})$. Then, for any $t\geq 0$, $\mathbb{P}(|Z-\mathbb{E}[Z]|\geq 16L+t)\leq 4\mathrm{e}^{-t^{2}/2L^{2}}$. For $\mathbf{x}=(x_{ij})$, $(i,j)\in[n]\times[p]$, we apply this result to the function

[TABLE]

which is Lipschitz with Lipschitz constant $\sqrt{\pi^{-1}T}$:

[TABLE]

where we have used $||a|-|b||\leq|a-b|$,$\|X\|_{\infty}\leq 1$ and $\mathbb{E}\left[\Sigma(\omega,X)\right]\leq T$. Thus, Talagrand’s inequality and the identity $\sqrt{\pi^{-1}T}\leq T/(2\times 96)+96/(2\pi)$ give

[TABLE]

Taking $t=T/6$ we get

[TABLE]

Now we bound the expectation $\mathbb{E}[Z_{T}]$ using a symmetrization argument (Ledoux, 2001, Section 7.2). Let $(\epsilon_{ij})$ be an i.i.d. Rademacher sequence. We have

[TABLE]

Then, the contraction inequality (see Koltchinskii (2011), Theorem 2.2) yields

[TABLE]

where $\Sigma_{R}$ is defined in (12). For $X\in\mathcal{C}(\eta,\rho,T)$ we have that $\|X\|_{*}\leq\sqrt{\rho\pi^{-1}T}$. Then by duality between the nuclear and operator norms we obtain

[TABLE]

Combined with (31) and using $8\sqrt{\rho\pi^{-1}T}\mathbb{E}\|\Sigma_{R}\|\leq\frac{T}{2\times 3}+\frac{3\times 8^{2}\rho\pi^{-1}}{2}(\mathbb{E}\|\Sigma_{R}\|)^{2}$ we finally obtain (30) using the definition of $\varsigma$ in (27). ∎

Lemma 3 implies that

[TABLE]

which concludes the proof. ∎

Case 1

If $\sum_{(i,j)\in[n]\times[p]}\pi_{ij}(\hat{X}_{ij}-X^{*}_{ij})^{2}\leq\eta$, then $\|\hat{X}-X^{*}\|_{2}^{2}\leq\eta/\pi$ and the result of Theorem 1 (15) is proved.

Case 2

If $\sum_{(i,j)\in[n]\times[p]}\pi_{ij}(\hat{X}_{ij}-X^{*}_{ij})^{2}>\eta$. Let us show that $(\hat{X}-X^{*})/2\gamma\in\mathcal{C}(\eta,64\operatorname{rank}(X^{*}))$. Using (22), $\sigma_{-}^{2}\sum_{(i,j)\in\Omega}(\hat{X}_{ij}-X^{*}_{ij})^{2}\geq 0$ and $\|\nabla\mathcal{L}(X^{*})\|\leq\lambda/2$, we obtain that

[TABLE]

On the other hand,

[TABLE]

Thus, Lemma 2 implies that with probability at least $1-8(n+p)^{-1}$,

[TABLE]

Combining (33) and (24) we obtain

[TABLE]

Finally, using the identity $ab\leq a^{2}+b^{2}/4$ and $\operatorname{rank}(X^{*})\leq\operatorname{rank}(\Theta^{*})+r$ we obtain

[TABLE]

8.2 Proof of Theorem 2

Theorem 2 derives from Theorem 1 and combining the two following steps: 1) computing a value of $\lambda$ such that the condition $\lambda\geq 2\|\nabla\mathcal{L}(X^{*})\|$ holds with high probability and 2) controlling $\mathbb{E}\|\Sigma_{R}\|$. Let us start with 1). Define the random matrices $Z_{ij}=\omega_{ij}(-Y_{ij}+\exp(X^{*}_{ij}))E_{ij}$ and the quantity

[TABLE]

Lemma 4**.**

Under the assumptions of Theorem 2,

[TABLE]

Proof.

For all $(i,j)\in[n]\times[p]$, $Z_{ij}Z_{ij}^{T}=\omega_{ij}(-Y_{ij}+\exp(X_{ij}^{*}))^{2}E_{ij}E_{ij}^{\top}$, and $\mathbb{E}[Z_{ij}Z_{ij}^{\top}]=\mathbb{E}[\omega_{ij}]\mathbb{E}[(-Y_{ij}+\exp(X_{ij}^{*}))^{2}]E_{ij}E_{ij}^{\top}$, which is a diagonal matrix with [math] everywhere except on the $i$-th element of its diagonal, where its value is $\mathbb{E}[\omega_{ij}]\mathbb{E}[(-Y_{ij}+\exp(X_{ij}^{*}))^{2}]$. Thus,

[TABLE]

is also a diagonal matrix, and the $i$-th element of its diagonal is $\sum_{j=1}^{p}\mathbb{E}[\omega_{ij}]\mathbb{E}[(-Y_{ij}+\exp(X_{ij}^{*}))^{2}]$. We obtain that

[TABLE]

Using $\mathbb{E}[Y_{ij}]=\exp(X_{ij}^{*})$ and $\sigma_{-}^{2}\leq\text{var}(Y_{ij})\leq\sigma_{+}^{2}$, we obtain:

[TABLE]

Using the same arguments, we also obtain

[TABLE]

Combining (37) and (37), we obtain that

[TABLE]

which concludes the proof. ∎

Note that $\mathbb{E}\left[Z_{ij}\right]=0$ for all $(i,j)\in[n]\times[p]$ and $\nabla\mathcal{L}(X^{*})=\sum_{i=1}^{n}\sum_{j=1}^{p}Z_{ij}.$ We use an extension of Theorem 4 in Koltchinskii (2013) to rectangular matrices via self-adjoint dilation (cf., for example, 2.6 in Tropp (2012)). Let $\Xi_{1},\ldots,\Xi_{m}$ be $m$ independent $(n\times p)$-matrices satisfying $\mathbb{E}[\Xi_{i}]=0$ and

[TABLE]

for some constant $M$ and for all $i\in\{1,\ldots,m\}$. Define

[TABLE]

and $\bar{U}=M\log(1+2\frac{M^{2}}{\sigma^{2}})$. Then, for $t\bar{U}\leq 2(e-1)\sigma^{2}m$,

[TABLE]

and for $t\bar{U}>2(e-1)\sigma^{2}m$,

[TABLE]

Under Assumption 3 we may apply this result with $m=np$, $(\Xi_{1},\ldots,\Xi_{m})=(Z_{11},\ldots,Z_{np})$, $M=2\delta$, $\sigma^{2}=\sigma_{Z}^{2}$ and $\bar{U}=2\delta\log(1+8\delta^{2}/\sigma_{Z}^{2})$. Taking

[TABLE]

and using Lemma 4, we get that with probability at least $1-(n+p)^{-1}$,

[TABLE]

Thus, taking $\lambda$ as in Theorem 2 ensures that $\lambda\geq 2\|\nabla\mathcal{L}(X^{*})\|$ with probability at least $1-(n+p)^{-1}$.

We now control $\mathbb{E}\|\Sigma_{R}\|$ with the following lemma.

Lemma 5**.**

There exists an absolute constant $C^{*}$ such that the two following inequality holds

[TABLE]

Proof.

We use an extension to rectangular matrices via self-adjoint dilation of Corollary 3.3 in Bandeira and van Handel (2016).

Proposition 1**.**

Let $A$ be an $n\times p$ rectangular matrix with $A_{ij}$ independent centered bounded random variables. then, there exists a universal constant $C^{*}$ such that

[TABLE]

[TABLE]

Applying Proposition 1 to $\Sigma_{R}$ with $\sigma_{1}\vee\sigma_{2}\leq\sqrt{\beta}/|\Omega|$ and $\sigma_{*}\leq 1$ we obtain

[TABLE]

∎

Combining 1) and 2) with (34) and a union bound argument, we obtain the result of Theorem 2.

8.3 Proof of Theorem 3

In what follows we denote for $X_{0}\in\mathcal{X}_{0}$ and $\Theta\in\mathcal{T}$ $\mathcal{F}^{\lambda}(X_{0},\Theta)=\mathcal{L}(X_{0}+\Theta)+\lambda\|\Theta\|_{*}.$ We establish below that $\lambda_{0}(Y)$ defined in (17) is equal to

[TABLE]

where for $\mathcal{K}\subset\mathbb{R}^{n\times p}$ ,$\chi_{\mathcal{K}}(X)$ is the characteristic function of the set $\mathcal{K}$, equal to [math] on $\mathcal{K}$ and $+\infty$ elsewhere, and $\hat{X}_{0}=\underset{X\in\mathcal{X}_{0}}{\operatorname{argmin}}\mathcal{L}(X)$ (see (17)). The subdifferential of the objective function $\mathcal{F}^{\lambda}$ with respect to $\Theta$ is given by

[TABLE]

$0\in\partial_{\Theta}\chi_{\mathcal{T}}(\Theta)\mid_{\Theta=0}$. Lemma 6 ensures that $0\in\partial\mathcal{F}^{\lambda}(\Theta)\mid_{\Theta=0}$ if and only if

[TABLE]

This is equivalent to $\lambda\geq\left\|\mathcal{P}_{\mathcal{T}}(\nabla\mathcal{L}(\hat{X}_{0}))\right\|.$ Additionally, at the optimum $\hat{X}_{0}$, we have $\mathcal{P}_{\mathcal{T}}(\nabla\mathcal{L}(\hat{X}_{0}))=\nabla\mathcal{L}(\hat{X}_{0})$, which concludes the proof.

Lemma 6**.**

Let $g:\mathcal{T}\rightarrow\mathbb{R}_{+}$ be the function defined by $g(A)=\left\|A\right\|_{*}$ for $A\in\mathcal{T}$. $\partial g(0)=\left\{W\in\mathbb{R}^{n\times p},\left\|\mathcal{P}_{\mathcal{T}}(W)\right\|\leq 1\right\}$.

Proof.

By definition of the subdifferential we need to prove that for all $W\in\mathbb{R}^{n\times p}$, $\left\|\mathcal{P}_{\mathcal{T}}(W)\right\|<1$, and for all $B\in\mathcal{T}$, $g(B)\geq g(0)+\left\langle W,B-0\right\rangle$. First $B\in\mathcal{T}$ implies $\langle W,B\rangle=\langle\mathcal{P}_{\mathcal{T}}(W),B\rangle$, therefore $\left\|\mathcal{P}_{\mathcal{T}}(W)\right\|\leq 1$ is a sufficient condition for $W\in\partial g(0)$. Now assume $\left\|\mathcal{P}_{\mathcal{T}}(W)\right\|>1$ and let $\mathcal{P}_{\mathcal{T}}(W)=U\Sigma V^{T}$, where $U$ and $V$ are orthogonal matrices of left and right singular vectors, and $\Sigma_{11}=\left\|\mathcal{P}_{\mathcal{T}}(W)\right\|>1$. Let us define $B=U\tilde{\Sigma}V^{T}$, $\tilde{\Sigma}_{11}=1$ and $\tilde{\Sigma}_{ij}=0$ elsewhere; note that with this definition $B\in\mathcal{T}$. We have $g(B)=1$ and $\langle\mathcal{P}_{\mathcal{T}}(W),B\rangle=\Sigma_{11}>g(B)$. Therefore $\left\|\mathcal{P}_{\mathcal{T}}(W)\right\|>1\Rightarrow W\notin\partial g(0)$, from which we conclude

[TABLE]

∎