General theme

Estimation problem

\[ \hat \theta \in \underset{\theta \in \Theta}{\arg \min} \, \ell(\theta; X, y) \]

Penalised estimation

\[ \tilde \theta \in \underset{\theta \in \Theta}{\arg \min} \, \ell(\theta; y, X) + P(\theta) \]

Penalty function

\[ \lim\limits_{n \to \infty} P(\theta_n) = \infty \]

for \(\{\theta_n\}_{n \in \mathbb{N}}\) such that either

1. \(\lim\limits_{n \to \infty} \| \theta_n \| = \infty\), or
2. \(\lim\limits_{n \to \infty} \theta_n \in \partial \Theta\), and

Improved penalised estimation

Maximum softly penalized likelihood

• Soft scaling of penalty function to make distortions of penalty asymptotically negligible

• Preserve potential optimal properties of original estimator

Logistic regression model

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 | x_i \sim \mathop{\mathrm{Bernoulli}}(\pi_i)\,, \quad \log \frac{\pi_i}{1 - \pi_i} = \eta_i = x_i^\top \beta \]

Maximum likelihood estimation

Log-likelihood

\(\displaystyle l(\beta; y, X) = \sum_{i = 1}^n \left\{y_i \eta_i - \log\left(1 + e^{\eta_i}\right) \right\}\)

Maximum likelihood (ML) estimator

\(\hat{\beta} = \arg \max l(\beta; y, X)\)

Artificial example

\[ \begin{array}{l} n = 2000 \\ p = 400 \\ x_{ij} \stackrel{\text{iid}}{\sim} \mathop{\mathrm{N}}(0, 1/n) \\ \\ \beta_0 = (b_1, \ldots, b_{200}, 0, \ldots, 0)^\top \\ b_j \stackrel{\text{iid}}{\sim} \mathop{\mathrm{N}}(6.78, 1) \\ \end{array} \implies \begin{array}{l} \gamma^2 = \mathop{\mathrm{var}}(x_i^\top \beta_0) \approx 4.7 \\ \kappa = 0.2 \\ \end{array} \]

\[ y_i \stackrel{\text{ind}}{\sim} \mathop{\mathrm{Bernoulli}}\left(\pi_i\right), \quad \log\frac{\pi_i}{1 - \pi_i} = x_i^\top \beta_0 \]

Frequentist performance

Recent developments

Candès & Sur (2020)

sharp phase transition about when the ML estimate does not exist, when \(\eta_i = \delta + x_i^\top \beta\), \(x_i \sim \mathop{\mathrm{N}}(0, \Sigma)\), \(p/n \to \kappa \in (0, 1)\), \(\mathop{\mathrm{var}}(x_i^\top\beta_0) \to \gamma^2\)

Phase transition, \(\hat\beta\) (Candès & Sur, 2020)

Phase transition, \(\hat\beta / {\mu}_{\star}\) (Sur & Candès, 2019)

Diaconis-Ylvisaker prior

Prior

