Efficient Truncated Linear Regression with Unknown Noise Variance

Abstract:
Truncated linear regression is a classical challenge in statistics, wherein a label, y = wx + ε, and its corresponding feature vector, x ∈ R, are only observed if the label falls in some subset S ⊆ R; otherwise the existence of the pair (x, y) is hidden from observation. Linear regression with truncated observations has remained a challenge, in its general form, since the early works of Tobin (1958); Amemiya (1973). When the distribution of the error is normal with known variance, recent work of Daskalakis et al. (2019) provides computationally and statistically efficient estimators of the linear model, w. In this paper, we provide the first computationally and statistically efficient estimators for truncated linear regression when the noise variance is unknown, estimating both the linear model and the variance of the noise. Our estimator is based on an efficient implementation of Projected Stochastic Gradient Descent on the negative log-likelihood of the truncated sample. Importantly, we show that the error of our estimates is asymptotically normal, and we use this to provide explicit confidence regions for our estimates.

heading: 1 Introduction
text: A common challenge facing statistical estimation is the systematic omission of relevant data from the data used to train models. Data omission is a prominent form of dataset bias, which may occur for a variety of reasons, including incorrect experimental design or data collection campaigns, which prevent certain sub-populations from being observed, instrument errors or saturation, which make certain measurements unreliable, societal biases, which may suppress the realization of certain samples, as well as legal or privacy constraints, which might prevent the use of some data. In turn, it is well-understood that failure to account for the systematic omission of data in statistical inference may lead to incorrect models, and this has motivated a long line of research in statistics, econometrics, and a range of other theoretical and applied fields, targeting statistical inference that is robust to missing data; see e.g. Galton (1897); Pearson (1902); Pearson & Lee (1908); Fisher (1931); Tukey (1949); Tobin (1958); Amemiya (1973); Hausman & Wise (1977); Heckman (1979); Breen et al. (1996); Hajivassiliou & McFadden (1998); Mohan et al. (2013); Balakrishnan & Cramer (2014).
In this paper, we revisit the classical challenge of truncated linear regression. As in the standard linear regression setting, we assume that the world produces pairs (𝑥(𝑖), 𝑦(𝑖))𝑖, where each 𝑥(𝑖) ∈ R𝑘
*CD and PS were supported by NSF Awards IIS-1741137, CCF-1617730, and CCF-1901292, by a Simons Investigator Award, by the Simons Collaboration on the Theory of Algorithmic Fairness, by a DSTA grant, 11 and by the DOE PhILMs project (No. DE-AC05-76RL01830)
†MZ was supported by NSF DMS-2023505 (FODSI)
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
is a feature vector drawn from some distribution 𝒫 and its corresponding label 𝑦(𝑖) is linearly related to 𝑥(𝑖) according to 𝑦(𝑖) = 𝑤*𝑇𝑥(𝑖) + 𝜀(𝑖), for both some unknown coefficient vector 𝑤* ∈ R𝑘 and random noise 𝜀(𝑖) ∼ 𝒩 (0, 𝜎*2), where 𝜎*2 is unknown and 𝜀(𝑖) is 𝑖.𝑖.𝑑 over all 𝑖. The difference to the standard setting is that we do not get to observe all pairs (𝑥(𝑖), 𝑦(𝑖)) that the world produces. Rather, a pair (𝑥(𝑖), 𝑦(𝑖)) is only observed if 𝑦(𝑖) ∈ 𝑆, where 𝑆 ⊂ R is a fixed “observation set,” that is known or given to us via oracle access, i.e. we can query an oracle about whether or not some point belongs to 𝑆. If 𝑦(𝑖) /∈ 𝑆, then (𝑥(𝑖), 𝑦(𝑖)) is removed from observation. Our goal is to estimate 𝑤* and 𝜎*2 from the pairs that survived removal, called “truncated samples.”
Overview of our results. While both truncated and censored linear regression are fundamental problems, which commonly arise in practice and have been studied in the literature since at least the 1950s (Tobin, 1958; Amemiya, 1973), there are no known statistical inference methods that are computationally efficient and enjoy statistical rates that scale with the dimension 𝑘 of the problem. Recent work of Daskalakis et al. (2019) has obtained such methods in the simpler setting where the noise variance 𝜎*2 is known. Our goal here is to extend those results in two important ways. First, we study settings where 𝜎*2 is unknown, providing computationally and statistically efficient estimators for the joint estimation of 𝑤* and 𝜎*2. Moreover, we establish the asymptotic normality of our estimators, and use this to derive confidence regions for our estimators. To achieve these stronger results we make some additional assumptions compared to Daskalakis et al. (2019). Namely, we assume that the covariates are sampled from some prior distributions, and we also make a stronger assumption about the survival probability of our model as we discuss below. Other works that consider learning and regression problems with censored or truncated samples are Plevrakis (2021); Fotakis et al. (2020); Moitra et al. (2021); Ilyas et al. (2020).
We next provide informal statements of our results, postponing formal statements to Sections 3 and 5. It is important to note that, as shown in prior work (Daskalakis et al., 2018, 2019), for the model parameters to be identifiable from truncated samples, we need to place some assumptions that limit how aggressively the data is truncated. The assumptions we make are listed below, and stated more precisely in Section 2. They combine normalization assumptions (Assumption 1.1), assumptions that are needed even without truncation (Assumption 1.3), and assumptions that are needed for identifiability purposes in the presence of truncation (Assumption 1.2).
Assumption 1.1. We assume that ‖𝑤*‖22 ≤ 𝛽 and ⃦⃦ 𝑥(𝑖) ⃦⃦2 2 ≤ 1 for all 𝑖. Assumption 1.2. For all (𝑥(𝑖), 𝑦(𝑖)) in the truncated dataset, conditioning on 𝑥(𝑖), the probability that 𝑤*𝑇𝑥(𝑖) + 𝜀 ∈ 𝑆, with respect to the randomness in 𝜀 ∼ 𝒩 (0, 𝜎*2), is larger than some 𝑎. Assumption 1.3. If (𝑥(𝑖), 𝑦(𝑖))𝑛𝑖=1 is the truncated dataset, we assume that the feature covariance matrix, 𝑋 = 1𝑛 ∑︀𝑛 𝑖=1 𝑥 (𝑖)𝑥(𝑖) 𝑇 , satisfies: E𝑥∼𝒫 [𝑋] = E𝑥∼𝒫 [𝑥𝑥𝑇 ] ⪰ 𝑏 · 𝐼 for some 𝑏.
With these assumptions we prove our main results for the joint estimation of 𝑤* and 𝜎*2, listed as Informal Theorems 1.4–1.5 below, and respectively Theorems 3.2 and 5.1 in later sections. Informal Theorem 1.4 (Estimation). Suppose we are given 𝑛 samples (𝑥(𝑖), 𝑦(𝑖))𝑛𝑖=1 from the truncated linear regression model with parameters 𝑤*, 𝜎*2, and suppose Assumptions 1.1–1.3 hold. Then running Projected Stochastic Gradient Descent (Algorithm 1) on the conditional negative log-likelihood of the sample, using the projection set described in Definition 3.1, we obtain with probability at least 2/3 (which can be boosted to as large a constant as desirable) estimations ?̂?, ?̂?2
such that ‖?̂? − 𝑤*‖2 + ⃒⃒ ?̂?2 − 𝜎*2 ⃒⃒ ≤ 𝜀, as long as 𝑛 ≥ 𝐶1 log 4(1/𝜀) 𝜀4 or 𝑛 ≥ 𝐶2 𝜀2 (depending on the learning rate), where 𝐶1 depends polynomially on 𝛽, 1/𝑎 and 1/𝑏, and 𝐶2 depends linearly in 1/𝑏 and exponentially in 𝛽/𝑎, where 𝛽, 𝑎 and 𝑏 are respectively the constants in Assumptions 1.1, 1.2 and 1.3. Informal Theorem 1.5 (Asymptotic Normality). Let ?̂?, ?̂?2 be the estimates of Informal Theorem 1.4. There exists a matrix 𝐴 that depends on ?̂?, ?̂?2 (and is explicitly provided in Theorem 5.1) such that the random vector √ 𝑛( (︀ ?̂? ?̂?2 )︀ − (︀ 𝑤* 𝜎*2 )︀ ) converges asymptotically to a normal distribution with zero mean and covariance matrix 𝐴. Hence the 1−𝛼 confidence region of (?̂?, ?̂?2) can be computed using the matrix 𝐴 and the quantiles of the chi-squared distribution.
Techniques. Our estimates are computed by running Projected Stochastic Gradient Descent (PSGD) on the population conditional negative log-likelihood of the truncated sample; a convex function of
the parameters, after an appropriate reparameterization. We provide two analyses of PSGD resulting in two different bounds that are more effective in a different range of parameters, as we discuss next.
Our first analysis hinges on establishing a lower bound on the convexity of the conditional negative log-likelihood function. We can’t however, expect that a strong bound holds over the whole parameter domain. For this reason, we add a projection step to our method, identifying a projection set wherein we show that the true parameters lie and wherein the conditional negative log-likelihood is strongly convex. This approach is similar to that of Daskalakis et al. (2019), except that the unknown noise variance introduces many challenges in identifying an appropriate projection set. The main result bounding the strong convexity of the conditional negative log likelihood function within the projection set is stated as Theorem 4.1. The main disadvantage of this analysis is that, when the variance is unknown, the bound on the strong convexity of the conditional negative log-likelihood becomes small and the resulting dependence of the sample size on the parameters 𝛽 and 1/𝑎 of the problem is exponential. On the other hand, the dependence of the sample size on 1/𝑏 is linear and its dependence on the required estimation accuracy is Θ(1/𝜀2), which is optimal. Our other analysis follows an approach that has not been pursued in prior work, and in particular does not rely on a global lower bound on the strong convexity of the conditional negative log-likelihood. Instead, we only bound the strong convexity at the optimum and use an upper bound on the smoothness of the objective function to prove the convergence of PSGD in terms of function value. Again, it is impossible to prove a bound on the smoothness over the whole domain, because large noise variance can lead to very bad smoothness bounds. Thus we use the same projection set as above and show that, within this set, we can get an upper bound on the smoothness that depends polynomially on all parameters, 𝛽, 1/𝑎, 1/𝑏, of the problem. The main drawback of this method is that the dependence on the accuracy is Θ(log4(1/𝜀)/𝜀4). It is worth noting that obtaining an algorithm whose sample complexity has an optimal Θ(1/𝜀2) dependence on the estimation accuracy, and also scales polynomially in 1/𝑎 is a fundamental bottleneck, which is already present in the much simpler special case of our problem of estimating the variance of a 1-dimensional truncated Gaussian, studied in Daskalakis et al. (2018).3
Finally, to help build intuition around the intricacies in estimating truncated linear regression models with unknown noise variance, let’s consider why natural approaches of reducing this problem to estimating truncated linear regression models with known noise variance – which can be done using the algorithm of Daskalakis et al. (2019) – fail. A simple such reduction might compute the empirical variance 𝜎20 of the truncated sample, and plug that into the algorithm of Daskalakis et al. (2019) to estimate the linear model, pretending that 𝜎20 is the true noise variance. A more sophisticated approach might grid over candidate noise variances, plug each of these candidates into the algorithm of Daskalakis et al. (2019) to estimate a linear model, and then select the (linear model, noise variance) pair that is most consistent, e.g. the pair whose noise variance is closest to the empirical variance of the residuals between the linear model’s predictions and the observed labels. Both of these approaches fail, as suggested by the example below. Suppose that the ground truth is 𝑦 = 𝑤*𝑥 + 𝑤*0 + 𝜀, where 𝑥 is one-dimensional, (𝑤*, 𝑤*0) = (0, 0), and the noise variance is 𝜎
*2 = 1. Suppose also that we create a truncated dataset by, repeatedly, sampling 𝑥𝑖 from 𝒩 (0, 1), generating 𝑦𝑖 according to the ground truth, and truncating the resulting pair (𝑥𝑖, 𝑦𝑖) unless 𝑦𝑖 ≥ 𝜏 , where 𝜏 is some threshold. It is clear that, if we guess the true noise variance to be the empirical one, then as 𝜏 gets bigger, the empirical variance gets smaller, thus the linear model that would be estimated by treating the empirical variance as the true noise variance would have intercept that is close to 𝜏 rather than 0. So this naive approach would fail. Let us consider the more sophisticated one that grids over candidate noise variances. In particular, let us compare what happens when (1) we guess the noise variance to be 𝜎21 = 𝜎
*2, and run the algorithm of Daskalakis et al. (2019) with this guess to get an estimate (?̂?, ?̂?0); and (2) when we guess the noise variance to be 𝜎22 ≡ 𝜎20 (the empirical variance) and run the algorithm of Daskalakis et al. (2019) with this guess to get an estimate (?̂?′, ?̂?′0). It is easy to see that, as 𝜏 gets large, ⃒⃒⃒ 𝜎*2− 1𝑁 ∑︀ 𝑖(𝑦𝑖− ?̂?𝑥𝑖− ?̂?0)2 ⃒⃒⃒ ≫ ⃒⃒⃒ 𝜎20 − 1𝑁 ∑︀ 𝑖(𝑦𝑖− ?̂?′𝑥𝑖− ?̂?′0)2 ⃒⃒⃒ (here, 𝑁 is the number of samples. This holds even if 𝑁 goes to infinity). Thus the selection criterion proposed above will fail to prefer (1) over (2), even though the guess of the noise variance under (1) is the correct one. Our approach in this paper amounts to using a different criterion, namely selecting the pair with the largest conditional log-likelihood. But analyzing when and how well this works, as well as obtaining an efficient algorithm to find the best pair is the main contribution of this work.
3While Daskalakis et al. (2018) does not explicitly analyze the dependence of the sample complexity on 𝛼; by tracing through their bounds, one can see that the dependence of the sample complexity on 𝛼 is exponential.
Roadmap. We present the truncated linear regression setting and our assumptions in Section 2. We state and prove our main estimation results Section 3 and 4, and prove the asymptotic normality of our estimators in Section 5. Finally, in Section 6 we assess the performance of our methods.

