Diaconis–Ylvisaker prior penalized likelihood
for \(p/n \to \kappa \in (0,1)\) logistic regression
LSE Fellow in Statistics
Department of Statistics
London School of Economics
p.sterzinger@lse.ac.uk
psterzinger.at psterzinger
Greek Stochastics 2026
Zakynthos, Greece
28 July 2026
Joint with
This talk
Sterzinger P (2026+). Proportional-limit asymptotics for Diaconis–Ylvisaker-penalised logistic regression with fitted intercept. ArXiV coming soon.
Related work
Sterzinger P, Kosmidis I (2026). Diaconis–Ylvisaker prior penalized likelihood for \(p / n \to \kappa \in (0, 1)\) logistic regression. Biometrika. 113(2), asag014 DOI: 10.1093/biomet/asag014
Logistic regression
UCI’s Multiple Features
200 handwritten numeral patterns per class, from Dutch utility maps, digitized as binary images
- 76 Fourier coefficients of the character shapes, rotation invariant, per digit
- 64 Karhunen-Loève coefficients per digit
1000 digits for training + 1000 digits for testing
Aim
Explain the contribution of the feature sets in describing the digit “7”
Considerations
Font, digitization
noise, downscaling,
complicate
discrimination
If only rotation
invariant features are
used, “7” and “2” can
look similar
Logistic regression
Data
\(y_1, \ldots, y_n \in \{0, 1\}\)
\(x_1, \ldots, x_n \in \Re^p\)
Model
\(Y_{1}, \ldots, Y_{n}\) conditionally independent with
\[ Y_i \mid {x_i} \sim \mathrm{Bernoulli} \left( \rho'(x_i^\top \beta_0) \right), \quad \rho(\eta) = \log(1 + \exp\{\eta\}) \]
Maximum likelihood
Log-likelihood
\(\displaystyle \ell(\beta; y, X) = \log(L(\beta; y, X)) = \sum_{i = 1}^n \left\{y_i (x_i^\top \beta) - \rho\left(x_i^\top \beta \right) \right\}\)
Maximum likelihood (ML) estimator
\(\hat{\beta} = \underset{\beta \in \Re^p}{\arg \max} \, \ell(\beta; y, X)\)
Classification
\(\phi(x) = \mathbb{1}\{ \rho'(x^\top \hat \beta) > 1/2 \}\)
Predicting “7”
fou + kar vs fou
| Residual Df | Log-likelihood | Df | LR | pvalue |
|---|---|---|---|---|
| 923 | 0.00 | |||
| 859 | 0.00 | 64 | 0.00 | 1 |
| Model | Separation | Infinite estimates |
|---|---|---|
| fou + kar | true | true |
| fou | true | true |
ML estimates do not exist (cf. Albert & Anderson (1984))
Diaconis-Ylvisaker prior
Diaconis-Ylvisaker prior
Prior
\(\displaystyle \log p(\beta; X) = \frac{1 - {\color[]{#E16A86}{\alpha}}}{{\color[]{#E16A86}{\alpha}}} \sum_{i=1}^n \left\{ \rho'\left(x_i^\top {\color[]{#E16A86}{\beta}_P}\right) x_i^\top \beta - \rho\left(x_i^\top \beta \right) \right\} + C\)
Penalized log-likelihood
\(\displaystyle \ell^*(\beta; y, X) = \frac{1}{\alpha} \sum_{i = 1}^n \left\{y_i^* (x_i^\top\beta) - \rho'\left(x_i^\top \beta\right) \right\}, \quad y_i^* = \alpha y_i + (1 - \alpha) \rho'\left(x_i^\top \beta_P\right)\)
Maximum DY prior penalized likelihood
Maximum DY prior penalized likelihood (MDYPL)
\(\hat{\beta}^{\textrm{\small DY}}= \arg \max \, \ell^*(\beta; y, X)\)
MDYPL implementation using ML procedures
MDYPL is ML with pseudo-responses: \(\ell(\beta; y^*, X) / \alpha = \ell^*(\beta; y, X)\)
Existence and uniqueness
\(\hat{\beta}^{\textrm{\small DY}}\) is unique and exists if \(X\) has full column rank (Rigon & Aliverti, 2023, Theorem 1)
Shrinkage
Predicting “7”
MDYPL
Predicting “7”
Penalised likelihood ratio test statistic
| Df | PLR |
|---|---|
| 64 | 63.02 |
\(\displaystyle \qquad \Longrightarrow \Pr\left(\chi^2_{64} > 63.02\right)\approx 0.511\). Fail to reject \(H_0: \beta_{\texttt{kar}} = 0\)
\(p/n \to \kappa \in (0,1)\)
Tall data
FIxed \(p\)-asymptotics: number of parameters is small relative to sample size
Approximate with the asymptotic regime: \(\displaystyle p/n \to 0 \)
Modern data
dim(uci_fou_kar[train,])
[1] 1000 140Number of parameters not negligible relative to the sample size
Approximate with the asymptotic regime: \(\displaystyle p/n \to \kappa \)
Artificial example
\[ \begin{array}{l} n \in \{100, 200, \ldots, 6400 \} \\ \kappa = p/n = 0.2 \\ x_{ij} \stackrel{\text{iid}}{\sim} \mathop{\mathrm{N}}(0, 1/n) \\ \beta_0 =\sqrt{\frac{10}{\kappa}}\, (b_1,\ldots,b_{p/2},0,\ldots,0)^\top, \quad b_j = 1 \end{array} \qquad \Longrightarrow \qquad \gamma^2 = \mathop{\mathrm{var}}(x_i^\top\beta_0) = 5 \]
\[ Y_i \mid {x_i} \stackrel{\text{ind}}{\sim} \mathop{\mathrm{Bernoulli}}\left(\rho'(x_i^\top \beta_0)\right) \]
Performance
Recent developments
Candès & Sur (2020)
sharp phase transition about when the ML estimate does not exist, when \(\eta_i = \theta_0 + x_i^\top \beta_0\), \(x_i \sim \mathop{\mathrm{N}}(0, \Sigma)\), \(p/n \to \kappa \in (0, 1)\), \(\mathop{\mathrm{var}}(\eta_i) \to \gamma^2\)
Sur & Candès (2019), Zhao et al. (2022)
a method, based on approximate message passing, that for \(\theta_0 = 0\), recovers estimation and inferential performance by appropriately rescaling \(\hat\beta\), whenever that exists
Rigon & Aliverti (2023)
study maximum Diaconis-Ylvisaker prior penalized likelihood estimator, which always exists, and report solid empirical performance in high-dimensional setting
Phase transition
Phase transition, \(\hat\beta\) (Candès & Sur, 2020)
Phase transition, \(\hat\beta / v^*\) (Sur & Candès, 2019)
Crossing the phase transition
?
Proportional-limit asymptotics for MDYPL
Setup
Data
\(\displaystyle x_i \sim \mathrm{N}(0_p, p^{-1}I_p), \quad p / n \to \kappa \in (0,1), \quad Y_i = \mathbb{1}\{ \varepsilon_i < \rho'({\color[]{#E16A86}{\theta}_0} + x^\top \beta_0) \}, \quad \varepsilon_i \sim \mathrm{U}[0,1]\)
Signal
\(\displaystyle \Gamma_p = \begin{bmatrix} \frac{\| \beta_0 \|_2^2}{p} & \frac{\langle \beta_0, \beta_P \rangle}{p} \\ \frac{\langle \beta_0, \beta_P \rangle}{p} & \frac{\| \beta_P \|_2^2}{p} \end{bmatrix} \longrightarrow \begin{bmatrix} \gamma^2 &{ \varphi} \\ { \varphi }& { \delta}^2 \end{bmatrix} = \Gamma, \quad \begin{bmatrix} \theta_0 \\ \theta_P \end{bmatrix} \longrightarrow \begin{bmatrix} \theta_0^* \\ \theta_P^* \end{bmatrix}\)
\(\displaystyle \pi_p = \frac{1}{p}\sum_{j = 1}^{p} \delta_{(\beta_{0,j}, \beta_{P,j})} \overset{W_2}{\longrightarrow} \pi_{\{\bar{\beta}_0, \bar{\beta}_P\}}\)
MDYPL
\(\displaystyle Y_i^* = {\color[]{#E16A86}{\alpha}} Y_i + (1-{\color[]{#E16A86}{\alpha}}) \rho'({\color[]{#E16A86}{\theta}_P} + x_i^\top {\color[]{#E16A86}{\beta}_P})\)
CGMT (Thrampoulidis et al., 2018)
PO–AO pair
\[\begin{align} \Phi( G) &= \underset{w\in\mathcal S_{w}}{\min} \underset{u\in\mathcal S_{u}}{\max} \left\{ u^\top {\color[]{#E16A86}{G}} w+\psi_0(w,u) \right\} \tag{PO}\label{eq:cgmt-pair-main} \\ \phi(g,h) &= \underset{w\in\mathcal S_{w}}{\min} \underset{u\in\mathcal S_{u}}{\max} \left\{ \| w \|_2{\color[]{#E16A86}{g}}^\top u + \| u \|_2{\color[]{#E16A86}{h}}^\top w + \psi_0(w,u) \right\} \tag{AO} \end{align}\]
- \(G \in \Re^{n \times p}\), \(g \in \Re^n\), \(h \in \Re^p\) all i.i.d. \(\mathrm{N}(0,1)\)
- \(\psi(w,u)\) convex–concave
- \(\mathcal{S}_u, \mathcal{S}_w\) compact
Value link
\[ \begin{aligned} \Pr(\Phi( G)<c) &\leq 2\Pr(\phi(g,h)\leq c), \quad \Pr(\Phi( G)>c) \leq 2\Pr(\phi(g,h)\geq c) \end{aligned} \]
MDYPL PO–AO pair
\[\begin{align} \Phi(H_2) &= \underset{w}{\min} \underset{\lambda\in\mathcal B_1^n}{\max} % \left\{ -\frac{1}{n\sqrt p}\lambda^\top H_2w +\psi(\lambda) % \right\} \tag{PO} \\[0.8em] \phi(g,h) &= \underset{w\in W_p}{\min} \underset{\lambda\in\mathcal B_1^n}{\max} \, - \frac{1}{n \sqrt p} \left(\| w \|_2g^\top\lambda + \| \lambda \|_2w^\top h\right) + \psi(\lambda) \tag{AO} \end{align}\]
for \(B = [\beta_0, \beta_P]\), and \(\beta \in \Re^p\), \(X = H / \sqrt p\), \(H\beta = H_1 v + H_2 w\) and
\[ \psi(\lambda) = \frac1n \underset{(v,\theta,\eta)\in\mathcal K}{\min} \, % \left\{ 1^\top\rho(\eta)-\eta^\top{Y^*} +\lambda^\top \left( \eta-\theta1-\frac1{\sqrt p}H_1v \right) % \right\} \]
CGMT (Thrampoulidis et al., 2018)
Limiting AO
\[ \begin{aligned} \underset{\substack{ \sigma\geq0\\ \theta\in\Re\\ t>0\\ u \in \operatorname{range}(\Gamma) }}{\min} \underset{r\in(0,1)}{\max} \, % \left\{ \mathop{\mathrm{E}} \left[ m_\rho \left( \xi +\sigma G +\frac{t}{r}Y^*, \frac{t}{r} \right) \right] - \sigma r\sqrt\kappa + \frac{rt}{2} - \mathop{\mathrm{E}}\left[Y^* \left(\xi u\right)\right] - \frac{t}{2r}\mathop{\mathrm{E}}\left[\left(Y^*\right)^2\right] \,, % \right\}\,, \end{aligned} \] \(\xi = \theta+z^\top u, z \sim \mathrm{N}(0_2,I_2), [Q_1, Q_2] = \Gamma^{1/2} z, Y^* = \alpha \mathbb{1} \{ \epsilon < \rho'(\theta_0^* + Q_1) \} + (1 - \alpha) \rho'(\theta_P^* + Q_2)\)
FOCs
\[ \begin{aligned} \mathop{\mathrm{E}} \left[ \rho' \left( \textrm{prox}_{\lambda\rho}\left(\xi + \sigma G +\lambda Y^*\right) \right) - Y^* \right] &= 0, \\ P_{\operatorname{range}(\Gamma)} \mathop{\mathrm{E}} \left[ z \left\{ \rho' \left( \textrm{prox}_{\lambda\rho}\left(\xi + \sigma G+\lambda Y^*\right) \right) - Y^* \right\} \right] &= 0, \\ \lambda \mathop{\mathrm{E}} \left[ \frac{ \rho'' \left( \textrm{prox}_{\lambda\rho}\left(\xi + \sigma G+\lambda Y^*\right) \right) }{ 1 + \lambda \rho'' \left( \textrm{prox}_{\lambda\rho}\left(\xi + \sigma G+\lambda Y^*\right) \right) } \right] &= \kappa, \\ \mathop{\mathrm{E}} \left[ \left\{ \rho' \left( \textrm{prox}_{\lambda\rho}\left(\xi + \sigma G+\lambda Y^*\right) \right) - Y^* \right\}^2 \right]&= \frac{\sigma^2\kappa}{\lambda^2} \end{aligned} \]
Scalar characterisation of MDYPL limit
Idea
For MDYPL estimates \((\hat{\theta}^{\textrm{\small DY}}, \hat{\beta}^{\textrm{\small DY}})\), decompose
\[ \hat{\beta}^{\textrm{\small DY}}= B \hat{v} + E \hat w \,, \] with \(E\) spanning \(\operatorname{range}(B)^\perp\)
Limiting parameters \({\color[]{#E16A86}{(}\theta^*, v^*, \sigma^*, \lambda^*)}\) give low-dimensional characterisation of MDYPL through CGMT value comparisons on localising sets
\[ \left\{ \left| \frac{\hat w}{\sqrt p} - {\color[]{#E16A86}{\sigma}^*} \right| > \epsilon \right\}, \quad \left\{ \left| \Gamma_p^{+1/2} \hat v - {\color[]{#E16A86}{u}^*} \right| > \epsilon \right\}, \quad \left\{ \left| \hat{\theta}^{\textrm{\small DY}} - {\color[]{#E16A86}{\theta}^* }\right| > \epsilon \right\} \]
MDYPL asymptotics
Informal result
\(\displaystyle \hat{\beta}^{\textrm{\small DY}}\approx {\color[]{#E16A86}{v}_1^*}\beta_0 + {\color[]{#E16A86}{v}_2^*}\beta_P + {\color[]{#E16A86}{{\sigma}_{\star}}}\varepsilon,\quad \hat{\theta}^{\textrm{\small DY}} \to {\color[]{#E16A86}{\theta}^*},\quad \text{as } n \to \infty\)
Theorem 1 (Convergence of the MDYPL estimator) Assume the conditions of the setup, and for \(\beta^* = B v^*\), the MDYPL estimator satisfies:
(i) \[ \hat{\theta}^{\textrm{\small DY}} \overset{\textrm{a.s.}}{\longrightarrow} \theta^*. \]
(ii) For every pseudo-Lipschitz function \(\psi:\Re^3\to\Re\) of order two1, \[ \frac{1}{p} \sum_{j=1}^{p} \psi\left( \hat{\beta}^{\textrm{\small DY}}_j-v_1^* \beta_{0,j} - v_2^* \beta_{P,j}, {\beta}_{0,j}, {\beta}_{P,j} \right) \overset{\textrm{a.s.}}{\longrightarrow} \mathop{\mathrm{E}}\left[ \psi(\sigma^*G,\bar\beta_0,\bar\beta_P) \right]. \] Here \((\bar\beta_0,\bar\beta_P)\sim \pi_{\{\bar\beta_0,\bar\beta_P\}}\) and \(G\sim\mathop{\mathrm{N}}(0,1)\) is independent of \((\bar\beta_0,\bar\beta_P)\).
1 \(\psi:\Re^d\to\Re\) is pseudo-Lipschitz of order \(a\geq1\) if, for some \(L_\psi>0\), \(\lvert\psi(x)-\psi(y)\rvert\leq L_\psi(1+\lVert x\rVert_2^{a-1}+\lVert y\rVert_2^{a-1}) \lVert x-y\rVert_2\) for all \(x,y\in\Re^d\).
Aggregate asymptotics
| \(\psi(a,b_0,b_P)\) | Statistic | Quantity | Almost-sure limit |
|---|---|---|---|
| \(a\boldsymbol b\) | \(\displaystyle \frac{1}{p}\boldsymbol B^\top(\hat{\beta}^{\textrm{\small DY}}-\boldsymbol\beta^*)\) | Residual alignment | \(\boldsymbol 0_2\) |
| \(a^2\) | \(\displaystyle \frac{1}{p}\left\|\hat{\beta}^{\textrm{\small DY}}-\boldsymbol\beta^*\right\|_2^2\) | Residual variance | \((\sigma^*)^2\) |
| \(\left(a+\boldsymbol v^{*\top}\boldsymbol b\right)^2\) | \(\displaystyle \frac{1}{p}\left\|\hat{\beta}^{\textrm{\small DY}}\right\|_2^2\) | Squared slope norm | \(\displaystyle (\sigma^*)^2+\boldsymbol v^{*\top}\boldsymbol\Gamma\boldsymbol v^*\) |
| \(\left\{a+(\boldsymbol v^*-\boldsymbol e_1)^\top\boldsymbol b\right\}^2\) | \(\displaystyle \frac{1}{p}\left\|\hat{\beta}^{\textrm{\small DY}}-\boldsymbol\beta_0\right\|_2^2\) | Slope MSE | \(\displaystyle (\sigma^*)^2+(\boldsymbol v^*-\boldsymbol e_1)^\top\boldsymbol\Gamma(\boldsymbol v^*-\boldsymbol e_1)\) |
| \(\displaystyle a^2/(v_1^*)^2\) | \(\displaystyle \frac{1}{p}\left\|\frac{\hat{\beta}^{\textrm{\small DY}}-v_2^*\boldsymbol\beta_P}{v_1^*}-\boldsymbol\beta_0\right\|_2^2\) | Oracle-adjusted slope MSE | \(\displaystyle (\sigma^*)^2/(v_1^*)^2\) |
\(\hat{\beta}^{\textrm{\small DY}}\) for \(\alpha = 1 / (1 + \kappa)\)
\(\hat{\beta}^{\textrm{\small DY}}\) for \(\alpha = 1 / (1 + \kappa)\)
\(\hat{\beta}^{\textrm{\small DY}}/ v^*\) for \(\alpha = 1 / (1 + \kappa)\)
Testing
Let \(I=I_p\subseteq\{1,\ldots,p\}\), with \(|I|=k\) fixed, and consider \(H_0: \, \beta_{0,I}=b_I^0\)
Adjusted \(Z\)-statistics
For \(v_1^*\neq0\), let \(\displaystyle \check{\beta} = (\hat{\beta}^{\textrm{\small DY}}-v_2^*\beta_P) /v_1^*\)
\[ \frac{\check{\beta}_I-b_I^0}{(\sigma^* / v_1^*)} \overset{\mathrm d}{\longrightarrow} \mathrm N(0_k, I_k) \]
Penalized likelihood-ratio statistic
Oracle null center: \(m_I^0=v_1^*b_I^0+v_2^*\beta_{P,I}\)
Test statistic: \(\displaystyle \Lambda_I = \max_{\theta,\beta} \, \ell(\theta,\beta; Y^*, X) - \max_{\substack{\theta,\beta:\\ \beta_I=m_I^0}} \, \ell(\theta,\beta;Y^*, X)\)
\[ \frac{2\lambda^*}{(\sigma^*)^2}\Lambda_I \overset{\mathrm d}{\longrightarrow} \chi_k^2 \]
Prediction
For fixed \(\alpha\) and \(j\in\{\mathrm{DY},\mathrm{adj}\}\,,\) define \[ \hat\eta^{\mathrm{DY}}(x) = \theta^{\textrm{\small DY}}+x^\top\hat{\beta}^{\textrm{\small DY}}\,, \qquad \hat\eta^{\mathrm{adj}}(x) = \theta^{\textrm{\small DY}}+x^\top\check\beta \,, \] where the adjusted score requires \(v_1^*\neq0\,.\) For a fresh observation, \[ \bigl\{ \eta_0(x_{\rm new}), \hat\eta^j(x_{\rm new}) \bigr\} \overset{\mathrm d}{\longrightarrow} (S_0,S^j) \] where \((S_0,S^j)\) are bivariate Gaussian with known mean and covariance
Classification
explicit limiting misclassification risk \(\mathcal E_\alpha^j(h)\) for every threshold \(h\)
closed-form optimal threshold and an oracle criterion for tuning \(\alpha\)
Probability prediction
explicit expression for limiting oracle map \(\pi_{S^j}(s)=\Pr(Y=1\mid S^j=s)\)
calibrated and Bregman-optimal among all probability predictions based on \(S^j\)
MDYPL performance
Predicting “7”
MDYPL
| Method | Statistic | p-value |
|---|---|---|
| ML | 0.00 | 1.000 |
| MDYPL | 63.02 | 0.511 |
| Rescaled MDYPL | 163.59 | 1.14e-10 |
Extensions
Covariate structure
Arbitrary covariance: \(x_{i} \sim \mathrm{N}(0, \Sigma)\)
Subgaussian covariate distributions
Estimating unknowns
- Consistent estimation procedure of \(\gamma, \varphi, \theta_0\)
Outlook
Feedforward neural networks
Results
References
Thank you!
Sterzinger P, Kosmidis I (2026). Diaconis-Ylvisaker prior penalized likelihood for \(p / n \to \kappa \in (0, 1)\) logistic regression. Biometrika. 113(2), asag014 DOI: 10.1093/biomet/asag014
Philipp Sterzinger – High-dimensional logistic regression