\[ \log p(\beta; X) = \frac{1 - {\color[]{#E16A86}{\alpha}}}{{\color[]{#E16A86}{\alpha}}} \sum_{i=1}^n \left\{ \zeta'\left(x_i^\top {\color[]{#E16A86}{\beta}_P}\right) x_i^\top \beta - \zeta\left(x_i^\top \beta\right) \right\} + C \] with \(\zeta(\eta) = \log(1 + e^\eta)\)

Penalized log-likelihood

\[ l^*(\beta; y, X) = \frac{1}{\alpha} \sum_{i = 1}^n \left\{y_i^* x_i^\top\beta - \zeta\left(x_i^\top \beta\right) \right\} \] with \(y_i^* = \alpha y_i + (1 - \alpha) \zeta'\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 l^*(\beta; y, X)\)

MDYPL implementation using ML procedures

MDYPL is ML with pseudo-responses: \(l(\beta; y^*, X) / \alpha = \ell^*(\beta; y, X)\)

Existence and uniqueness

Rigon & Aliverti (2023, Theorem 1): \(\hat{\beta}^{\textrm{\small DY}}\) is unique and exists for all \(\{y, X\}\) configurations

Shrinkage

\(\displaystyle \hat{\beta}^{\textrm{\small DY}}\longrightarrow \left\{ \begin{array}{ll} \hat\beta\,, & \alpha \to 1 \\ \beta_P\,, & \alpha \to 0 \end{array}\right.\)

Rigon & Aliverti (2023) suggest adaptive shrinkage with \(\alpha = 1 / (1 + \kappa)\)

Shrinkage direction

Fixed \(p\) (Cordeiro & McCullagh, 1991)

ML estimator of a logistic regression model is biased away from \(0\)

\(p / n \to \kappa \in (0, 1)\) (Sur & Candès, 2019; Zhao et al., 2022)

Empirical evidence that this continues to hold in high dimensions

\[\LARGE\Downarrow\]

\[\beta_{P} = 0\]

Setting

Signal \(\beta_0\)

The empirical distribution of \(\beta_0\) elements converges to \(\pi_{\bar{\beta}}\) and \(\sum_{j = 1}^p \beta_{0,j}^2 / p \to \mathop{\mathrm{E}}(\bar{\beta}^2) < \infty\)

Covariate distribution, dimension, and signal strength

\(x_{ij} \stackrel{iid}{\sim} \mathop{\mathrm{N}}(0, 1 / n)\)

As \(p / n \to \kappa \in (0, 1)\), \(\mathop{\mathrm{var}}(x_i^\top \beta_0) \to \gamma^2 > 0\)

Idea

Develop a generalized AMP recursion (see Feng et al., 2022 for an overview of AMP), whose iterates have a known asymptotic distribution, with stationary point \(\hat{\beta}^{\textrm{\small DY}}\)

AMP recursion

\[ \begin{aligned} \hat{\beta}^t &= \hat \beta^{t - 1} + \kappa^{-1} \Psi^{t-1}(y^*, S^{t-1}) \\ S^t &= X \hat \beta^t - \Psi^{t-1}(y^*, S^{t-1}) ,, \end{aligned} \] for \[ \Psi^t(y_i^*, S^t_i) = b_t (y^*_i - \zeta'(\textrm{prox}_{b_t \zeta}\left(b_t y_i^* + S_i\right))) \] and where \(\textrm{prox}_{b\zeta}\left(x\right)=\arg\min_{u} \left\{ b\zeta(u) + \frac{1}{2} (x-u)^2 \right\}\) is the proximal operator

Fixed point satisfies \[ X^\top (y^* - X\beta^*) = 0 \]

State evolution of AMP recursion

Solution \(({\mu}_{\star}, {b}_{\star}, {\sigma}_{\star})\) to the system of nonlinear equations in \((\mu, b, \sigma)\)

Solvers for nonlinear systems, after cubature approximation of expectations

MDYPL asymptotics

Behaviour of \(\hat{\beta}^{\textrm{\small DY}}\)

Behaviour of \(\hat{\beta}^{\textrm{\small DY}}/ {\mu}_{\star}\)

\(\hat{\beta}^{\textrm{\small DY}}/ {\mu}_{\star}\) for \(\alpha = 1 / (1 + \kappa)\)

\(\hat{\beta}^{\textrm{\small DY}}/ {\mu}_{\star}\) for \(\alpha = 1 / (1 + \kappa)\)

Abritrary covariate covariance

\(\hat{\beta}^{\textrm{\small DY}}\) is the mDYPL estimator with covariate vectors \(x_{i} \stackrel{iid}{\sim} \mathop{\mathrm{N}}(0, \Sigma)\) and signal \(\beta_0\)

Dimension and signal strength

As \(p / n \to \kappa \in (0, 1)\), \(\mathop{\mathrm{var}}(x_i^\top \beta_0) = \beta_0^\top \Sigma \beta_0 \to \gamma^2 > 0\)

Aggregate asymptotics

Under the assumption \(\underset{n \to \infty}{\limsup } \, \lambda_{\max}(\Sigma) / \lambda_{\min}(\Sigma) < \infty\), Theorem 1 still holds with \[ \frac{1}{p} \sum_{j=1}^p \psi \left( {\color[]{#E16A86}{\sqrt{n} \tau_j}} \left( \hat{\beta}^{\textrm{\small DY}}_{j}-{\mu}_{\star}\beta_{0,j}\right), {\color[]{#E16A86}{\sqrt{n} \tau_j}} \beta_{0,j} \right) \overset{\textrm{a.s.}}{\longrightarrow} \mathop{\mathrm{E}}\left[\psi({\sigma}_{\star}G, \bar{\beta})\right] \] where \(\tau_j^2 = \mathop{\mathrm{var}}(x_{ij} \mid x_{i, -j})\)

Inference

Adjusted \(Z\)-statistics

DY prior penalized likelihood ratio test statistics

Frequentist performance

Current work

Inclusion of intercept

\[ Y_i | x_i \sim \mathop{\mathrm{Bernoulli}}(\pi_i)\,, \quad \log \frac{\pi_i}{1 - \pi_i} = \eta_i ={ \color[]{#E16A86}{\theta}_{0}} + x_i^\top \beta_0 \]

Nonzero prior signal

\[ y_i^* = \alpha y_i + (1-\alpha) \zeta'({\color[]{#E16A86}{\theta}_P} + x_i^\top {\color[]{#E16A86}{\beta_P}}) \]

with

\[ \operatorname{cov}\left(\theta_0 + x_i^\top \beta_0, \theta_P + x_i^\top \beta_P\right) \overset{\text{a.s.}}{\longrightarrow} \begin{bmatrix} \gamma^2 & \varphi \\ \varphi & \delta^2 \end{bmatrix} \]

CGMT (Thrampoulidis et al., 2018)

Let

\[ \begin{aligned} \Phi(G) &= \underset{w \in \mathcal{S}_w}{\min} \, \underset{u \in \mathcal S_u}{\max} \, w^\top G u + \psi(u, w) && (\text{PO})\\ \phi(g,h) &= \underset{w \in \mathcal{S}_w}{\min} \, \underset{u \in \mathcal S_u}{\max} \, \|u\|_2 g^\top w + \|w\|_2 h^\top u + \psi(u, w) && (\text{AO}) \end{aligned} \] where \(G \in \Re^{n \times m}\), \(g \in \Re^n\), \(h \in \Re^m\), \(S_w \subset \Re^n\), \(S_u \subset \Re^m\), and \(\psi : \Re^n \times \Re^m \rightarrow \Re\) is convex concave.

If \(\mathcal{S}_u, \mathcal{S}_w\) are convex, compact sets and \(G, g, h\) have i.i.d. \(\mathrm{N}(0,1)\) entries, then

1. Tail bounds of the PO in terms of the AO
2. Performance metrics of the PO optimiser converge to the AO-predicted limits

CGMT pipeline

Original problem

\[ \hat \theta,\, \hat \beta = \underset{\substack{\theta \in \Re \\ \beta \in \Re^p}}{\min} \, \frac{1}{n} \left\{1_n^\top \zeta(\theta 1_n + X\beta) - (y^*)^\top (\theta 1_n + X\beta) \right\} \]

Primary optimisation (PO)

\[ \underset{\substack{\theta \in \Re \\ \beta \in S_1}}{\min} \, \underset{\lambda \in S_2}{\max} \, \frac{1}{n} \left\{1_n^\top \zeta(\eta) - \eta^\top y^* + \lambda^\top \left(\eta - \theta 1_n - X\beta \right) \right\} = \underset{\beta \in S_1}{\min} \, \underset{\lambda \in S_2}{\max} \, \frac{1}{n} \lambda^\top X \beta + \psi(\lambda) \]

Auxiliary optimisation (AO)

\[ \underset{\beta \in S_1}{\min} \, \underset{\lambda \in S_2}{\max} \, \| \lambda \|_2 g^\top \beta + \|\beta\|_2 h^\top \lambda + \psi(\lambda) \]

CGMT pipeline

Limiting AO

CGMT link

\[ \begin{aligned} \hat \theta &\overset{\text{a.s.}}{\longrightarrow} \theta_{\star}\\ \hat \beta &\approx \mu_{1,\star} \beta_0 + \mu_{2,\star} \beta_P + {\sigma}_{\star}\varepsilon \end{aligned} \] where \((\mu_{1,\star}, \mu_{2,\star}, {\sigma}_{\star}, \lambda_{\star}, \theta_{\star})\) are the solution to Equation 2, and \(\varepsilon\) is some mean zero noise

Opens door to inference, subgaussian isotropic covariates