heading: 2 Models and assumptions
text: We study the truncated linear regression setting studied in Daskalakis et al. (2019), with the additional complexity that the variance of the noise model is unknown. In particular, let 𝑆 ⊂ R be a measurable subset of the real line, to which oracle access is provided. Assume we have 𝑁 = 3𝑛 truncated samples (𝑥(𝑖), 𝑦(𝑖))𝑁𝑖=1 (we will split the samples into three, see 4.5), each generated according to the following procedure, for some unknown 𝑤* ∈ R𝑘 and 𝜎* ∈ R that are fixed across different 𝑖’s:
1. sample 𝑥(𝑖) from some distribution 𝒫0 with support supp(𝒫0) ⊆ R𝑘; 2. generate 𝑦(𝑖) according to the following model, where 𝜀 ∼ 𝒩 (0, 𝜎*2);
𝑦(𝑖) = 𝑤*𝑇𝑥(𝑖) + 𝜀 (2.1)
3. if 𝑦(𝑖) ∈ 𝑆 then return (𝑥(𝑖), 𝑦(𝑖)) as the 𝑖th sample, otherwise repeat from step 1. After truncation, denote the distribution of 𝑥 to be 𝒫 with the supp(𝒫) = supp(𝒫0). The sampled 𝑥 is still 𝑖.𝑖.𝑑, so we could only consider the distribution 𝒫 . As already established in prior work (Daskalakis et al., 2018, 2019; Kontonis et al., 2019), some assumptions must be placed on 𝑆 and its interaction with the unknown parameters (𝑤*, 𝜎*2) and the 𝑥(𝑖)’s in order for the parameters to be identifiable. We state our assumptions after a useful definition. Definition 2.1 (Survival Probability). Let 𝑆 ⊆ R be measurable, 𝑥,𝑤 ∈ R𝑘, and 𝜎 ∈ R. We define the survival probability 𝛼(𝑤, 𝜎, 𝑥, 𝑆) as 𝛼(𝑤, 𝜎, 𝑥, 𝑆) = P {𝑌 ∈ 𝑆} ; where 𝑌 ∼ 𝒩 (𝑤𝑇𝑥, 𝜎2). When 𝑆 is clear from context we may omit 𝑆 from the arguments of 𝛼 and simply write 𝛼(𝑤, 𝜎, 𝑥).
We are now ready to state our assumptions. Our first assumption gives an upper bound on the norm of the true parameters and the feature vectors.
Assumption 2.2 (Normalization/Data Generation). Assume that ‖𝑤*‖22 ≤ 𝛽, and ‖𝑥‖ 2 2 ≤ 1 for all 𝑥 ∈ supp(𝒫). We also assume oracle access for set 𝑆, i.e. we have an oracle which determines, for any 𝑦, whether 𝑦 ∈ 𝑆.
Our next assumption is the survival probability lower bound following Daskalakis et al. (2018, 2019). Assumption 2.3 (Constant Survival Probability). For all 𝑥 ∈ supp(𝒫), we have 𝛼(𝑤*, 𝜎*, 𝑥) > 𝑎.
Our final assumption is about the spectrum of the covariance matrix of the features and is classical in some other linear regression settings, such as in Daskalakis et al. (2019).
Assumption 2.4 (Thickness of Feature Covariance). Let 𝑋 = 1𝑛 ∑︀𝑛 𝑖=1 𝑥 (𝑖)𝑥(𝑖) 𝑇 be the 𝑘 × 𝑘 feature vector covariance matrix. We assume that, for some positive real number 𝑏, we have E𝑥∼𝒫 [𝑋] = E𝑥∼𝒫 [𝑥𝑥𝑇 ] ⪰ 𝑏 · 𝐼 . Here we abuse the notation that the (𝑥(𝑖), 𝑦(𝑖))𝑛𝑖=1 are samples used for PSGD.

