Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
72,609 characters · 19 sections · 52 citation commands
Inference for Low-Rank Models
\address{Department of Economics, MIT, Cambridge, MA 02139} \email{[email removed]}
\address{Booth School of Business, University of Chicago, Chicago, IL 60637} \email{[email removed]}
\address{Department of Economics, Rutgers University, New Brunswick, NJ 08901} \email{[email removed]}
\address{Department of Economics, Brandeis University, 415 South St, Waltham, MA 02453 } \email{[email removed]}
Key words: nuclear norm penalization, singular value thresholding, sample splitting
\onehalfspacing
We study inference for linear low-rank models: $$ Y= X \circ \Theta + \mathcal E $$ where $(Y,X,\Theta,\mathcal E)$ are $n\times p$ matrices with both $n,p\to\infty$ and $\circ$ denotes the matrix element-wise product. We observe data $(X,Y)$, and $\mathcal E$ represents unobserved statistical noise. The model parameter is the matrix coefficient $\Theta.$ We assume $\Theta$ follows an approximate spiked low rank model: $\Theta$ can be well-approximated by a low-rank matrix whose rank $J$ is either fixed or grows slowly compared to $n,p$ and whose largest $J$ singular values diverge with $(n,p)$. Our main goal is performing statistical inference on both sparse and dense linear combinations of elements of $\Theta$.
Under the approximate spiked low rank model structure, nuclear-norm regularization provides a natural benchmark approach to estimating $\Theta$. There is a substantial literature that studies rates of convergence of nuclear-norm penalized estimators; see, e.g., koltchinskii2011nuclear and negahban2011estimation for prominent examples. Providing results in low-rank models that allow for ready construction of inferential objects such as confidence intervals has been a topic in the more recent literature. For example, xia2019statistical and chen2019inference study inference in settings where the matrix parameter of interest has an exact low-rank structure with fixed rank and elements of $X$ are i.i.d. copies from an unknown distribution.
We contribute to this literature by establishing asymptotic normality for low-rank estimators. Our method starts with an initial estimator of $\Theta$ obtained using nuclear-norm regularization from which we extract the right singular vectors. We then treat the extracted singular vectors as data and obtain estimates of the left singular vectors and updated estimates of the right singular vectors by applying additional least squares steps. The final estimator is then the product of the estimated left and right singular vectors. We make use of a rotation argument to show that, in terms of estimating the space of the singular vectors, the regularization bias of the first step nuclear-norm penalized estimation aligns with the space spanned by the true singular vectors. Thus, the regularization bias is effectively eliminated by the additional least squares steps.
We prove that our estimator for linear functionals of the low-rank matrix is asymptotically normal. We also establish the semiparametric efficiency bound and show that our estimator attains the efficiency bound. The notion of semiparametric efficiency in the presence of high-dimensional nuisance parameters is adopted from jankova2018semiparametric. Our result is novel relative to jankova2018semiparametric because they deal with sparse models while we look at linear combinations of a low-rank matrix.
Our conditions allow possible strong dependence within $X$, which is useful in many contexts. For example, in the matrix completion context, we can accommodate persistence in observed and missing entries rather than relying on independent missingness. In our treatment effects example, allowing strong dependence allows us to consider scenarios where units are first in the control state for a period of time and then enter the treated state and remain there until the end of the sample period.
We rely on two key technical conditions in establishing asymptotic normality of our proposed estimator. We first assume $\Theta$ has spiked singular values (SSV), which requires that the non-zero singular values are large. This condition ensures that the rank of $\Theta$ can be consistently estimated and that the singular vectors are estimated sufficiently well for use in Stage 2 of the procedure. In the inference context, to make entrywise inference for the low-rank matrix, the SSV condition on the singular values seem necessary.
The second condition relates to incoherent singular vectors as defined in, e.g., candes2009exact,candes2010matrix, keshavan2010matrix; and chen2019noisy. This condition requires that the signals on the singular vectors should be approximately evenly distributed across their entries. Under the incoherence condition, we use a “rotation" argument to show that our approach provides asymptotically unbiased estimates of the eigenvector space. We note that incoherence does rule out the setting of “sparse PCA,” which needs a separate treatment and often requires explicit debiasing steps as in, e.g., jankova2021biased.
We note that the SSV and incoherence conditions are strong and are often absent in the literature when probability bounds are derived. However, probability bounds for the Frobenius risk in general cannot imply the asymptotic distribution of estimators. In particular, one of the objects of interest in this paper is to make inference on sparse linear combinations of rows (or columns) of $\Theta$, including elementwise inference. Recent developments for perturbation bounds of entrywise eigen-analysis require SSV to make entrywise inference for a low-rank matrix; see, e.g., abbe2020entrywise.
To further illustrate that both conditions seem necessary for good performance, we provide new minimax theory on convergence rates without them. These results verify that it is impossible to guarantee entrywise consistency without SSV or incoherence and show that, for dense linear combinations, the optimal rates one can achieve without these conditions are potentially much worse than those available under them. Finally, as these minimax rates do show that inference for dense linear combinations may proceed without SSV or incoherence, we provide an alternative inference for dense functionals when these conditions are relaxed.
The literature
Low-rank regression has been extensively studied in the literature. Much of this work focuses on deriving sharp deviation bounds for low-rank estimators; see, for instance, recht2011simpler; gross2011recovering; rohde2011estimation; koltchinskii2011nuclear; dray2015principal; zhu2019high; candes2010matrix; hastie2015matrix; keshavan2010matrix; and sun2012calibrated. As with xia2019statistical and chen2019inference, our paper complements the literature by providing asymptotic distributional results. A key difference between our work and xia2019statistical and chen2019inference is that we adopt a different approach that does not rely on explicit debiasing steps to achieve asymptotic normality. Rather, we rely on the fact that $\Theta$ is a “product parameter” obtained by multiplication of left and right singular vectors to verify that our procedure produces sufficiently regular estimators for asymptotic normality to hold without explicit debiasing.
Our paper is related to chernozhukov2018inference which considers inference in linear panel data models with multivariate coefficient matrices that admit factor structures. There are several important differences between the two papers. Because chernozhukov2018inference consider estimation of multiple matrix parameters, they employ a complicated orthogonalization step to deal with the fact that regularization bias in any of the matrix parameters spills over and impacts estimation of all other matrix parameters. chernozhukov2018inference also rely on strong conditions on regressors while our conditions allow the regressor in our model to be strongly persistent. This generalization allows us to handle matrix completion problems with “systematic missingness." We also provide several new optimality results to establish semiparametric efficiency and minimaxity. Finally, we explicitly allow for the low-rank structure to be an approximation by accounting for approximation errors and allowing the rank of the approximating low-rank structure to increase with the sample size. Admitting these characteristics broadens the applicability of the method. For example, many nonparametric models cannot be formulated as exact low-rank models with fixed rank but can be approximated with low-rank models with slowly growing rank.
Throughout the paper, we denote the maximum and minimum singular values of a matrix $A$ as $\psi_{\max}(A)$ and $\psi_{\min}(A)$. We use $\psi_j(A)$ to denote the $j^{\text{th}}$ largest singular value of $A$. We use $\|A\|_F$, $\|A\|$ and $\|A\|_{(n)}=\sum_{k=1}^{\min\{n,p\}}\psi_k(A)$ to respectively denote the matrix Frobenius norm, operator norm, and nuclear norm. We let $\|A\|_{\max}=\max_{ij}|(A)_{ij}|$ be the elementwise norm. Let $\mathsf{vec}(A)$ denote the vector that stacks the columns of $A$. For two stochastic sequences, we write $a_n\asymp b_n$ if $a_n=O_P(b_n)$ and $b_n=O_P(a_n)$, which means $a_n/b_n=O_P(1)$. Finally, $a\vee b$ means $\max(a,b)$.
Consider the following model
where we observe data $(y_{ij}, x_{ij})$ and $\varepsilon_{ij}$ is the noise term. Let $(Y, X, \Theta, \mathcal E)$ denote the $n\times p$ matrices of $(y_{ij}, x_{ij}, \theta_{ij}, \varepsilon_{ij})$. Then the matrix form of (ref) is $$ Y= X \circ \Theta +\mathcal E $$ where $\circ$ denotes the matrix element-wise product. The goal is to make inference about linear combinations of elements of $\Theta$. Throughout the paper, we impose that $\Theta$ and its associated singular values/vectors are random. Suppose $\Theta$ can be decomposed as
where $\Theta_0$ and $R$ are $n\times p$ matrices satisfying the following conditions:
(i) $\Theta_0$ is a rank $J$ matrix where $J$ is either bounded or grows slowly compared to $(n,p)$. In addition, the nonzero singular values of $\Theta_0$ are “spiked": $$ \psi_1(\Theta_0)\geq...\geq\psi_J(\Theta_0)\geq \psi_{np},\quad \psi_j(\Theta_0)=0 \ \forall \ j>J $$ for some sequence $\psi_{np}\to\infty$.
(ii) $R$ is the low-rank approximation error whose entries $r_{ij}$ satisfy $$ \max_{i, j} |r_{ij}| \leq O_P( r_{np}) $$ for some sequence $r_{np}\to 0$.
(iii) Let $U_0=[u_1,...,u_n]'$ and $V_0=[v_1,...,v_p]' $ respectively denote the $n \times J$ and $p \times J$ matrices that collect the left singular vectors and right singular vectors of $\Theta_0$ corresponding to the nonzero singular values. We assume incoherent singular-vectors: $$ \max_{j\leq p}\|v_j\| =O_P(\sqrt{Jp^{-1}}),\quad \max_{i\leq n}\|u_i\| =O_P(\sqrt{Jn^{-1}}). $$
Given the approximate low-rank structure of $\Theta$, a natural estimation strategy is nuclear-norm penalized optimization:
where $\mathcal A=\{\Theta: \|\Theta\|_{\max}\leq M\}$ and $\nu$ is a tuning parameter. Imposing the $\max$-norm constraint with a large constant $M>0$ helps stabilize the solution; see, e.g., klopp2014noisy. Nuclear-norm penalized regression is natural as the solution is easy to compute. Another option would be to explicitly penalize the matrix rank; however, obtaining the solution to the rank penalized problem is in general difficult unless all elements of $X$ are equal to one. Statistical properties of ((ref)), focusing on the minimax rate for $\|\widetilde\Theta-\Theta\|_F$, have been well-studied in the literature; see, e.g., koltchinskii2011nuclear and negahban2011estimation. It is also well-known that the singular values of $\widetilde\Theta$ suffer from shrinkage biases, and dealing suitably with this shrinkage bias is an important component of our inferential procedure.
Finally, we assume that $J$, the rank of the low-rank component $\Theta_0$, is known for simplicity. For instance, in the treatment effect study where the parameter matrix is approximated by a low-rank structure via a sieve representation, the rank equals the sieve dimension which could be pre-specified. In cases where rank is unknown, it can be consistently estimated. For example, one can apply the singular value thresholding method where the cut-off value for “large singular values” can be chosen to dominate the noise level; see, e.g., onatski2010determining and fan2022estimating.
Let the singular value decomposition of $\Theta_0$ be $$ \Theta_0= U_0 D_0 V_0':= \Gamma_0V_0',\quad \Gamma_0:= U_0D_0. $$ Here, $D_0$ is the $J\times J$ diagonal matrix containing the nonzero singular values of $\Theta_0$, and $U_0$ and $V_0$ are respectively the $n\times J$ left singular vector matrix of $\Theta_0$ and $p\times J$ right singular vector matrix of $\Theta_0$ corresponding to the non-zero singular values. Let $\gamma_{0,i}'$ for $i = 1,...,n$ denote the rows of $\Gamma_0$, and let $v_{0,j}'$ for $j = 1,...,p$ denote the rows of $V_0$.
We only iterate least squares once in Step 2. The least squares steps following the use of nuclear-norm penalized estimation are analogous to approaches in the sparse regression setting that rely on refitting the least squares using selected regressors in a first step, such as post-lasso, e.g. belloni2013least. The motivation is similar in wanting to alleviate shrinkage biases induced in the initial penalized estimation step. In addition, we split the sample $ \{1,...,n\}\backslash \{i\}=\mathcal I \cup \mathcal I^c $ so that $i$ is excluded from both subsamples. Splitting in this way ensures that the $\varepsilon_{ij}$ for the $i$ of interest are independent of observations in both subsamples assuming independence across $i$.
Stage 2, which involves two least squares estimation steps, is the essential stage to alleviating shrinkage bias. It starts with treating $\widetilde V$ from the penalized regression as observed data. A key ingredient of the analysis is to establish that this step produces an approximately unbiased estimator $\widehat\Gamma$, which then allows construction of a well-behaved estimator $\widehat\Theta$ in the final step. Given its importance, we provide the intuition for this step in Section (ref).
We make use of a “rotation" argument and the structure of the low-rank matrix parameter to prove that eigenspace estimation is approximately unbiased if singular vectors are incoherent. Before turning to the matrix parameter setting, we introduce the main idea in the context of estimating a scalar parameter that is itself a product of two parameters.
Consider the problem of estimating a scalar parameter $\theta$ that can be written as the product of another two scalar parameters: $$ \theta= \gamma \beta,\quad \gamma, \beta\in\mathbb R. $$ Suppose some initial estimate $\widetilde\beta$ can be obtained for $\beta$, which is consistent but may have first-order bias. In addition, suppose that $\theta$ can be identified as the unique minimizer of a population loss function: $$ \theta =\arg\min_{a\in \mathcal A} Q(a) $$ where $Q(\cdot)$ is the loss function and $\mathcal A$ is the parameter space. Let $Q_n(\cdot)$ denote the sample version of $Q$ and suppose both $Q_n$ and $Q$ are twice continuously differentiable. Let $\dot{Q}_n(a)=\frac{d}{da}Q_n(a)$, $\dot{Q}(a)=\frac{d}{da}Q(a)$, and $\ddot{Q}(a)= \frac{d^2}{da^2}Q(a)$. Let $(\gamma,\beta,\theta)$ represent the true values of the parameters.
We consider an iterative procedure to estimate $\theta$ that mimics the approach we propose in the matrix parameter setting:
(i) Obtain $\widehat\gamma=\arg\min_{\gamma} Q_n(\gamma \widetilde\beta)$
(ii) Obtain $\widehat\beta=\arg\min_{\beta} Q_n(\widehat\gamma \beta)$
(iii) Set $\widehat\theta= \widehat\gamma\widehat\beta $.
In step (i), standard analysis based on Taylor expansion leads to
where $G=-\partial^2_{\gamma,\gamma}Q_n(\gamma\beta)$. The first term in the expansion is the score which leads to asymptotic normality in usual cases. The second term reflects the effect of the initial estimate $\widetilde\beta$.
{In general, the second term will lead to poor performance of $\widehat\gamma$ if the initial estimator $\widetilde\beta$ is ill-behaved. One approach, dating back to at least Neyman59, is to rely on estimation strategies make use of appropriately “orthogonalized” scores. This property would correspond to basing estimation on an objective function that satisfied $\partial^2_{\gamma,\beta} Q(\gamma\beta)=0$ at the population level in the present case.} See, e.g., CHS:AnnRev for a review of such approaches.
The fact that the “product parameter” $\theta$, rather than $\gamma$ itself, is the object of interest allows a new argument in this paper. The key is that the loss function depends on $\theta$ only through the product of $(\gamma,\beta)$. It is straightforward to verify that $$\partial_{\gamma,\beta}^2 Q(\gamma\beta) =\gamma\ddot{Q}(\theta) \beta + \underbrace{\dot{Q}(\theta)}_{\text{score}=0} =\gamma\ddot{Q}(\theta) \beta. $$ Substituting this expression for $\partial_{\gamma,\beta}^2 Q(\gamma\beta)$ into ((ref)) then produces $$ \widehat\gamma -\gamma =G^{-1} \beta\dot{Q}_n(\theta) +G^{-1} \gamma \ddot{Q}(\theta)\beta(\widetilde\beta-\beta) +o(|\widehat\gamma-\gamma|). $$ An important observation is that the second term $G^{-1} \gamma \ddot{Q}(\theta)\beta(\widetilde\beta-\beta)$ is proportional to $\gamma$. We can move it to the left-hand-side of the expansion for $\widehat\gamma$ to obtain $$ \widehat\gamma - H\gamma = G^{-1} \beta\dot{Q}_n(\theta) +o(|\widehat\gamma-\gamma|) $$ for $H:= 1+ G^{-1} \ddot{Q}(\theta)\beta(\widetilde\beta-\beta).$ Hence, $\widehat\gamma$ estimates a “rotated” version of $\gamma$ with no first-order bias. As such, in the sense of estimating the “space” of $\gamma$, the effect $\widetilde\beta-\beta$ is negligible as it is “absorbed" by the rotation matrix. In addition, $H$ is asymptotically invertible since $H\to^P 1$.
Moving on to step (ii), it is clear that $\widehat\beta$ estimated in this step will be an approximately unbiased estimator for $H^{-1}\beta$. The rotation matrices will then cancel in estimating the parameter of interest: $$\widehat\theta:=\widehat\gamma\widehat\beta =\gamma HH^{-1}\beta +o_P(1) = \theta + o_P(1). $$ After appropriate scaling, the leading term hidden in the $o_P(1)$ in the final expression will also be asymptotically normal. It is this cancellation of rotation matrices that underlies our “rotation-unbiasedness". Furthermore, in models where $\sqrt{n}$-consistency is attainable, $\sqrt{n}(\widehat\theta-\theta)$ is asymptotically normal as long as the initial estimator satisfies $|\widetilde\beta-\beta|=o_P(n^{-1/4})$.
The intuition of “rotation-unbiasedness" as described above has also been observed previously in the literature. keshavan2010matrix studied local geometric properties in Grassmann manifold and related optimization algorithms. sun2016guaranteed examined the local geometry of the loss $f(\Gamma, V)=\|Y-\Gamma V'\|_F^2$ in the matrix completion context. Our observation aligns with theirs, but we use this observation in the context of estimation bias. In our setting, the geometry of product-parameter $\gamma\beta$ ensures that the effect of first-step estimation error $\widetilde\beta-\beta$ is aligned with the space of the true $\gamma$. This alignment results in our ability to establish asymptotic normality of our final estimator without relying on any additional debiasing schemes beyond the use of a single set of least squares steps in Step 2 of our algorithm.
In the low-rank inference context, recall that $\Theta_0=\Gamma_0 V_0',$ which is the product of two parameters. Related to the simple example in the previous section, we think about $V_0$ as $\beta$ and use the singular vectors $\widetilde V$ extracted from the nuclear-norm regularized estimator as its initial estimate.
Write $\widehat\Gamma=(\widehat\gamma_1,...,\widehat\gamma_n)'$ and $\widetilde V=(\widetilde v_1,...,\widetilde v_p)'$. Then for each $i\leq n$, $$ \widehat\gamma_i=\arg\min_{\gamma} Q_{i}(\gamma, \widetilde V),\quad Q_{i}(\gamma, \widetilde V):=\sum_{j=1}^p[y_{ij}- x_{ij}\cdot \gamma'\widetilde v_j]^2. $$ Then for some $J\times J$ matrix $G^{-1}$, Taylor expansion leads to $$ \widehat\gamma_i-\gamma_i = G^{-1} \partial_\gamma Q_{i}(\gamma_i, V_0) + \frac{\partial^2Q_{i}(\gamma_i, V_0)}{\partial\gamma\partial\mathsf{vec}(V)} \mathsf{vec}(\widetilde V-V_0) + \text{higher order terms}. $$ The leading term $G^{-1} \partial_\gamma Q_{i}(\gamma_i, V)$ is asymptotically normal if $V_0$ is incoherent. The second term satisfies
for some rotation matrix $H_1$ and higher order term $\Delta_{i}$.
The term $H_1 \gamma_i$ is a rotated version of $\gamma_i$. Defining $H:= I+ H_1$ and moving $H_1\gamma_i$ to the left side then yields the matrix form expansion: $$ \widehat\Gamma-\Gamma_0 H= \partial_\Gamma Q_{p}(\Gamma, V_0) G^{-1} + \text{higher order terms} $$ where $\partial_\Gamma Q_{p}(\Gamma, V_0) $ is an $n\times J$ matrix whose $i^{\text{th}}$ row is the transpose of $\partial_\gamma Q_{i}(\gamma_i, V_0)$. Following the logic outlined in Section (ref), we have that the follow-up estimator $\widehat V$ will recover an appropriately rotated version of $V$ to cancel with $H$. Consequently, $$ \widehat\Theta = \widehat\Gamma\widehat V' \approx \Gamma_0 HH^{-1}V_0' =\Gamma_0V_0'=\Theta_0. $$
It will then follow that $\widehat\Theta$ is approximately unbiased with sampling distribution that can be approximated by a centered Gaussian distribution. As in the simpler scalar case, the key feature we take advantage of is that we only need the estimated $V_0$ to have the same span as the actual $V_0$ if our goal is inference about $\Theta$ or the space spanned by the singular vectors.
Our argument for demonstrating that the higher-order term, $\Delta_i$ is asymptotically negligible relies on sample splitting. The structure of $\Delta_i$ is $$ \Delta_{i}= B \sum_{j=1}^p(\widetilde v_j-v_j)\varepsilon_{ij} x_{ij}. $$ for some matrix $B$.
For a fixed $i$, let $\mathcal I\subset\{1,...,n\}\backslash i $ be a subset of unit indexes that does not include $i$; and let $$ \mathcal D_{\mathcal I} =\{(y_{kj}, x_{kj}): k\in \mathcal I, j\leq p\}. $$ Our approach uses only data $\mathcal D_{\mathcal I}$, rather than making use of the full data set, for the initial nuclear-norm penalized regression from which we extract singular vectors for the subsequent OLS rotation-debiasing step. Maintaining independence across $i$, estimation errors in the initial estimator of the singular vectors are then independent of variables indexed by $i$ because $i\notin \mathcal I$. Assuming $ \varepsilon_{ij} $ is independent across subjects $i=1,..., n$, we then have that $\varepsilon_{ij} x_{ij} $ is independent of estimation error in the singular vectors, $\widetilde v_j-v_j$. We can then easily argue that $ \Delta_{i}$ has no impact on the asymptotic distribution of the final estimator.
We now present our main results. In Section (ref), we lay out key conditions and state our result on asymptotic normality. We then provide a brief discussion of semiparametric efficiency in Section (ref) and then highlight the role of the key SSV and incoherence conditions in Section (ref) where we present novel minimax results. Finally, we present an alternative estimation scheme for dense linear combinations in Section (ref).
The goal is to establish inferential theory for the linear functional $\theta_i'g$. Here $\theta_i'$ denotes the $i^{\text{th}}$ row of $\Theta$, and $g=(g_1,...,g_p)'\in\mathbb R^p$ is a vector of weights of interest with non-zero weights collected in $$ \mathcal G=\{j\leq p: g_j\neq 0\}. $$ Inference on a linear combination of a column of $\Theta$ can be carried out similarly by switching the roles of $i$ and $j$. Two examples of $g$ are of particular interest:
Sparse weights: $g$ is a sparse vector with a bounded number of non-zero elements:
$\theta_i'g$ thus corresponds to a linear combination of a small number of elements and may be used when we are particularly interested in just a few components of $\theta_i$. The sparse $g$ scenario includes $g=e_j$ where $e_j=(0,...,0,1,0,...,0)$ is the $j^{\text{th}}$ standard vector for a particular $j$ in which case $\theta_i'g=\theta_{ij}$.
Dense weights: $g$ is a dense vector, in the sense that $|\mathcal G| = O(p)$, but
In this case, $\theta_i'g$ typically represents a weighted average of all components of $\theta_i$ and includes $g=(\frac{1}{p},...,\frac{1}{p})'$ as a special case.
The following assumption formally quantifies the requirement of $g$. Consider the matrix of standardized right singular vectors: $$ \bar V'= \sqrt{p}V_0'. $$
The next assumption restricts the noise data generating process (DGP).
Assumption (ref) (ii) is stated in a way that specifically covers the well-known matrix completion problem: $$ y_{ij}^*= \theta_{ij} + e_{ij} $$ where $y_{ij}^*$ may not be observable, and $x_{ij}$ indicates the observability for each element. Then $\varepsilon_{ij}= e_{ij} x_{ij}$.
The assumption below restricts the DGP of the design variable $x_{ij}$. The restrictions imposed are mild, and the assumption is stated so as to cover a variety of cases. Specifically, conditions (a)-(c) in Assumption (ref) allow for various types of dependence among the $x_{ij}$.
Under Condition (ii).a, we can accommodate both conventional factor models by setting $x_{ij} = 1$ for all $i,j$ as well as conditional empirical factor models, where $x_{ij} = x_j$, with varying coefficients. An example of the latter is an asset pricing model with risk premia that vary across assets and over time where $x_j$ represents the common time-varying market factor.
Condition (ii).b could cover examples of PCA with missing data under heterogeneous missing probabilities as in zhu2019high. In this case, we may take $j$ to represent subjects and $i$ to represent the index of repeated sampling within subject. The condition also accommodates scenarios where $x_{ij}$ represents a treatment indicator where random assignment of subjects $i$ to treatment states occurs independently in each period $j$. Such a structure may approximate some digital experimentation settings.
Condition (iii).c allows for some types of strong dependence in $x_{ij}$ across both $i$ and $j$ but restricts $x_{ij}$ to be binary as would be appropriate in missing data, matrix completion, and treatment assignment settings. In this condition, the set $\mathcal B_i$ represents unit-specific “observation times" for unit $i$; and the set $\bar{\mathcal B}$ is common to all units. The quantity $\max_{i\leq n}\sum_{j=1}^n1\{j\in \bar{\mathcal B} \vartriangle \mathcal B_i\}$ thus measures the difference between the “unit specific" observation times and the “common" observation times. Condition (iii).c requires that these differences should be negligible. Hence, all units should be observed at approximately the same time. For instance, suppose every unit is observed most of the time in the sense that $$ \max_{i\leq n} \sum_{j=1}^p1\{j: x_{ij}=0\} = o_P\left(d_{n,p}\right). $$ Then Condition (iii).c holds with $\bar{\mathcal B}=\{1,...,p\}$.
Next, recall that $v_j$ and $u_i$ are respectively the $j^{\text{th}}$ right singular vector and the $i^{\text{th}}$ left singular vector of $\Theta_0$.
The incoherence condition ensures that information regarding the eigenspace accumulates as the dimension increases and allows us to apply our “rotation" argument to argue that estimating the eigenvector space is asymptotically unbiased. We provide low-level conditions that are sufficient for the incoherence condition in a treatment effects context where the low-rank matrix is formulated using nonparametric sieve representations in equation ((ref)).
The next assumption places restrictions on various moments.
Finally, we present the required conditions on $\psi_{np}$, the signal strength of the non-zero singular values. Recall that $\psi_j(A)$ denotes the $j^{\textnormal{th}}$ largest singular value of $A$. We allow the eigengap to change with $J$, depending on constants $b, d \geq 0$. This generality complicates statement of the condition but is needed to accommodate settings where the rank $J$ is allowed to increase with sample sizes. We provide low-level conditions that are sufficient for the following assumption in the context of a treatment effect example in Lemma (ref).
To estimate the asymptotic variance, we need to preserve the rotation invariance property of the asymptotic variance. We therefore estimate $\sigma_{np}^2$ separately within subsamples and produce the final asymptotic variance estimator by averaging the results across subsamples. We consider the homoskedastic case where $\mathsf{Var}(\varepsilon_{ij}|\Theta, X)= \sigma_j^2 $ for some constant $\sigma_j^2$, $j=1,...,p.$ In this case, standard errors can be estimated as
where $\widehat L_{j,\mathcal I}=\sum_{t\notin\mathcal I} x_{tj}^2\widehat \gamma_t\widehat\gamma_t'$, and $ \widehat B_{\mathcal I}=\sum_{j=1}^p x_{ij} ^2\widetilde v_{j,\mathcal I}\widetilde v_{j,\mathcal I}' $, and $\widehat L_{j,\mathcal I^c}$ and $\widehat B_{\mathcal I^c}$ are defined similarly.
It is interesting to note that $\sigma^2_{np}:=s_{np,1}^2+s_{np,2}^2=O_P(\frac{1}{n}\|g\|^2+\frac{1}{p})$ in the case of fixed $J$. Thus, in this setting, the scaling of the asymptotic variance depends heavily on $\|g\|^2$.
The semiparametric efficiency bound for the case of sparse $g$ was established by chen2019inference (Lemma 2) in matrix completion settings and by iwakura2014asymptotic (Theorem 4.5) in pure factor models. Our asymptotic variance attains these previously established bounds if $e_{ij} $ is i.i.d. homoskedastic Gaussian, so we do not further discuss semiparametric efficiency in the sparse setting.
We now provide a semiparametric efficiency bound in the case of dense $g$ and verify that our estimator achieves this bound. For concreteness, suppose we are interested in $h(\Theta)=\theta_1'g$ where $\theta_1'$ is the first row of $\Theta$ and $g$ is dense. In providing our result, we will allow for a wide range of distributions for $X_{1j}$ while maintaining the assumption that the error term is Gaussian to make calculation tractable.
Specifically, we suppose that $x_{ij}$ follows the distribution $f$ and $e_{ij}\sim N(0,\sigma^{2})$ are independent across $(i,j)$. Let $\mu_f=\mathsf{E} _{f} X_{1j}^2$. Under Assumptions (ref)-(ref), the dominant term in the asymptotic variance is
and we also have $s_{np,1}^2=o_P(s_{np,2}^2)$. Hence, $\widehat\theta_i'g-\theta_i'g=O_P(\|V_0'g\|)$ with asymptotic variance $$ s_{np,1}^{2}+s_{np,2}^{2}=s_{*}^{2}(\Theta,f,\sigma)(1+o_{P}(1)) $$ in this case.
The following result verifies that $s_{*}^{2}(\Theta,f,\sigma)$ matches with the semiparametric efficiency bound. The notion of semiparametric efficiency in the presence of high-dimensional nuisance parameters is adopted from jankova2018semiparametric. The idea is to derive the asymptotic Cram\'er-Rao bound for asymptotically unbiased estimators, and needs to be formally established in the high-dimensional setting. Our result is novel relative to jankova2018semiparametric because they deal with sparse models and our setting has low-rank matrices as the nuisance parameters.
Two key conditions that underlie our main results are the incoherence condition, Assumption (ref), and the spiked singular-value (SSV) condition, Assumption (ref). We demonstrate the role of these conditions by providing minimax theory for estimating $\theta_i'g$ for a sparse or dense $g$ in a simple matrix completion problem where the missing indicators $x_{ij}$ are independent Bernoulli random variables without imposing SSV or incoherence.
Define the following set of low-rank matrices \[ \mathcal{S}=\left\{ A\in\mathbb{R}^{n\times p}:\ \mathsf{rank}(A)\leq J\ {\rm and}\ \max_{1\leq i\leq n}\max_{1\leq j\leq p}|A_{ij}|\leq c_{1}\right\} \] for a constant $c_{1}>0$ and for $J\geq1$. Here, $J$ is allowed to be either a fixed constant or a sequence tending to infinity.
We prove the following result for matrix completion over the space $\mathcal{S}$. Let $y_{ij}=x_{ij}\theta_{ij}+e_{ij}$, where $x_{ij}\sim{\rm Bernoulli}(\rho_{j})$ and $e_{ij}\sim N(0,\sigma_{ij}^{2})$ are independent across $(i,j$). Suppose that there are constants $c_{2},...,c_{6}>0$ such that $\rho_{j}\in(c_{2},1-c_{2})$ and $\sigma_{ij}\in(c_{3},c_{4})$ for any $(i,j)$. Let $\rho=(\rho_{1},...,\rho_{p})'$ and $\sigma=\{\sigma_{ij}\}_{1\leq i\leq n,\,1\leq j\leq p}$. In the theorem below, $T$ represents any measurable function of the data, typically regarded as an “estimator" for $h(\Theta)=\theta_1'g.$
Theorem (ref) gives the minimax rate without SSV and the incoherence condition. It provides a similar intuition to koltchinskii2020efficient. For instance, ((ref)) shows that it is impossible to guarantee entrywise consistency for sparse $g$ in the considered setting without SSV or incoherence.
In addition, Equation ((ref)) implies that the rate $O_P(p^{-1/2})$ is minimax optimal for estimating dense averages in the absence of SSV and incoherence. This rate of convergence is slower than that obtained in Theorem (ref) which makes use of SSV and incoherence. For instance, in the factor model with a finite number of strong factors, Theorem (ref) implies that the rate of convergence can be as fast as $\frac{1}{p}\sum_j\widehat\theta_{ij}-\frac{1}{p}\sum_j \theta_{ij}= O_P\left(\frac{1}{\sqrt{np}}+\frac{1}{p}\right)$.\footnote{This rate holds if the factors have zero mean so that $V_0'g=\frac{1}{p}\sum_{j=1}^pv_j=O_P(p^{-1})$, which is the case for no-intercept factor models. Strictly speaking, this setting was ruled out by Assumption (ref), which requires $\|V_0'g\|\geq cp^{-1/2}$. However, Assumption (ref) is used only for obtaining the asymptotic distribution. This assumption can be relaxed when only the rate of convergence is of interest. }
These minimax results for estimating linear combinations of elements of a low-rank matrix {without SSV and incoherence} are new to the literature. The result closest to ours is koltchinskii2020efficient which provides minimax rates for estimating linear functionals of the eigenvectors of low-rank matrices. They show that the minimax optimal rate can be slow if the SSV condition does not hold. Other results on the minimax bounds for learning an eigenspace can be found in berthet2013optimal, birnbaum2013minimax, and cai2013sparse.
When $g$ is a vector of dense weights, the second minimax result in Theorem (ref) suggests that consistency can be achieved without the SSV and incoherence conditions at the cost of a slower rate of convergence. For completeness, we introduce an alternative estimator that could be used when one does not with to impose these assumptions.
Specifically, suppose $g=(g_1,...,g_p)'\in\mathbb R^p$ is a vector of dense weights as defined in (ref), and we are interested in the functional $h_i(\Theta):=\theta_i'g$. We propose the following estimator in the spirit of inverse probability weighting: $$ \widehat{h_i(\Theta)} = \sum_{j=1}^p\frac{g_j y_{ij}x_{ij}}{\widehat \mu_{j,i}^2},\quad \widehat\mu_{j,i}^2=\frac{1}{n-1}\sum_{k\neq i}x_{kj}^2. $$ Note that this estimator does not require knowing the rank or even that the rank is consistently estimable. It is defined as the weighted average of the $i^{\text{th}}$ row of $Y$ and $X$ with weight proportional to a leave-one-out estimator of the inverse of $\mu_j^2:=\mathsf{E} x_{ij}^2$.
Let $$W_{ij}: =x_{ij}\varepsilon_{ij}+ x_{ij}^2\theta_{ij}.$$
As an important illustration, we show how to apply our framework in a treatment effects setting. Suppose that, for each time $j=1,...,p$ and each unit $i=1,..., n$, there is a pair of potential outcomes
Here $m$ is denotes treatment $(m=1)$ or control ($m=0$) state. In any time period $j$ and for any unit $i$, we observe either $Y_{ij}(1)$ or $Y_{ij}(0)$, but not both, depending on the unit's realized treatment state in that period. The treatment effect depends on time-varying functions $h_{j,m}(.)$ of unit specific state variable $\eta_i$; both $h_{j, m}(\cdot)$ and $\eta_i$ may be unobservable and random. {For clarity, we focus on the scenario where the goal is to perform statistical inference on a long-run treatment effect for a given unit $i$: $$ \tau_i:= \frac{1}{p}\sum_{j=1}^p\nu_{ij} $$ where $\nu_{ij}= h_{j,1}(\eta_i)- h_{j,0}(\eta_i)$ is the treatment effect for unit $i$ at time $j$.}
Define the treatment status indicator $$ x_{ij}(m)=1\{\text{unit $i$ at period $j$ is in state $m$} \}= 1\{Y_{ij}(m) \text{ is observable} \}. $$ Consider the following treatment scenario. Suppose the entire time span $\{1,2,...,p\}$ is divided into two periods, $$ T_0=\{1,...,p_0\} \ \text{and} \ T_1=\{ p_0+1,...,p\}, $$ where both $p_0$ and $p_1:=p-p_0$ are large and both periods are known. We assume
where $d_{n,p_0}$ and $d_{n,p_1}$ are slowly growing sequences defined in ((ref)). That is, each unit is in the control state during most periods in $T_0$, and each unit is the treatment state during most periods in $T_1$. We thus refer to $T_0$ and $T_1$ respectively as the “control period” and the “treatment period”. We refer to this treatment scenario as “systematic treatment,” and note that treatment assignments are strongly dependent in this setting, which results in an important difference from much of the literature on inference in matrix completion settings. In terms of our formal conditions, this scenario corresponds to the case of Assumption (ref) (ii).c.
Let $ \theta_{ij}(m):=h_{j,m}(\eta_i).$ We can then rewrite the model for potential outcomes (ref) as
where $y_{ij}(m) = Y_{ij}(m) x_{ij}(m)$, and $\varepsilon_{ij}(m)=e_{ij}(m)x_{ij}(m).$ Let $\Theta(m)$ denote the $n\times p$ matrix of $(\theta_{ij}(m))_{n\times p}$. As, e.g., previously note by athey2018matrix, it is then clear that recovering elements of $\Theta(m)$ is equivalent to solving a matrix completion problem.
In Section (ref) we provide sufficient conditions to establish that $\Theta(m)$ is an approximate low-rank matrix that satisfies both the SSV and incoherence conditions. Under these conditions, we can then estimate treatment effects by simply applying Algorithm (ref) twice -- once using the data from period $T_0$ and once using the data from period $T_1$.
Step 1: Apply Algorithm (ref) to ((ref)) to estimate $\Theta(0)$.
Step 2: Apply Algorithm (ref) to ((ref)) to estimate $\Theta(1)$.
Step 3: Make inference on the treatment effects from the estimated $\Theta(1)-\Theta(0)$.
Let $\widehat\theta_{ij}(m)$ denote the $(i,j)$ element of the estimated matrix $\Theta(m)$. The average treatment effect estimator is then given by $$ \widehat\tau_i:= \frac{1}{p_1}\sum_{j\in T_1} \widehat\theta_{ij}(1) - \frac{1}{p_0}\sum_{j\in T_0} \widehat\theta_{ij}(0). $$ It is straightforward to extend Theorem (ref) to this context, which leads to the asymptotic distribution of the estimated treatment effects. Formal results are to be presented in Section (ref).
We show that the matrix formed from elements $h_{j, m}(\eta_i)$ can be approximated by a low-rank matrix with slowly growing rank. To aid in focusing on the main idea, we suppress the notation “$m$” throughout this section.
Consider a family of time-varying functions $h_{j}(\cdot)$ of subject-specific latent variables $\eta_i$. Let $\Theta$ be the $n \times p$ matrix obtained by setting the $(i,j)$ element of $\Theta$ to $h_j(\eta_i)$. Suppose $h_j(\cdot)$ has a sieve approximation:
where $\Phi_i:=(\phi_1(\eta_i),...,\phi_J(\eta_i))'\in\mathbb R^J$ is a set of sieve transformations of $\eta_i$ using $\phi_k(\cdot)$ as the basis functions, $\lambda_j=(\lambda_{j,1},...,\lambda_{j,J})'$ is the vector of sieve coefficients for $h_j(\cdot)$, and $r_{ij}$ is the sieve approximation error. Write $\Phi$ as the $n\times J$ matrix of $\Phi_i$, $\Lambda$ as the $p\times J$ matrix of $\lambda_j$, and $R$ as the $n\times p$ matrix of $r_{ij}$. Then the matrix form of ((ref)) is $$ \Theta = \underbrace{\Phi\Lambda' }_{\Theta_0}+ R. $$ Clearly, $\mathsf{rank}(\Theta_0)\leq J$, and there is a rotation matrix $H$ so that columns of $\Lambda H$ are the right singular-vectors of $\Theta_0$. The error $R$ is naturally present as the sieve approximation error which will decrease as more elements are considered in the sieve approximation. It is then natural to consider sequences where $J$ increases slowly with $(n,p)$.
We now illustrate how both the SSV and incoherence conditions can hold in this setting under sensible conditions on the functional space and the sieve bases. Suppose $h_j$ belongs to a H\"{o}lder class: For some $C, \beta,\alpha>0$, $$ \{h: \max_{\gamma_1+...+\gamma_{d}=\beta}\left|\frac{\partial^\beta h(x)}{\partial x_1^{\gamma_1}...\partial x_{d}^{\gamma_{d}} }-\frac{\partial^\beta h(y)}{\partial y_1^{\gamma_1}...\partial y_{d}^{\gamma_{d}} }\right|\leq C\|x-y\|^\alpha, \text{for all }x,y\}. $$ Further suppose that a common basis, such as polynomials or B-splines are considered. We will then have $$ \max_{ij}|r_{ij}|\leq C J^{-a},\quad a= (\beta+\alpha)/\dim(\eta_i), $$ which can be made arbitrarily small for sufficiently smooth functions even if $J$ grows slowly.
Now, suppose there exists a $b \geq 0$ such that $\psi_{J}(\Theta) \leq \psi_1(\Theta) \leq C J^b\psi_J(\Theta)$ for some $C>1$. It is then easy to show that the sequence $\psi_{np}$ can be taken as $$ \psi_{np}\asymp \sqrt{J^{-(2b+1)}\sum_{i=1}^n\sum_{j=1}^p h_j(\eta_i)^2}. $$ We then have that the top $J$ singular values grow at this rate which leads to the SSV condition.
Finally, write $S_\Lambda=\frac{1}{p} \Lambda'\Lambda$, $S_\Phi= \frac{1}{n}\Phi'\Phi$, and $A =S_{\Phi}^{1/2} S_{\Lambda} S_{\Phi}^{1/2}.$ Also let $G_\Phi$ be a $J\times J$ matrix whose columns are the eigenvectors of $A$, and let $T$ be the diagonal matrix of corresponding eigenvalues. Letting $H_\Phi:= S_\Phi^{-1/2} G_\Phi$, it can be verified that $$ \Theta_0\Theta_0' \Phi H_\Phi = pn \Phi H_\Phi T \ \text{and} \ \frac{1}{n} ( \Phi H_\Phi)' \Phi H_\Phi= I. $$ Thus, the columns of $\frac{1}{\sqrt{n}} \Phi H_\Phi$ are the left singular-vectors of $\Theta_0$, and the eigenvalues of $ npA $ equal the first $J$ eigenvalues of $ \Theta_0'\Theta_0$. Similarly, we can define $ H_\Lambda= S_\Lambda^{-1/2} G_\Lambda$ where $G_\Lambda$ is a $J\times J$ matrix whose columns are the eigenvectors of $S_{\Lambda}^{1/2} S_{\Phi} S_{\Lambda}^{1/2}$. Hence, we have
Thus,
It then follows that the incoherence condition holds as long as we can obtain proper upper bounds for $ \max_{j\leq p} \|\lambda_j\|$ and $ \max_{i\leq n} \|\Phi_i\|$. For example, if $\{h_j(\cdot): j\leq p\}$ is further restricted to a Hilbert space with a uniform $L_2$- bound, $$ \max_{j\leq p}\sum_{k=1}^{\infty}\lambda_{j, k}^2<\infty, $$ then $ \max_{j\leq p} \|\lambda_j\|<C$.
We formalize the preceding discussion in the following assumption and lemma.
We now verify the eigengap condition: Let $A= \frac{1}{pn}\Theta\Theta'.$ There are constants $b,d\geq0$ such that
Below we verify the above conditions when the treatment functions are generated from a Gaussian process.
Suppose $\eta_i$ are uniformly generated from $[0,1]$, and functions $h_{j}(\cdot)$ are independently generated from a Gaussian process with covariance kernel $$ K(\eta_1,\eta_2) = \mathsf{Cov}(h_j(\eta_1), h_j(\eta_2)), $$ where $K(\cdot,\cdot)$ is a continuous positive semi-definite kernel function supported on a compact set. In addition, suppose the associated integral operator $$(Tf)(\cdot)=\int K(\cdot, \eta) f(\eta) d\eta$$ is positive semi-definite. Let $\{\bar\phi_k(\cdot)\}$ and $\nu_k\geq0$ be the eigenfunctions and eigenvalues of $T$. Then by Mercer's theorem, $\{\bar\phi_k(\cdot)\}$ is an orthonormal basis so that $K$ has the following representation: $$ K(\eta_1,\eta_2)=\sum_{k=1}^{\infty} \nu_k\bar\phi_k(\eta_1)\bar\phi_k(\eta_2), $$ where the infinite sum can be approximated arbitrarily well by finite truncation $J$ as $J\to\infty.$
Now consider the $n\times n$ matrix $\frac{1}{p}\Theta\Theta'$, whose $(i,l)$ element is $$ \frac{1}{p} \sum_{j=1}^ph_j(\eta_i)h_j(\eta_l) = K(\eta_i,\eta_l) +o_P(1) =\bar\Phi_i' D_{\lambda} \bar\Phi_l +o_P(1) $$ where $\bar\Phi_i'=(\bar\phi_1(\eta_i),...,\bar\phi_J(\eta_i))$ and $D_{\lambda}$ is a diagonal matrix of $(\nu_1,...,\nu_J)$. Also, because the $h_j$ are independently generated from the Gaussian process, the $o_P(1)$ terms are uniform over all elements. Thus, we have an approximate low-rank representation of $\Theta\Theta'$: $$ \Theta\Theta' =\left[\sum_{j}h_j(\eta_i)h_j(\eta_l)\right]_{n\times n}\approx p \bar\Phi D_{\lambda} \bar\Phi' . $$
Because the columns of $\bar\Phi$ are formed from eigenfunctions, its columns are approximately orthonormal bases as eigenvectors of $\Theta\Theta'$. Hence the diagonals of $D_{\lambda}$ are also approximately the top $J$ eigenvalues of $\frac{1}{np}\Theta\Theta'$. This observation heuristically shows that the top eigenvalues of $\frac{1}{np}\Theta\Theta'$ are approximately the same as those of the integral operator $T$ associated with the reproducing kernel function.
Rigorously, we can verify this condition as follows. The conditions of the lemma below are required to hold for both $m\in\{0,1\}$ in our treatment effect setting.
Building on the previous subsections, suppose $h_{j, m}(\eta_i) $ has the following sieve representation: $$ h_{j, m}(\eta_i) = \sum_{k=1}^J\lambda_{j, k, m}\phi_k(\eta_i) + r_{ij}(m),\quad m=0,1. $$ We then have that the matrix $\Theta(m):=(\theta_{ij}(m))_{n\times p_m}$ admits an approximate low-rank structure for each $m\in\{0,1\}$:
where $\Lambda_m$ is the $p\times J$ matrix of $\lambda_{j,k,m}$.
Note that $\widehat\tau_i$ estimating a sensible average treatment effect relies on an additional stability assumption. Define, for $m\in\{0,1\}$, $$ \zeta_{ij}(m):= x_{ij}(m) v_j(m)'\bar B(m)^{-1} \frac{1}{p_m}\sum_{j\in T_m} v_j(m) $$ where $\bar B(m)= \sum_{j\in T_m}x_{ij}(m) v_j(m) v_j(m)'. $ Applying the analysis of Theorem (ref), we have
This expansion yields the asymptotic distribution of $\widehat\tau_i$ under the condition that the second line on the right-hand-side is bounded by $o_P(\min\{p_0, p_1\}^{-1/2})$. That is, we need stability of treatment and control averages in the sense that the average of $\theta_{ij}(0)$ and $\theta_{ij}(1)$ obtained over the respective subsamples does not deviate too far from the infeasible average that would be obtained looking over the entire sample period.
We now illustrate the performance of our inferential approach through a small simulation study in the systematic treatment assignment setting. We report results for $n=p=400$.
To generate data, we first divide the period of observation $\{1,...,p\}$ equally into two periods $T_0$ and $T_1$ each consisting of $p_m = p/2$ observation times. To generate $x_{ij}(m)$, we generate $n_i$ integers $j_1...j_{n_i}$ without replacement to form a set $ A_i(m)=\{j_1,..., j_{n_i}\} \subset T_m$. The number $n_i\leq N_0$ is uniformly generated to be less than a predetermined number $N_0 \in \{p_m^{1/2},p_m^{1/3},p_m^{1/4}\}$. We then set $$ x_{ij}(m)=
. $$ Hence, for each unit $i$, $x_{ij}(m)=1$, \textit{with up to $N_0$ exceptions}, throughout period $T_m$ whose total length is $p_m$. In addition, we generate the noise $\varepsilon_{ij}$ independently across both $(i,j)$ and $\varepsilon_{ij}(m)\sim \mathcal N(0,\sigma_{e}^2)$ for $\sigma_{e}=1$.
One of the key conditions in this scenario is that the treatment effect should be stable in the sense that $\frac{1}{p } \sum_{j=1}^p\theta_{ij}(m)$ can be well approximated by $\frac{1}{p_m} \sum_{j\in T_m}\theta_{ij}(m)$. We thus consider the simplest possible setting where this condition holds by generating time invariant treatment functions: $$ h_{0}(\eta_i) = \sum_{k=1}^{\infty} \frac{|W_{k}|}{k^a}\sin(k\eta_i),\quad h_{1}(\eta_i) = \sum_{k=1}^{\infty} \frac{(|W_{k}|+2)}{k^a}\sin(k\eta_i). $$ Here $\eta_i\sim $ Uniform$[-1,1]$, $W_{k}\sim \mathcal N(0,1)$, and the noise is $e_{ij}\sim \mathcal N(0,1)$. The power parameter $a>1$ quantifies the decay speed of the sieve coefficients.
In terms of implementation of our procedure, we also need $J$ and $\nu$. We do not attempt to infer the rank $J$ from the data. Rather, we look at estimates based on four pre-specified values of the rank: $J=1,...,4$. We set the parameter $\nu$ for the nuclear-norm penalized optimization through a simple plug-in procedure. Specifically, we set
where $\bar Q(W ; q)$ denotes the $q^{\textnormal{th}}$ quantile of a random variable $W$ and $Z$ is an $n\times p_m$ matrix whose elements $z_{ij}$ are generated as $\mathcal N(0,\widehat\sigma_{e}^2)$ independent across $(i,j)$ for some estimated $\widehat\sigma_{e}^2$.\footnote{We set $\widehat\sigma^2_{e}$ by obtaining an initial guess, $\tilde\sigma^2_{e}$, from estimating the simple model $y_{ij} = x_{ij} \theta_i +\sigma_{e}^{-1}u_{ij}$ where $\mathsf{Var}(u_{ij})=1$. We then obtain an initial solution to the nuclear-norm regularized optimization problem with tuning parameter set as in (ref) with $z_{ij} \sim N(0,\tilde\sigma_{e}^2)$. Letting $\widetilde\theta_{ij}$ denote the nuclear-norm regularized estimator obtained with this initial tuning. We then set $\widehat\sigma^2_{e}=\frac{1}{np}\sum_{ij} \widetilde \varepsilon_{ij}^2,$ where $\widetilde \varepsilon_{ij}=y_{ij}-x_{ij}\widetilde\theta_{ij}$.} This choice can be motivated as in belloni2013least and chernozhukov2018inference.
We report simulation coverage probabilities of 95% confidence intervals for $\tau_1$ formed using estimated standard errors based on 1000 simulation replications in Table (ref). Overall, the derived asymptotic distributions seem to provide reasonable approximations to the finite sample distributions under our simulation settings, and the good performance appears quite robust to the choice of $J$ in this simulation.