Finite-Sample Performance of Gradient Descent in Logistic Regression with Gaussian Design¶

Conference: NeurIPS2026
arXiv: 2606.21683
Authors: Junren Chen, Arya Mazumdar
Version: arXiv v2, 2026-09-26
Area: Optimization & Theory
Keywords: logistic regression, gradient descent, Gaussian design, finite-sample analysis, separate direction and norm estimation

TL;DR¶

Under well-specified Gaussian logistic regression, the paper combines population curvature analysis with approximate invertibility of the empirical gradient to prove linear convergence of gradient descent to a statistical-error neighborhood, and obtains a sharper high-dimensional error bound by estimating direction and norm separately; large-step acceleration is only local, and the stated near-optimal regime contains a condition discrepancy that must be retained.

Background & Motivation¶

The negative log-likelihood of logistic regression is convex, but decreasing the loss does not directly establish how close a finite-data iterate is to the true parameter. Existing optimization studies often consider bounded covariates, separable data, and the limiting classification direction. Statistical studies analyze the maximum likelihood estimator (MLE), without necessarily identifying how many gradient descent (GD) iterations attain its accuracy. This paper uses standard Gaussian covariates as a tractable benchmark and explicitly tracks sample size, dimension, and the norm of the true parameter.

Write the true parameter norm as \(B=\|\theta^*\|_2\geq1\). A larger \(B\) makes labels resemble deterministic halfspace classification: direction becomes easier to identify, while norm becomes harder. Once labels mostly indicate which side of a hyperplane contains an observation, further increasing parameter magnitude barely changes the observed probabilities. This statistical direction–norm asymmetry also appears in the loss curvature. A single global smoothness constant therefore cannot explain optimization speed, and a stronger signal does not automatically imply more accurate recovery of the full parameter.

The paper first analyzes ordinary GD repeatedly using the same observations, then establishes local large-step acceleration, and finally constructs a separate estimator that is not the MLE. Core idea: analyze anisotropic contraction of the population gradient separately from uniform empirical-gradient perturbations, and calibrate the norm through a one-dimensional projection on independent samples so that the full noise of a high-dimensional vector mean does not enter norm estimation.

Method¶

Overall Architecture¶

The input consists of \(n\) independent and identically distributed observations with covariates \(x_i\sim N(0,I_d)\) and conditional labels \(y_i\mid x_i\sim\operatorname{Bernoulli}(s(x_i^\top\theta^*))\), where \(s(a)=1/(1+e^{-a})\) and the true parameter is fixed rather than random. The main output is a parameter estimate with a Euclidean-error guarantee, not classification accuracy.

There are two distinct algorithmic routes. Algorithm 1 applies full-batch GD to the empirical logistic loss. Its analysis identifies radial and orthogonal curvature in the population Hessian, then transfers contraction to finite-sample iterates through an approximate invertibility condition (AIC). Algorithm 2 first estimates the unit direction through normalized direction iteration, then computes a label-weighted mean on independent held-out observations, projects it onto the estimated direction, inverts a scalar function, and multiplies the estimated norm by the unit direction.

These routes must not be conflated: Algorithm 2's error bound is not an improved guarantee for ordinary GD, nor does it establish that large-step GD attains optimal statistical accuracy. The first two key designs below explain Algorithm 1's analytical mechanisms; the last two describe Algorithm 2's actual estimation procedure. A proof outline is not presented as a network architecture diagram.

Key Designs¶

1. Anisotropic curvature: explaining slow small-step convergence and the locality of large-step acceleration