heading: 3 Main Result
text: When the variance parameter 𝜎*2 is unknown, the conditional log-likelihood function of the plain linear regression becomes non-concave, and is therefore not straightforward to optimize. For this reason, we follow the common practice of reparameterizing the problem with respect to the parameters 𝑣 = 𝑤/𝜎2 and 𝜆 = 1/𝜎2. The algorithm that we use is Projected SGD on the population negative log-likelihood function ℓ̄(𝑣, 𝜆) as presented in Algorithm 1. We formally define the loss function ℓ̄ in Section 4.1.
Our goal is to apply the above algorithm to the conditional negative log-likelihood function of the truncated linear regression model. To execute the above algorithm in polynomial time, we need to solve the following algorithmic problems.
(a) initial feasible point: compute an initial feasible point in some projection set 𝐷,
Algorithm 1 Projected SGD on ℓ̄
Input: (𝑥(𝑖), 𝑦(𝑖)) for 𝑖 ∈ [𝑛], learning rates 𝜂𝑡 1: (𝑤0, 𝜎20)← OLS estimates using other 𝑛 samples; (𝑞, 𝜆(0))← (𝑤0/𝜎20 , 1/𝜎20)
◁ (a) Initialization part 1 2: In terms of 𝜎20 define convex set 𝐷 over pairs (𝑣, 𝜆) as described in Definition 3.1 3: 𝑣(0) ← argmin𝑣:(𝑣,𝜆(0))∈𝐷‖𝑣 − 𝑞‖2 ◁ (a) Initialization part 2 4: for 𝑡 = 1, . . . , 𝑛 do 5: Sample gradient 𝑢(𝑡) such that E [︀ 𝑢(𝑡)|(𝑣, 𝜆)(𝑡−1) ]︀ = ∇ℓ̄(𝑣(𝑡−1), 𝜆(𝑡−1)) ◁ (b) Unbiased Gradient Estimate, sampling method see Section B.1 6: (𝑣(𝑡), 𝜆(𝑡))← (𝑣(𝑡−1), 𝜆(𝑡−1))− 𝜂𝑡 · 𝑢(𝑡) 7: (𝑣(𝑡), 𝜆(𝑡))← arg min(𝑣,𝜆)∈𝐷 ⃦⃦ (𝑣, 𝜆)− (𝑣(𝑡), 𝜆(𝑡)) ⃦⃦ 2
◁ (c) Efficient Projection 8: end for 9: (𝑣, ?̂?)← (︀ 𝑣(𝑛), 𝜆(𝑛) )︀ ; return (?̂?, ?̂?2)← (︁ 𝑣/?̂?, 1/?̂? )︁
(b) unbiased gradient estimation: sample an unbiased estimate of∇ℓ̄(𝑣(𝑡−1), 𝜆(𝑡−1)), (c) efficient projection: design an algorithm to project to the set 𝐷.
Solving the computational problems (a)-(c) is the first step in the proof of our main result. First we define the appropriate projection set to use in our algorithm. Definition 3.1 (Projection Set). We define the projection set to be the (𝑣, 𝜆) satisfying
𝐷 = {︂ (𝑣, 𝜆) = (︂ 𝑤
𝜎2 ,
1
𝜎2
)︂ ∈ R𝑘 × R ⃒⃒⃒ 1 8(5 + 2 log(1/𝑎)) ≤ 𝜎 2 𝜎20 ≤ 96 𝑎2 , ‖𝑤‖22 ≤ 𝛽 }︂ where 𝜎20 is the estimated variance from solving the classical linear regression problem (using 𝑛 samples) on the data ignoring truncation.
This projection set contains two parts: the first part that restricts the weights and the second part that restricts the variance. For the second part, we use the variance predicted using ordinary least squares.
Now we are ready to formally state our main theorem about the estimation of the parameters 𝑤 and 𝜎 from truncated samples under the Assumptions 2.2- 2.4.
Theorem 3.2. Let (𝑥(1), 𝑦(1)), · · · , (𝑥(𝑛), 𝑦(𝑛)) be 𝑛 samples from the linear regression model (2.1) with ground truth 𝑤*, 𝜎*2. If Assumptions 2.2- 2.4 hold, then we can instantiate the learning rates used in Algorithm 1 such that, with success probability ≥ 2/3, the output estimates ?̂?, ?̂?2 satisfy either one of the following bounds:
‖?̂? − 𝑤*‖2 + ⃒⃒⃒ ?̂?2 − 𝜎*2 ⃒⃒⃒ ≤ poly(𝜎 *·𝛽 𝑎·𝑏 , 1 𝑎·𝜎* ) · log(𝑛)
𝑛1/4 (3.1)
‖?̂? − 𝑤*‖2 + ⃒⃒⃒ ?̂?2 − 𝜎*2 ⃒⃒⃒ ≤
poly(𝜎*, 1𝑏 ) · exp (︁ poly( 𝛽𝑎·𝜎* ) )︁
√ 𝑛
(3.2)
for some polynomials that we make explicit in the detailed proof of the theorem (Section B.7). Remark 3.3. It is easy to see that the probability of success in Theorem 3.2 can be boosted to 1− 𝛿 with an additional cost of order log(1/𝛿) in the rates. We can achieve this using a folklore boosting technique for parameter estimation problems. In particular, we can run the algorithm of Theorem 3.2 log(1/𝛿) times independently and then we can pick the estimate that is 𝜀 close to at least the half of the rest of the estimates, where 𝜀 is the target estimation error. It is easy to see that with this boosting technique we can get error at most 2𝜀 with probability of success at least 1− 𝛿.
To prove Theorem 3.2 we provide two different analyses of the Algorithm 1. The first analysis uses only the fact that the stochastic gradients of the population conditional negative log-likelihood function have bounded second moment and is based on the following theorem. Theorem 3.4 (Theorem 2 of Shamir & Zhang (2013)). Let ℓ̄ be a convex function with a minimizer (𝑣*, 𝜆*) and suppose that there exists a constant 𝜌 such that
(i) bounded step variance : E [︁⃦⃦ 𝑢(𝑡) ⃦⃦2 2 ]︁ ≤ 𝜌2 where 𝑢(𝑡) is the sampled step variance in
algorithm 1,
(ii) bounded domain: max(𝑣,𝜆)∈𝐷‖(𝑣 − 𝑣*, 𝜆− 𝜆*)‖ 2 2 ≤ 𝜌2,
(iii) feasibility of optimal: the minimizer (𝑣*, 𝜆*) ∈ 𝐷
then if we apply Algorithm 1 on ℓ̄ with learning rates 𝜂𝑡 = 𝑐/ √ 𝑡 where 𝑐 is an absolute constant then for every 𝑡 ∈ [𝑛] it holds that E [︀ ℓ̄(𝑣(𝑡), 𝜆(𝑡))− ℓ̄(𝑣*, 𝜆*) ]︀ ≤ poly(𝜌)·log(𝑡)√
𝑡 .
For the second analysis of Algorithm 1 we also need to bound the strong convexity of ℓ̄. This is the reason that we get the exponential dependence on some of the parameters of the problem, but we can also get a better consistency rate. For this second analysis we use the following theorem. Theorem 3.5 (Lemma 1 of Rakhlin et al. (2011)). Let ℓ̄ be a convex function with a minimizer (𝑣*, 𝜆*) and suppose that there exist constants 𝜌, 𝜁 such that
(i) bounded step variance: E [︁⃦⃦ 𝑢(𝑡) ⃦⃦2 2 ]︁ ≤ 𝜌2 where 𝑢(𝑡) is the sampled step variance in
algorithm 1
(iii) feasibility of optimal: the minimizer (𝑣*, 𝜆*) ∈ 𝐷
(iv) strong convexity: ℓ̄ is 𝜁-strongly convex
then if we apply Algorithm 1 on ℓ̄ with learning rates 𝜂𝑡 = 1/(𝜁 · 𝑡), then for every 𝑡 ∈ [𝑛] it holds that E [︂⃦⃦⃦ (𝑣(𝑡), 𝜆(𝑡))− (𝑣*, 𝜆*) ⃦⃦⃦2 2 ]︂ ≤ 4𝜌 2 𝜁2 · 𝑡 .
In sum, the main challenges for proving Theorem 3.2 are solving computational problems (a)–(c) and proving mathematical properties (i)–(iv).