Gaussian rotational symmetry and Stein's identity express the population gradient as the difference between radial maps evaluated at the candidate and true parameters. Define \(m(\tau)=\mathbb E[s'(\tau g)]\) and \(q(\tau)=\mathbb E[s(\tau g)g]=\tau m(\tau)\), where \(g\sim N(0,1)\). At the true parameter with norm \(B\), the population Hessian has eigenvalue \(q'(B)\asymp B^{-3}\) in the true direction and \(m(B)\asymp B^{-1}\) in orthogonal directions. The slowest direction under a strong signal is therefore the norm direction, rather than every direction being equally difficult.

The small-step analysis integrates the Hessian along the segment joining the candidate and true parameters. When the candidate norm is at most \(2B\), this averaged curvature retains a minimum eigenvalue of order \(B^{-3}\) and a maximum eigenvalue at most \(1/4\), yielding contraction for constant stepsizes. Theorem 1 explicitly uses zero initialization. Although the text describes the result as essentially global, its formal bound cannot simply be restated as holding for arbitrary initialization.

The large-step route linearizes around the true parameter. A stepsize of order \(B\) eliminates orthogonal error faster and improves the radial contraction gap from \(B^{-3}\) to \(B^{-2}\). However, the larger step also amplifies the nonlinear remainder, so the initial error must be at most \(c_0/B\) for that remainder to be absorbed into contraction. Theorem 2 specifies a strict interval, not every stepsize of the same order or an arbitrarily larger one:

\[ \frac{c}{B^2q'(B)}<\eta<\frac{2-cB^{-2}}{m(B)}. \]

Here \(q'(B)=\mathbb E[s'(Bg)g^2]\), and one admissible choice is \(\eta=1/m(B)\). This depends on the unknown signal norm; the paper does not provide a global large-step deployment procedure that dispenses with accurate initialization and norm information.

2. Uniform approximate invertibility: converting statistical perturbations into an iteration-error plateau

Concentration of the gradient at a fixed parameter is insufficient because every GD iterate depends on the same data. The authors therefore control the empirical–population gradient difference uniformly over the relevant parameter region, separating label noise at the true parameter from the empirical process induced by changing the candidate parameter.

The first component uses the conditional variance \(\mathbb E[(s(x_i^\top\theta^*)-y_i)^2\mid x_i]=s'(x_i^\top\theta^*)\). Bernstein's inequality under moment conditions and a sphere-covering argument yield a noise term that decreases as \(B\) increases. For the second component, the authors peel the region into distance shells around the true parameter, construct a covering within each shell, and use the sample covariance to extend the bound from net points to the continuous region. This preserves the fact that the parameter-change term becomes smaller near the truth, instead of replacing it with a coarse bound that remains constant throughout iteration.

Let \(h_{\theta^*}(u)=\nabla L(u)\). Under Theorem 1's sample condition, the analysis yields the uniform one-step relation:

\[ \|u-\theta^*-\eta h_{\theta^*}(u)\|_2 \leq\left(1-\frac{c}{B^3}\right)\|u-\theta^*\|_2 +C\left(\sqrt{\frac{d}{nB}}+\frac dn\right),\qquad \|u\|_2\leq2B. \]

This is the concrete meaning of AIC: the actual descent step does not exactly cancel parameter error, but the remaining error is controlled by contraction plus a random perturbation. Accumulated perturbations are divided by the contraction gap, turning \(\sqrt{d/(nB)}\) into \(\sqrt{B^5d/n}\). The proof also initially produces \(B^3d/n\), which the sample condition absorbs. Although a larger step contracts faster, it also amplifies one-step noise, leaving the final statistical-error plateau unchanged.

The induction in Appendix A.3 additionally ensures that iterates stay inside the region where AIC holds. For the local version, Appendix B.3 requires \(n\gtrsim B^7d\) so that the plateau remains within the small \(c_0/B\) neighborhood. This additional condition is essential for continued applicability of local contraction, not a decorative term.

3. Normalized direction iteration: learning only the more identifiable unit direction

Algorithm 2 does not begin by finding the MLE. Instead, it reuses the direction estimator of Matsumoto and Mazumdar. On the first \(\nu n\) observations, it computes discrepancies between predicted halfspace labels and observed labels, takes a fixed-scale subgradient step, and normalizes the vector to the unit sphere. Each step retains direction while discarding irrelevant length drift:

\[ \widetilde r_{t+1}=\widehat r_t-\frac{\sqrt{2\pi}}{\nu n}\sum_{i=1}^{\nu n}\left(\frac{\operatorname{sign}(x_i^\top\widehat r_t)+1}{2}-y_i\right)x_i, \qquad \widehat r_{t+1}=\frac{\widetilde r_{t+1}}{\|\widetilde r_{t+1}\|_2}. \]

Appendix E.1 explains that this is a subgradient method for a ReLU loss, not simply normalization of the logistic-loss gradient. The direction error is \(\operatorname{polylog}(n)(\sqrt{d/(nB)}+d/n)\). Multiplying it by the true norm contributes \(\sqrt{Bd/n}\) and \(Bd/n\) to full-parameter error. Direction estimation requires no prior knowledge of \(B\), and its first term improves with a stronger signal.

Theorem 3 uses the fixed initialization \(\widehat r_0=e_1\), requires \(T_0\geq\log_2\log_2(n/d)\), and takes \(\nu\in[0.1,0.9]\) with integer \(\nu n\). Although the direction-estimator lemma permits any unit initialization, this note retains the specific choice in the formal theorem for the complete estimator rather than enlarging its claim.

4. Independent projection and inversion: reducing norm calibration to one dimension while retaining saturation risk

The remaining observations produce a label-weighted mean with expectation \(q(B)\theta^*/B\). Estimating \(q(B)\) through the norm of this vector would let the Euclidean noise of a high-dimensional mean introduce dimension into the leading error term. Instead, the authors project onto the estimated unit direction so that norm calibration primarily depends on a scalar mean:

\[ \widehat v=\frac{1}{(1-\nu)n}\sum_{i=\nu n+1}^{n}y_ix_i, \qquad \widehat\lambda=q^{-1}(\langle\widehat r_{T_0},\widehat v\rangle), \qquad \widehat\theta=\widehat\lambda\widehat r_{T_0}. \]

Sample splitting makes the direction estimate independent of the held-out observations, enabling conditional one-dimensional concentration. Projection bias between two unit directions is quadratic in their distance, so direction error does not directly enter norm estimation as a first-order bias. Nevertheless, \(q'(B)\asymp B^{-3}\) means inversion amplifies scalar error. The final norm error still contains \(B^3/\sqrt n\) and \(B^2d/n\); strong-signal difficulty is not eliminated.

Inversion also has a practical domain restriction: \(q^{-1}\) is defined only on \((0,1/\sqrt{2\pi})\). Lemma 9 in Appendix C.1 uses \(1/\sqrt{2\pi}-q(B)\asymp B^{-2}\) and the sample condition to control projection error, ensuring meaningful inversion on the high-probability event. The pseudocode does not specify clipping or an out-of-range fallback valid for every finite sample. An implementation needs a domain check; added clipping must not be attributed to the paper, and a high-probability guarantee must not be described as unconditional validity.

Loss & Training¶

Algorithm 1 uses the objective and update:

\[ L(\theta)=\frac1n\sum_{i=1}^{n}\left[\log(1+e^{x_i^\top\theta})-y_ix_i^\top\theta\right], \qquad \theta_{t+1}=\theta_t-\eta\nabla L(\theta_t). \]

There is no regularization, stochastic mini-batching, or learning-rate schedule here. Theorem 1 bounds every iterate by \(\|\theta_t-\theta^*\|_2\leq(1-c/B^3)^tB+\widetilde C\sqrt{B^5d/n}\). Theorem 2 replaces the decaying term with \((1-c/B^2)^tc_0/B\), keeping the same plateau. Both provide finite-sample parameter guarantees for a neighborhood of the truth, not exact parameter convergence to the truth.

Iteration counts to reach the plateau are respectively \(\widetilde O(B^3)\) and locally \(\widetilde O(B^2)\). These are theoretical orders, not universal stopping rules inferred from the 100, 200, or 400 iterations used in simulations. Algorithm 2's direction iteration and norm inversion likewise cannot be replaced by further optimization of the logistic loss.

Key Experimental Results¶

Main Results¶

The following table first records directly checkable theoretical guarantees. \(C,c,c_0,\widetilde C\) are universal constants, and \(\operatorname{polylog}(n)\) hides logarithmic factors. This is not a table of measured errors.

Result and source	Sample-size condition	Initialization and stepsize/parameters	Error or contraction guarantee	Success probability
Theorem 1, (2)–(3)	\(n\geq C(B^6d\log n+B^6\log B)\)	\(\theta_0=0\), \(\eta\in[0.1,7.9]\)	Contraction factor \(1-c/B^3\); plateau \(\widetilde C\sqrt{B^5d/n}\)	\(1-2e^{-d}\)
Theorem 2, (5)–(7)	\(n\geq C(B^6d\log n+B^7d)\)	Initial error at most \(c_0/B\); strict stepsize interval above, including \(1/m(B)\)	Contraction factor \(1-c/B^2\); same plateau \(\widetilde C\sqrt{B^5d/n}\)	\(1-2e^{-d}\)
Theorem 3, (12)	\(n\geq\operatorname{polylog}(n)(Bd+B^4)\)	\(\widehat r_0=e_1\), \(T_0\geq\log_2\log_2(n/d)\), \(\nu\in[0.1,0.9]\)	\(\operatorname{polylog}(n)(B^3/\sqrt n+\sqrt{Bd/n}+B^2d/n)\)	\(1-n^{-1}\)

The paper also includes numerical simulations; it is not purely theoretical. All results average 50 independent trials, implemented in Matlab R2022a on a laptop with an Intel CPU up to 2.5 GHz and 32 GB RAM. The principal metric is Euclidean parameter error. The cache reports plotted trends without exact error values, so no decimals or percentage improvements are invented from the curves.

Simulation and source	Configuration	Actual comparison	Observation reported in the paper
GD estimation error, Figure 1(a)	\(n\in\{3000,6000,12000,24000\}\); \((d,B)\in\{(200,2),(400,2),(200,3)\}\)	Zero initialization, \(\eta=4\), estimator \(\theta_{100}\)	Log–log trends agree with \(n^{-1/2}\) decay; error increases with dimension and norm
Small-step convergence, Figure 1(b)	\((n,d,B)=(5000,200,4)\); first 200 iterations	Zero initialization, \(\eta=1\) versus \(\eta=4\)	Early log-error curves are approximately straight, supporting the qualitative linear-contraction trend
Local large-step comparison, Figure 1(c)	\((n,d,B)=(80000,100,8)\); first 40 iterations	\(\theta_0=\theta^*+u\), with a random unit direction \(u\); \(\eta=4\) versus \(1/m(8)\approx20.63\)	The larger step has faster initial decay and reaches the statistical plateau earlier; no exact acceleration factor is reported

Figure 1(c) constructs initialization using the known truth, with error exactly 1. Since the theorem does not give a numerical \(c_0\), the experiment cannot be certified to satisfy the formal radius \(c_0/B\). It illustrates the mechanism rather than providing a deployable initialization procedure or numerical certification of the theorem's sample threshold.

Ablation Study¶

There is no conventional module-removal ablation. Appendix E.3 provides algorithm comparisons and sensitivity analyses for sample size and signal norm. The actual settings are retained here without inventing ablation results.

Analysis and source	Shared configuration	Variation	Reported conclusion and boundary
Sample-size sensitivity, Figure 2(a)	\((d,B)=(1000,2)\); Algorithm 2 uses \(T_0=30\); GD uses zero initialization, \(\eta=4\), and 400 iterations	\(n\in\{3000,5000,8000,10000,15000,20000,30000\}\)	For \(n<10000\), Algorithm 2 is more accurate; for \(n>10000\), GD is slightly better; no explicit winner is stated at 10000
Signal-norm sensitivity, Figure 2(b)	\((n,d)=(5000,1000)\); the same two algorithms	\(B\in\{1,2,4,6,8\}\)	Algorithm 2 can substantially outperform GD with high dimension and moderate sample size; no exact errors or ratios are supplied
Theory–implementation difference, §4 and E.3	Both stages of Algorithm 2 use all \(n\) observations	Sample splitting is removed in simulations	The numerical advantage concerns the unsplit implementation; Theorem 3's independent-held-out analysis is not directly a proof for that implementation

Key Findings¶

A larger stepsize improves the speed of reaching the statistical plateau, not its theoretical height. Optimization acceleration and statistical improvement are separate contributions.
Algorithm 2 does not always outperform GD. Figure 2(a) explicitly reports a change in advantage with sample size, preventing a claim of universal empirical superiority over the MLE.
The 400-step GD estimator is used as an approximation to the MLE. No MLE optimality-residual certification is reported, so the paper's qualification that Algorithm 2 possibly also outperforms the MLE is retained.

Highlights & Insights¶

Connecting finite-sample accuracy to the optimization trajectory: AIC controls both one-step error and the final plateau, while uniform concentration handles dependence between iterates and data. This is closer to parameter recovery than loss descent alone.
Curvature reveals both sides of a stronger signal: Orthogonal and radial curvature scale as \(B^{-1}\) and \(B^{-3}\), respectively. Easier classification does not imply easier recovery of parameter length.
Projection is preferable to the norm of a high-dimensional mean: Direction estimation followed by scalar calibration on independent samples reduces the dimensional burden of leading norm noise. The idea may transfer to models with identifiable directions and monotone radial moment maps, but their moment relationships and independence arguments must be verified anew.

Limitations & Future Work¶

Strong model assumptions: Covariates must be standard Gaussian and labels must follow the correctly specified sigmoid model. The results do not automatically cover general sub-Gaussian or correlated features, misspecification, or real datasets. The sixth- and seventh-power dependence on \(B\) in GD's sample conditions is also conservative.
Large steps still require warm starts and tuning: Theorem 2 assumes initialization within \(c_0/B\) of the unknown truth and selects stepsizes using \(B\). Practical adaptive stepsizes and verifiable initialization remain open directions rather than completed solutions here.
Inversion is not an unconditionally safe program: Algorithm 2 does not specify an out-of-range fallback. The theorem handles the saturation endpoint through a high-probability event, while the supplied cache does not explain how the numerical implementation handles exceptional projections.
Remark 5 has an algebraic condition discrepancy: The paper simplifies its three-term bound to \(\widetilde O(\sqrt{Bd/n})\) under the stated condition \(n\gtrsim B^3d+B^5\). Direct comparison of \(B^3/\sqrt n\) with \(\sqrt{Bd/n}\) instead yields \(d\gtrsim B^5\), independently of \(n\); comparing the third term gives \(n\gtrsim B^3d\). Ignoring logarithmic factors, both independent conditions are needed, together with Theorem 3's prerequisites. The discrepancy is retained rather than treating the stated sample-size condition alone as a universal near-optimality conclusion.
Remark 6's sparse extension is not independently verified: Its first norm-related term is \(B/\sqrt n\), inconsistent with \(B^3/\sqrt n\) in the dense theorem and Appendix C; it also mixes sample symbols \(m\) and \(n\). No independent sparse proof is provided in the appendices, so this stronger expression is not treated as an established result.
Simulation dimension notation is questionable: Section 4 places the uniformly sampled true direction on \(\mathbb S^{n-1}\), inconsistent with the parameter belonging to \(\mathbb R^d\). This note does not infer how the actual code samples directions and records results only under the explicit \((n,d,B)\) configurations.

Compared with finite-sample MLE analysis: The existing bound of Chardon, Lerasle, and Mourtada is \(O(\sqrt{B^3d/n})\), sharper than this paper's ordinary-GD guarantee. The main addition here is trajectory and iteration-complexity analysis; a looser upper bound does not establish that GD actually performs statistically worse than the MLE.
Compared with direction recovery: Hsu and Mazumdar identify statistical scales for direction estimation in Gaussian logistic models, while Matsumoto and Mazumdar provide efficient normalized iteration. Algorithm 2 borrows its direction stage from the latter; its added mechanism is projected norm calibration on independent samples and analysis of full-parameter error.
Compared with edge-of-stability GD: Previous large-step studies often use bounded, separable data and parameters tending to infinity. This paper gives a local positive result around a finite true parameter under Gaussian, nonseparable data satisfying its sample conditions. Different assumptions and metrics prevent a unified ranking of convergence rates.
Follow-up work: A footnote in Section 3 of v2 explicitly states that the authors' follow-up paper 2608.17260 substantially improves that section's results. Its full text was not read for this note, so its stronger conclusions are not imported into the present paper.

Rating¶

Novelty: 4/5 — Combines non-asymptotic parameter error for Gaussian GD, contraction speed, and separate direction–norm estimation.
Experimental Thoroughness: 3/5 — Multiple synthetic settings with 50 trials, but no real datasets, exceptional-inversion treatment, or rigorous MLE residual certification.
Writing Quality: 3/5 — Clear main argument and complete proofs, with checkable discrepancies in the near-optimality condition, sparse error term, and simulation dimension notation.
Value: 4/5 — Useful for understanding strong-signal curvature and statistical error; practical extensions remain limited by Gaussian assumptions and conservative sample thresholds.