heading: 4 Overview of the proof of Theorem 3.2.
text: We provide an overview of the proof of Theorem 3.2, postponing most technical details to the supplementary material. The main steps of the proof are the following.
1. In Section 4.1 we derive the negative population conditional log-likelihood, its gradient, and its Hessian matrix.
2. In Section 4.2 we discuss the computational problems (a) - (c).
3. In Section 4.3 we formally state our results that prove the properties (i), (ii), and (iv).
4. In Section 4.4 we state our results for property (iii), the feasibility of the optimal solution.
5. Finally in Section 4.5 we use all the above results to complete the proof of 3.2.

heading: 4.1 The Negative Population Conditional Log-Likelihood Function of Truncated Linear Regression
text: In this section we define the objective function that we use to apply our PSGD algorithm. This objective function is derived by taking the negative expected value of the log-likelihood of 𝑦 conditional on the value of 𝑥. Given a sample (𝑥, 𝑦), the negative conditional log-likelihood that 𝑦 is a sampled from the truncated linear regression with parameters 𝑤, 𝜎2 given the value of 𝑥 is equal to
ℓ(𝑤, 𝜎2;𝑥, 𝑦) = 1
2𝜎2 (𝑦2 − 2𝑦𝑤𝑇𝑥) + log (︂∫︁ 𝑆 exp(− 1 2𝜎2 (𝑧2 − 2𝑧𝑤𝑇𝑥))d𝑧 )︂
If we reparameterize ℓ with 𝜆 = 1/𝜎2 and 𝑣 = 𝑤/𝜎2, we have
ℓ(𝑣, 𝜆;𝑥, 𝑦) = 1
2 (𝜆𝑦2 − 2𝑦𝑣𝑇𝑥) + log (︂∫︁ 𝑆 exp(−1 2 (𝜆𝑧2 − 2𝑧𝑣𝑇𝑥))d𝑧 )︂
(4.1)
We define the distributions 𝐹𝑥 = 𝒩 (𝑤*𝑇𝑥, 𝜎*2, 𝑆) and 𝑄𝑥 = 𝒩 (𝑤𝑇𝑥, 𝜎2, 𝑆), and the negative population conditional log-likelihood function ℓ̄ as follows
ℓ̄(𝑣, 𝜆) = E 𝑥∼𝒫
[︂ E
𝑦∼𝐹𝑥 [ℓ(𝑣, 𝜆;𝑥, 𝑦)]
]︂ . (4.2)
We can then easily compute
𝜕ℓ̄ 𝜕𝑣 = E 𝑥∼𝒫
[︂ E
𝑧∼𝑄𝑥 [𝑧 · 𝑥]− E 𝑦∼𝐹𝑥 [𝑦 · 𝑥]
]︂ 𝜕ℓ̄
𝜕𝜆 =
1 2 E 𝑥∼𝒫
[︂ E
𝑦∼𝐹𝑥
[︀ 𝑦2 ]︀ − E
𝑧∼𝑄𝑥
[︀ 𝑧2 ]︀]︂ ,
From these we also compute the Hessian matrix which is
H(𝑣, 𝜆) = E 𝑥∼𝒫
[︂(︂ Var𝑧∼𝑄𝑥 [𝑧]𝑥𝑥 𝑇 −Cov𝑧∼𝑄𝑥 [︀ 1 2𝑧 2, 𝑧 ]︀ 𝑥
−Cov𝑧∼𝑄𝑥 [︀ 1 2𝑧 2, 𝑧 ]︀ 𝑥𝑇 Var𝑧∼𝑄𝑥 [ 1 2𝑧 2]
)︂]︂ (4.3)
Here, H is the expectation of covariance matrices, which is positive definite. Thus, ℓ̄ is convex.

heading: 4.2 Computational Problems
text: The computational problems (a) - (c) can be tackled in the following way
(a) We start with computing (𝑤0, 𝜎0) by applying OLS to the truncated data. (b) As we see from the expression of the gradient of ℓ̄, we can compute an unbiased estimate of
the gradient using one of our samples (𝑥(𝑖), 𝑦(𝑖)), our sample access to the distribution 𝒫 of 𝑥’s, and by doing rejection sampling until we find a point with (𝑥, 𝑧) with 𝑧 ∈ 𝑆. To bound the time needed for this step, we need to understand the mass of the survival set 𝑆 for some 𝑥 ∼ 𝒫 and some vector of parameters (𝑣, 𝜆) ∈ 𝐷. See Supplementary Material B.2.
(c) We show in the supplementary material that 𝐷 is a convex set and define an efficient algorithm to project our estimates to 𝐷. It is elaborated in Supplementary Material A.2
For more details on the way we solve the computational problems (a) - (c) we refer to the supplementary material.

heading: 4.3 Bounded Step Variance, Bounded Domain, and Strong Convexity
text: In this section we present the statements that establish the properties (i), (ii), (iv) that we need in order to apply Theorem 3.4 and Theorem 3.5. Theorem 4.1. For every (𝑣, 𝜆) ∈ 𝐷, we have
E [︁ ‖(𝑦 − 𝑧)𝑥‖22 ]︁ ≤ poly(1/𝑎)(1+𝛽2+𝜎40), H(𝑣, 𝜆) ⪰ 𝑏·exp (︂ −poly(1/𝑎)(1 + 𝛽
𝜎20 ) )︂(︂ 𝜎20𝐼 0 0 𝜎40 )︂ where, 𝑥 ∼ 𝒫 , 𝑦 ∼ 𝐹𝑥 and 𝑧 ∼ 𝑄𝑥 and H from (4.3).
The proof of Theorem 4.1 can be found in Supplementary Material B.3. Also, it is easy to see that the diameter of 𝐷 is bounded: ‖𝑣‖2 < ‖𝑤‖2/𝜎2 ≤ 8(5 − 2 log 𝑎) √ 𝛽/𝜎0
2 and |𝜆| < 1/𝜎2 ≤

heading: 8(5− 2 log 𝑎)/𝜎20

heading: 4.4 Feasibility of Optimal Solution
text: We need to calculate the probability for the optimal point to be in the projection set, i.e., 𝑤*, 𝜎* ∈ 𝐷. Lemma 4.2. Under Assumptions 2.3 and 2.2 we have that, P((𝑤*, 𝜎*) ∈ 𝐷) ≥ 13/16.
The proof of this Lemma can be found in Supplementary Material B.4.

heading: 4.5 Proof of Theorem 3.2
text: The algorithm that proves Theorem 3.2 starts with splitting the samples 𝑛 into three parts of equal size in order to avoid dependence. The first part of 𝑛/3 samples is used to compute the OLS estimates
𝑤0, where 𝑤0 is used in the initialization step. The second part of 𝑛/3 samples is used to compute the OLS estimate of 𝜎0, where 𝜎0 is used both to initialize and define the projection set 𝐷, and is used at every step where a projection to 𝐷 is needed. The number of sample to define 𝐷 is poly(𝜎
*,𝛽 𝑎,𝑏 , 1 𝜎* ),
and our 𝑛 itself will reach this bound (see Section B.4.) Finally, the third part of the 𝑛/3 samples is used in our main PSGD algorithm. For the rest of the proof we abuse the notation and we use 𝑛← 𝑛/3 for simplicity. The next thing that we need to decide is the learning rates that we are going to use for the PSGD. We will analyze PSGD under two possible choices of its learning rate schedule, namely 𝜂𝑡 = 𝑐/ √ 𝑡 and 𝜂𝑡 = 1/(𝜁 · 𝑡), and will show that it produces estimates for the ground-truth
parameters that scale as 𝐴1 = poly(
𝜎0·𝛽 𝑎·𝑏 , 1 𝑎·𝜎0
)·log(𝑛) 𝑛1/4 and 𝐴2 = poly(𝜎0)·exp
(︁ poly( 𝛽𝑎·𝜎0 ) )︁
𝑏 √ 𝑛 respectively, for some polynomials that we will make explicit in Section B.7. Analysis for learning schedule 𝜂𝑡 = 𝑐/ √ 𝑡. First, by combining Theorem 3.4 and Theorem 4.1, and noticing the relationship between 𝜎0 and 𝜎* in the projection set, we get the following corollary whose proof can be found in Section B.7. Corollary 4.3. Suppose that (𝑣*, 𝜆*) = (𝑤*/𝜎*2, 1/𝜎*2) ∈ 𝐷 and ℓ̄ is defined as in (4.2). Then, the output (𝑣, ?̂?) of Algorithm 1 satisfies E [︁ ℓ̄(𝑣, ?̂?)− ℓ̄(𝑣*, 𝜆*) ]︁ ≤ poly( 𝜎*·𝛽 𝑎·𝑏 , 1 𝜎* )·log(𝑛)
𝑛1/2 .
Our next step is utilize the above bound to so that our estimates achieve small errors in the parameter space. For this we use the strong convexity of ℓ̄ at the optimum.
Lemma 4.4. The Hessian at the true (𝑣*, 𝜆*) satisfies H(𝑣*, 𝜆*) ⪰ 𝑏×poly(𝑎)poly(1/𝑎,𝛽/𝜎*2)
(︂ 𝜎*2𝐼 0
0 𝜎*4
)︂ .
The proof of Lemma 4.4 can be found in Supplementary Material B.3.1. We also utilize the following lemma to calculate an upper bound on ℓ̄(𝑣*, 𝜆*): Lemma 4.5. (Lemma 3 in Daskalakis et al. (2019)) Let 𝑥,𝑤 ∈ R𝑘, 𝜎 > 0, then E𝑦∼𝒩 (𝑤𝑇 𝑥,𝜎2,𝑆)(𝑦− 𝑤𝑇𝑥)2 ≤ (4− 2 log𝛼(𝑤, 𝜎, 𝑥))𝜎2.
By the above lemma, we have:
ℓ̄(𝑣*, 𝜆*) = E 𝑥∈𝒫 E 𝑦
[︃ (𝑦 − 𝑤*𝑇𝑥)2
2𝜎*2 + log𝛼(𝑤*, 𝜎*, 𝑥)
]︃ ≤ 4− 2 log min
𝑥∈𝒫 𝛼(𝑤*, 𝜎*, 𝑥) ≤ 4− 2 log 𝑎
Applying Markov’s inequality we have P (︁ ℓ̄(𝑣, ?̂?)− ℓ̄(𝑣*, 𝜆*) > 7E [︁ ℓ̄(𝑣, ?̂?)− ℓ̄(𝑣*, 𝜆*) ]︁)︁ < 1
7 . (4.4)
The inequality implies that, with at least 6/7 probability, the actual difference is at most 7 times its expectation. Also, the probability that (𝑣*, 𝜆*) ∈ 𝐷 is at least 13/16 by Lemma 4.2. Combining these two, we have at least 13/16− 1/7 > 2/3 probability to hold:
ℓ̄(𝑣, ?̂?)− ℓ̄(𝑣*, 𝜆*) < poly(𝜎
*·𝛽 𝑎𝑏 , 1 𝜎* ) · log(𝑛)
𝑛1/2 (4.5)
Now applying Lemma 4.4, since 𝑣, ?̂? converge to 𝑣*, 𝜆*, when they have small distance, we have
‖𝑣 − 𝑣*‖22 + (?̂?− 𝜆 *)2 ≤ 2
𝛿
(︁ ℓ̄(𝑣, ?̂?)− ℓ̄(𝑣*, 𝜆*) )︁ ≤ poly(𝜎 *·𝛽 𝑎𝑏 , 1 𝜎* ) · log(𝑛)
𝛿 · 𝑛1/2 (4.6) where H(𝑣*, 𝜆*) ⪰ 𝛿𝐼 , and 𝛿−1 = poly (︁
𝜎*·𝛽 𝑎·𝑏 , 1 𝜎*
)︁ .
Finally, doing some simple calculations for the reparameterization —see Section B.7 in supplementary material— we have that
‖?̂? − 𝑤*‖2+ ⃒⃒⃒ ?̂?2 − 𝜎*2 ⃒⃒⃒ ≤ √︂ 2 (︁ ‖?̂? − 𝑤*‖22 + ⃒⃒ ?̂?2 − 𝜎*2 ⃒⃒2)︁ ≤ poly(1/𝑎, 𝜎*) ⃦⃦⃦(𝑣, ?̂?)− (𝑣*, 𝜆*)⃦⃦⃦ 2 . (4.7) Combining inequalities (4.4), (4.5), (4.6), (4.7), inequality (3.1) of Theorem 3.2 follows.
Analysis for learning schedule 𝜂𝑡 = 1/(𝜁𝑡). We combine Theorem 3.5 and Theorem 4.1 to get the following corollary whose proof can be found in Section B.6.
Corollary 4.6. Suppose the true parameters (𝑣*, 𝜆*) = (𝑤*/𝜎*2, 1/𝜎*2) ∈ 𝐷 and ℓ̄ is defined as equation (4.2). Then, for 𝑣, ?̂? being the output of Algorithm 1, it holds that
E [︂⃦⃦⃦ (𝑣, ?̂?)− (𝑣*, 𝜆*) ⃦⃦⃦2 2 ]︂ ≤ poly (︂ 𝛽, 𝜎* + 1 𝜎* , 1 𝑎 )︂ × exp (︂ poly(1/𝑎) (︂ 1 + 𝛽 𝜎*2 )︂)︂ 1 𝑏2𝑛 .
Similarly —details can be found in Section B.7 of the supplementary material— to those used to prove inequality (3.1) above, we get that with probability ≥ 2/3 it holds that ⃦⃦⃦ (𝑣, ?̂?)− (𝑣*, 𝜆*) ⃦⃦⃦2 2 ≤
poly(𝜎*)·exp (︁ poly( 𝛽 𝑎·(𝜎*) ) )︁
𝑏2·𝑛 . Then, we can prove (3.2) using the definition of the projection set 𝐷.

heading: 5 Inference and Confidence Regions
text: We now prove the asymptotic normality of our estimates (Theorem 5.1), which we use to obtain their confidence regions (Corollary 5.2). For the proofs we refer to the Section B.8.
Theorem 5.1. Let (?̂?, ?̂?2) be the output of Algorithm 1 with 𝜂𝑡 = 1/(𝜁 · 𝑡). Then, the vector√ 𝑛𝑆−1/2( (︀ ?̂? ?̂?2 )︀ )− (︀ 𝑤* 𝜎*2 )︀ ), asymptotically converges to 𝒩 (0, 𝐼𝑘+1), where
- 𝑆 = 1𝜁𝑅 𝑇 Σ𝑅, 𝑅 =
(︂ 1/?̂?2𝐼 −?̂?/?̂?4 0 −1/?̂?4 )︂ ,
- Σ is the matrix that satisfies (︁ Ĥ− 𝜁2𝐼 )︁ Σ + Σ (︁ Ĥ− 𝜁2𝐼 )︁ = Γ(𝑣, ?̂?) = Γ (︀ ?̂? ?̂?2 , 1 ?̂?2 )︀ ,
- Ĥ is the empirical estimation of the Hessian matrix defined by (4.3) at the point (𝑣, ?̂?), - Γ(𝑣, 𝜆) = 1𝑛 ∑︀𝑛 𝑖=1 𝜕𝑙(𝑣,𝜆;𝑥(𝑖),𝑦(𝑖)) 𝜕(𝑣,𝜆) 𝜕𝑙(𝑣,𝜆;𝑥(𝑖),𝑦(𝑖)) 𝜕(𝑣,𝜆) 𝑇 .
Corollary 5.2. Let 𝑞𝛼 be the 1 − 𝛼 quantile of distribution 𝜒2𝑘+1. Then the 1 − 𝛼 asymptotic confidence region of (𝑤, 𝜎2) is (︂(︂
𝑤 𝜎2 )︂ − (︂ ?̂? ?̂?2 )︂)︂𝑇 𝑅𝑇 Σ𝑅 (︂(︂ 𝑤 𝜎2 )︂ − (︂ ?̂? ?̂?2 )︂)︂ ≤ 𝑞𝛼𝜁·𝑛 .
6 Synthetic Experiments
(a)
(b)
vector of ones, sample 𝑛 (𝑥(𝑖), 𝑦(𝑖)) pairs, where 𝑥(𝑖) 𝑖.𝑖.𝑑.∼ 𝒩 (0, 𝐼) and 𝜀(𝑖) 𝑖.𝑖.𝑑∼ 𝒩 (0, 1), and right truncate at zero. We generate our data once. However, for each experiment, we randomly select a 𝑘 samples that exceed zero, and run our experiment. We vary 𝑘 over the interval [10, 5000]. Our results show that the error bound asymptotically approaches zero at a rate that reaffirms our theoretical analysis.

heading: 7 Semi-Synthetic Experiment
text: Here we show results for an additional experiment that we conducted on semi-synthetic data. For this experiment, we used the Aldrin (2004) dataset. The Aldrin (2004) dataset was originally collected by the Norwegian Public Roads Administration for a study of air pollution at a road in Oslo, Norway. The dataset consists of 500 observations. Interestingly enough, it is common for environmental data to be truncated because of problems in reliably measuring low concentrations.
For this experiment, we apply left truncation to the dataset, varying a truncation parameter 𝐶 over the interval [1.0, 4.0]. We run our procedure a total of 10 times and retain the trial that has the smallest gradient as our algorithm’s prediction. We conduct the experiment with the same hyperparameters as for the synthetic data experiments. Below, we report our results.
In Figure 3 (a), we show the 𝐿2 distance between the model parameters of each method and the ground-truth, which we take to be the OLS estimation without truncation. In this figure, we show that our method performs very comparably to Croissant & Zeileis (2018). We also note that, at high values of the truncation, both methods exhibit non-monotonicity in their estimation error, which is presumably due to model misspecification and the truncation being too aggressive.
In Figure 3 (b), we report the 𝐿1 distance between the predicted noise variance and the ground-truth noise variance. Again, our method performs very comparably to Croissant & Zeileis (2018), and we observe non-monotonicity at high values of the truncation.

heading: Acknowledgements
text: CD and PS were supported by NSF Awards CCF-1901292, DMS-2022448 (FODSI) and DMS2134108, by a Simons Investigator Award, by the Simons Collaboration on the Theory of Algorithmic Fairness, by a DSTA grant, and by the DOE PhILMs project (No. DE-AC05-76RL01830). MZ was supported by NSF DMS-2023505 (FODSI).

Review summaries:
# Summary  The paper studies truncated linear regression when the additive error follows a zero-mean normal distribution with unknown variance.  This setting generalizes the work of Daskalakis-Gouleakis-Tzamos-Zampetakis-2019, who assumed the knowledge of the variance of error. This work provides two statistical rates for the joint estimation of the regressor (w) and the noise variance (\sigma) in Theorem 3.2. The rates in Theorem 3.2 are achieved using a SGD-without-replacement algorithm on the negative log likelihood after a suitable re-parameterization, which is also used in Daskalakis-Gouleakis-Tzamos-Zampetakis-2019.  

In this paper, the authors investigate truncated linear regression, i.e., inferring the linear model parameters when only a subset of the data can be observed. The primary contribution is that they provide error bounds on the estimators for model parameters when the noise variance is unknown and show that there exists a polynomial time algorithm to find the estimate (noise variance included).   They validate their results on synthetic data. 

The paper extends the study of Daskalakis et al. (2019) into the setting where the noise variance is unknown, and gives two theoretical contributions. First, the paper proves the convergence of the projected SGD algorithm when applied to minimizing the population negative log-likelihood function. Second, the paper provides the asymptotic normality of their estimator, from which the confidence regions for the estimator can be computed.