Multinomial Logistic Regression: Asymptotic Normality on Null Covariates in High-Dimensions
Authors: Kai Tan
Pierre C Bellec
Venue: NeurIPS 2023 poster (2023)
Keywords: High-dimensional statistics
statistical inference
multi-class classification
asymptotic normality
multinomial logistic regression

Abstract: This paper investigates the asymptotic distribution of the maximum-likelihood estimate (MLE) in multinomial logistic models in the high-dimensional regime where dimension and sample size are of the same order. While classical large-sample theory provides asymptotic normality of the MLE under certain conditions, such classical results are expected to fail in high-dimensions as documented for the binary logistic case in the seminal work of Sur and Candès [2019]. We address this issue in classification problems with 3 or more classes, by developing asymptotic normality and asymptotic chi-square results for the multinomial logistic MLE (also known as cross-entropy minimizer) on null covariates. Our theory leads to a new methodology to test the significance of a given feature. Extensive simulation studies on synthetic data corroborate these asymptotic results and confirm the validity of proposed p-values for testing the significance of a given feature.

Multinomial Logistic Regression: Asymptotic Normality on Null Covariates in High-Dimensions
Kai Tan Department of Statistics Rutgers University
Pierre C. Bellec
Department of Statistics
Rutgers University
Piscataway, NJ 08854 kai, tangrutgers..... edu
Piscataway, NJ 08854
pierre. bea.reebbbleed@rugerrr...e
Abstract
This paper investigates the asymptotic distribution of the estimate (MLE) in multinomial logistic models in the high-dimensional regime where dimension and sample size are of the same order. While classical large- sample theory provides asymptotic normality of the MLE under certain conditions, such classical results are expected to fail in high- digh-dimesssos as documented for the binary logistic case in the seminal work of Sur and Candès [2019]. We address this issue classification problems with more classes, by developing asymptotic normality and asymptotic chi-square results for the multinomial logistic MILE (also known as cross-entropy minimizer) on null covariates. Our theory leads to a new methodology to test the significance of given feature. Extensive simulation studies on synthetic data corroborate these asymptotic results and confirm the validity of proposed p-values for testing the significance of given feature.
Introduction
Multinomial logistic modeling has become a cornerstone of classification problems in machine learning, as witnessed by the omnipresence of both the cross-entropy loss (multinomial logistic loss) and the softmax function (gradient of the multinomial logistic loss) in both applied and theoretical machine learning. We refer to Cramer (2002) for an account of the history and early developments of logistic modeling.
Throughout, we consider a classification problem with K possible labels where K a fixed constant. This paper tackles asymptotic distributions of multinomial logistic estimates (or cross- entropy minimizers) in generalized linear models with moderately high-dimensions, where sample size n and dimension p have the same order, for instance n.p - +oo simultaneously while the ratio p/n converges a finite constant. Throughout the paper, let [n] {1, ...,. n} for all n E N, and I [statement] be the 0-1 valued indicator function, equal to I if statement is true and 0 otherwise (yy: = 1} in the next paragraph equals y = holds and otherwise).
The case of binary logistic regression. Let p()) log(1 e²) be the logistic loss and ((( (e-) be its derivative, often referred to as the sigmoid function. In the current moderately- high dimensional regime where n + 400 with p/n -> 6. > () for some constant k. recent works [Candès and Sur 2020 (Sur and Candès] 2019, /Zhao al., 2022] provide detailed theoretical understanding of the behavior of the logistic Maximum Likelihood Estimate (MLE) in binary logistic regression models. Observing independent observations (ær, y4 ye)vet from a logistic model defined as P(yy lx;)) = p' (xo (8) where v; ~ N(O, n -11), and lim ||,11|/n = 22 for a constant y for the limiting squared norm of the unknown regression vector B. These works prove that the behavior the MLE arg minbere E p(x: 1 1)aef is summarized by the solution
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

(4x,0*, 02, de) of the system of three equations
(- + VKO /rrrriii//K//tt
+ Vna 22)) = //2/ xp"ro=xxsppsx
(1.1)
+ Vka 22)))
where (2,, are i..d. ((0, 1) random variables and the proximal operator is defined as (2) = arg min Ap(t) + (t /)] The system 1.1 characterize, among others, the following behavior of the MLE 3: for almost any (7, ki), the system admits a solution if and only if 3 exists with probability approaching one and in this case, |/ and || - Bll 2112/n both have finite limits that may be expressed as simple functions of (a..; o., A.). and for any feature j E [p] such that B3 = 0 (i.e.. j is a null covariate), the j-th coordinate of the MLE satisfies
B 4 NO,
The proofs in Sur and Candes 2019| based on approximate message passing (AMP) techniques; we refer to Berthier et al. [2020], Feng et al.] [2022], Gerbelot and Berthier [2021] and the references therein for recent surveys and general results. More recently, Zhao et al. 2022] extended the result Sur and Candès [019)] from isotropic design to Gaussian covariates with an arbitrary covariance structure: if now 24 (O, Σ) for some positive definite Σ and pim
K; null covariates j E [p] (in the sense that yi is independent of Tıj given (Fa))e())))) of the MLE satisfy
Σ, EB
(n/S2
4 N (0,
(1.2)
where o. is the same solution of 1 and Ω = E-1 Zhao et al. 2022 also obtained asymptotic normality results for non-null covariates, that is, features displays can be used to test the null hypothesis Ho Yt is independent of Tij given (Ii*) (GE))) M/) and develop the corresponding p-values if a. is known; in this binary logistic regression model the
( such that [3 # 0. The previous
ProbeFrontier (Sur and Candès [019] and SLOE Yadlowsky et al. [2021] give means to estimate
the solutions (a., ).) of system (I.I) without the knowledge of 7. Mai et al. [2019] studied the
performance of Ridge regularized binary logistic regression in mixture models, Salchi et al. [2019]
extended [Sur and Candès [2019] to separable penalty functions, [Bellec [2022] derived asymptotic
normality results similar to (1.2) 2) in single-indes models including binary logistic regression without
resorting to the system I.I showing that for a null covariate j E [p] in the unregularized case that
(6/A))
(1.3)
where 0 = + - - (f 3 E aip" (xf 9)af) scalar and so #2 (1 1 -1
(x 8))². In summary, in this high dimensional binary logistic model,
(i) The phase transition from Candès and Sur 2020 splits the (7, k) plane into two connected components: in one component the MLE does not exist with high probability, in the other component the MLE exists and / 12 is bounded with high probability (boundedness is
a consequence of the fact that | /112 or | 3)1 admit finite limits);
(ii) In the component of the (7, k) plane where the MLE exists, for any null covariate j [p]. the asymptotic normality results
1.3 holds.
Multiclass classification. The goal of this paper is to develop theory for the asymptotic normality of the multinomial logistic regression MLE (or cross-entropy minimizer) on null covariates when the number of classes, K + 1, is greater than and n, p are of the same order. In other words, we aim to generalize results such as (1.2) or I.3p for three or more classes. Classification datasets with 3 or more classes are ubiquitous in machine learning (MNIST, CIFAR name few), which calls for such multiclass generalizations. In Gaussian mixtures and logistic models, Thrampoulidis et al. 22000] derived characterizations of the performance of of geast-squares and class-averaging estimators. excluding cross coosspyn minimizers or minimizers of non- linear losses. Loureiro et al. [221]] extended (Sur and Candès 2019) Zhao et al. [2022) Salehi et al. 2019 to multiclass classification problems in Gaussian mixture model, and obtained the fixed-point equations that characterize the performance and empirical distribution of the minimizer of the cross-entropy loss plus convex
regularizer. In the same vein Loureiro et al. 2002]] Cornacchia et al. [002 studied the limiting

---

fixed-pont point equations in a multiclass teacher-studen, learning model where labels are generated by noiseless channel with response arg minke (K x; BR where B&. E IP is unknown for each class k. These two aforementioned works assume multiclass Gaussian mixture model, which is different than the normality assumption for used in the present paper. More importantly, these results cannot be readily used for the purpose testing significant covariates (ef. (1.10) below) since solving the fixed-point equations require the know lnowledge of several unknown parameters, including the limiting spectrum of the mixture covariances and empirical distributions of the mixture means (cf. for instance Corollary 3 in Loureiro et al. [202T]. )) In the following sections, we fill this gap with new methodology to test the significance of covariates. This is made possible by developing new asymptotic normality results for cross-entropy minimizers that generalize (1.3) without relying on the low-dimensional fixed-point equations.
Notation. Throughout, Ip E IRPP is the identity matrix, for matrix A Rmx", AA denotes the transpose of A, AT denotes the Moore-Penrose inverse of A. If A is psd, A1/2 denotes the unique symmetric square root, i.e., the unique positive semi-definite matrix such that (A1/2)2 A. The symbol @ denotes the Kronecker product of matrices. Given two matrices A E Rn*k B E R=*g with the same number or rows, (A, B) E R=*(k+*) is the matrix obtained by stacking the columns of A and B horizontally. If R" is a column vector with dimension equal to the number of rows in A, construct (A, v) E (k++(++) similarly. We use O,, and 1, denote the all-zeros vector and all-ones vector in R", respectively; we do not bold vectors and matrices other than O, and 1 We 11 may omit the subscript giving the dimension if clear from context; e.g., in /K+1 K+T the vector 1 in RK The Kronecker product between two matrices is denoted by @ and veo(M) E Rndd is the vectorization operator applied matrix M Ruxd For an integer K and a E (0, 1), the quantile x (ax) is the unique real number satisfying P(W > X (a)) where W has a chi-square distribution with degrees of freedom. The symbols and denote convergence in distribution and in probability.
Throughout, classical asymptotic regime refers to the scenario where the feature dimension p is fixed and the sample size n goes to infinity. In contrast, the term high-dimensional regime refers to the situation where n and p both tend to infinity with the ratio p/ n converging to a limit smaller than 1.
1.1 Multinomial logistic regression
Consider multinomial logistic regression model with K + classes. We have i.i.d. data samples y:) where x; E Rº is the feature vector and y; (ys1, Y.(k
R++1 is the response. Each response y; is the one-hot encoding of a single label, ie., y, E {0,1]K*1 with kk+1 Yik = 1 such that Yı = 1 if and only if the label for i-th observation is k. A commonly used generative model for y; the multinomial regression model, namely
cxp(x B ex)
P(ya = 1/2,)
E
K
(1.4)
Σ Σ exp(x: (
where B" E (k+( is an unknown logistic model parameter and e₁ E K+1 ew E RR+ are the k-th and w'tth canonical basis vectors. The MLE for B* the model (1.4 any solution that minimizes the cross-entropy loss,
B E arg minBeRxx(X11) Ex= L,(BT..) x.))
(1.5)
where L; R K+1 - R is defined L;(u) Lkm1 Uu log Et'1 exp(ugr). If the solution set in (15) is non-empty, we define for each observation i E [n] the vector of predicted probabilities with
p. (Pat, P((K+1))*
exp [z Bex) Σ
Pu
Plys = 1)
for each & E {L, K 1}.
(1.6)
E exp(x! Bexr
K+I
Our results will utilize the gradient and Hessian of L, evaluated at B x,, denoted by
g = VB, 24) yy + pi,
H = (BL 22) diag(p.) P.P.. (1.7)
3

---

The quantities (B, P.. H;) can be readily computed from the data {2),,))) To be spe- cific, the MLE B in (.5)
be obtained by invoking multinomial regression solver (e.g.. LogisticRegresseee
sklearn. linear wwieeannnoeellle model ties Ps: gr, H; can be further computed from eqs. (1.6) and (.7) by few matrix multiplications and application of the softmax function.
from Pedregosa et al. (201T]. and the quanti-
Log-odds model and reference class. The matrix B" in 1.4) is not identifiable since the conditional distribution of y; [z in the model (1.4) remains unchanged if we replace columns of B" by B -b1 T for any & E IR.. In order to obtain identifiable model, classical and natural remedy is to model the log-odds, here with the class K + + as the reference class:
K +11
P(ya 1/2;) = x! eRe log P(y:(++1)= 1/x4)
Vk E [K]
(1.8)
where 6& is the k-th canonical basis vector of RKK and A E RPXK is the unknown parameter. The matrix A E kxxX in log-odds model (L.8) is related to B E Rex(K+1) in the model (.44) by A = (IK. -Ix)². This log-odds model has two benefits: First it is identifiable since the unknown matrix A is uniquely defined. Second, the matrix A" lends itself well to interpretation as its k-th column represents the contrast coefficient between class is and the reference class + 1.
The MLE Å of A in (1.8) is Å arg min ACRPXX is non-empty, A is related to any solution in I.5 by A = B(Ik. -1x)". Equivalently,
L(((A,0p)² 2.). If the solution set in 1.5
Áςa Byk By(K+1)
(1.9)
for each E [p] and & [K]. If there are three classes (i.e. K + 1= 3), this parametrization allows us to draw scatter plots of realizations of ynez Á = (VπiÁ₁ı, vnAA22) as in Figure !
1.2 Hypothesis testing for the j-th feature and classical asymptotic normality for MLE
Hypothesis testing for the j-th feature. Our goal is to develop a methodology to test the sig- nificance of the j-th feature. Specifically, for desired confidence level (1 - a)) E (0, 1) (say. 1 a 0.95) and a given feature j E [p) of interest, our goal is to test
Hn : y: is conditionally independent of mij given (MG)))
(1.10)
Namely, we want to test whether the j-th variable is independent from the response given all other explanatory variables (Iay,.) JJ E [p) \ (i}). Assuming normally distributed x; and multinomial
model as in (1.4) (.8). is equivalent to test
Ho ee. A" 0%
versus
11 :r:eAAA* + 40
(L.11)
where Ej E IR the j-th canonical basis vector.
If the MLE B in I.5) exists in the sense that the solution set in (.5) is nonempty, the conjecture that rejecting Ho when ey B is far from 0k+1 1 is reasonable starting point. The important question, then, to determine quantitative statement for the informal "far from 0k+1", similarly to (.2) in binary logistic regression.
Classical theory with fixed. If is fixed and - in model (.8). likelihood theory [Van der Vaart 1998 Chapter 5] provides the asymptotic distribution of the MLE A, which can be further used to test (1.11). Briefly, if x has the same distribution as any 2,, the MLE Á in the multinomial logistic model asymptotically normal with
classical maximum
(v(cvc(AA veo(A*)) "x,))=~(
where I = E[((x"") @ (diag(m²) n ")
is the Fisher information matrix evaluated at the true parameter A", vec(.) ) is the usual vectorization operator, and x² RK has random entries Tेz each & E [K]. In particular, under Ho e A O%.
exp(x" Atex)/ + E#=1 exp(x" A Atew)) for
Vn'"e, d (O,
(1.12)

---

where S (e) Ix) 2 Ikk - (cov(x))) e,(E --="==")>>~~'"
When (..2
holds, by the delta method we also have vns,
1/2 ATe, d NO, Ix) and
s 1/2 A egll - -
(1.13)
where the limiting distribution is chi-square with K degrees of freedom. This further suggests the size test that rejects Ho when TX Y) > (a), where (x, Y) = n||s; 1/2 "/2 A'ee,ll is the test statistic, If (1.13) holds, this test is guaranteed to have type error converging a. The p-value
of this test is given by
(x) Fxi (t)dt,
(1.14)
where fx (.) is the density of the chi-square distribution with K degrees of freedom.
As discussed in the introduction, Sur and Candès [2019] showed that in binary logistic regression, classical normality results for the MLE such as (1.12 fail in the high-dimensional regime because
the variance (1.12 underestimates the variability of the MLL even for null covariates; see also
the discussion surrounding (1.2). Our goal is to develop, for classification problems with K + 1 2 3
classes, theory that correctly characterize the asymptotic distribution of A ej for null covariate j E [p] in the high-dimensional regime.
We present first some motivating simulations that demonstrate the failure of classical normal ap- proximation [I.2] in finite samples. These simulations are conducted for various configurations of (n.p) with K + = 3 classes. We fix the true parameter A and obtain 1000 realizations of (Az1, A32) by independently resampling the data { / 1000 times. If the result (1.12) holds, then PP(vnATe eg E CC.) + +1 - a, where C₁2- = fu E R ||5- S X2 (a)). Figure 1 displays scatter plots of Vπ(A,1, A52) along with the boundary of 95% confidence set C.. with 0.05. We observe that, across the three different configurations of (n. p). the 95% confidence sets from our theory (Theorem 2.2 presented in next section) cover around 95% of the realizations, while the set Că from classical theory only covers approximately 30% of the points, which is significantly lower than the desired coverage rate of 95%. Intuitively and by analogy with results in binary classification [Sur and Candès| 2019]. this is because the classical theory (1.12) underestimates the variation of the MLE in the high-dimensional regime. Motivated by this failure of classical MLE theory and the results in binary classification [Sur and Candès 2019] among others], the goal of this paper is to develop a theory for multinomial logistic regression that achieves the following objectives:
madem: conmngee24. 3%
modem: COARsgeEO..5
madem: connupec5..35 2%
(a) (n. p) (2000, 600)
(b) (n,p) (3500, 1000)
(c) (n. p) (5000, 1500)
Figure 1: Scatter plot of pairs (VT.A,1, vzzAz2 with K = 2 over 1000 repetitiios. The blue ellipsoid is the boundary of the 95% confidence set for VnATe, ej under H, from the classical MLE theory
(1122 (.13) based on the Fisher information, the dashed red ellipsoids are the boundaries of the 95% confidence set for vnA'e e, under Ho from this paper (cf. (2.3 below). Each of the 1000 repetition gives slightly different dashed ellipsoid. The solid red ellipsoid is the average of these 1000 dashed
ellipsoids. Each row of X is i.i.d. sampled from (( Σ) with Σ (0.5 ) The first [p/4] rows of A are i.i.d. sampled from N(0, Jx) while other rows are set to zeros. We further normalize A such that T = Ix. The last coordinate j = p is used as the null coordinate.
5

---

Establish asymptotic normality of the multinomial MLE A"ez for null covariates as p. +00 simultaneously with finite limit for n/pp
Develop valid methodology for hypothesis testing of (1.10 in this regime, i.c., testing for the presence of an effect of feature j E [p] on the multiclass response.
The contribution of this paper two-fold: (i) For null covariate j E [p]. we establish asymptotic normality results for ATe, that are valid in the high-dimensional regime where n and have the same order; (ii) propose user-friendly test for assessing the significance of feature multiclass classification problems.
2
Main result: asymptotic normality of BTe, and ÅTe, on null covariates
In this section, we present the main theoretical results of our work and discuss their significance. We work under the following assumptions.
Assumption 2.1. For constants 8 > 1, assume that n, - oo with p/n S 8-11 and that the design matrix X E R=*P has n i..d. rows (az)ie(n) ~ (N,, 2) for some invertible E RPxr, The observations (r; Yi)ie(m] are i.i.d. and yach is of the form y; = fU,, x B*) for some deterministic function f, deterministic matrix B* [Px(K+1) such that B=1 Op. and latent random variable U, independent of 24,
Assumption 2.2 (One-hot encoding). The response matrix Y is in [x Its i-th row is one-hot encoded vector, that is, valued in {0, 1]K+1 with K+1 Yıu 1 for each i E [n]]
The model y; 2( B*) for some deterministic ff and B" and latent random variable Ug in Assumption 2.1 is more general than a specific generative model such as the multinomial logistic conditional probabilities in (1.4). broad choices for are allowed. In words, the model y; <U B²) with B*1x+1 = Op means that y; only depends on x; through a K dimensional projection of v, (the projection on the row-space of B*). The assumption p/n S 8-1 more general than assuming fixed limit for the ratio p/n; this allows us to cover low-dimensional settings satisfying p/n n + as well.
The following assumption requires the labels to be "balanced": we observe each class at least yn times for some constant 7 >> (ye)eejn) are i.i.d. as in Assumption 2.1 with distribution independent of n, p, by the law of large numbers this assumption is equivalent to minke[K+1] PP( 1) 0.
Assumption 2.3. There exits a constant 7 E (0, R-1|. such that for each k E [K + 1], with probability approaching one at least yn observations i E [m) are such that yik 1. In other words, P Tyys = 1) 2 y + 1 for each & E [K' 1].
As discussed in item list (i) page 2, in binary logistic regression, Candès and Sur 2020), Sur and Candès [2019] show that the plane [ | [I) is split by smooth curve into two connected open components: in one component the MLE does not exist with high probability, while in the other
component, with high probability the MLE exists and is bounded in the sense that | < <'
or equivalently 1 ||x < T for constants T' independent of n, p. The next assumption requires the typical situation of the latter component, in the current multiclass setting: B in (I.5) exists in the sense that the minimization problem has solutions, and at least one solution is bounded.
Assumption 2.4. Assume P(B exists and |X - 11" I. < nT) + as n, +oo for some large enough constant т.
Note that the validity of Assumption 2.4 can be assessed using the data at hand; if multino- mial regression solver (e.g. sklearn. linear saarrrmmoell._ www.mmeee....iiiiooooo LogisticRegressi
converge and [ K- - || is no larger than a predetermined large constant then we know Assump- tion 2.4 holds. Otherwise the algorithm does not converge or produces an unbounded estimate: we know Assumption 2.4 fails to hold and we need collect more data.
11"
Herre we refer to standard convergence assessment methods for convex solvers, looking at the gradient/Hessian values at the current iterate, looking at the duality gap if available.

---

Our first main result, Theorem
2.1 provides the asymptotic distribution of Bye, where j E [p] is
null covariate, where B is any minimizer B of 1.5) Throughout, denote by 52 the precision matrix defined as Ω = --
Theorem 2.1. Let Assumptions 2. to 2.4 be fulfilled. Then for any [p) such that Hg 1.10)
holds, and any minimizer B of 1.5 we have
T
Big
Vi
eT NN(0, IK+1 K+T
11
(2.1)
1=1
1=1
scolar
((x+))x(++1 RR
con. +(x+1)x((++)
where g: = -v + Ps as in and V, = H - (H, @ DE H @ (zta) ( @ x2)
The proof of Theorem 2.1 is given in Supplementary Section Theorem 2.1 establishes that under Hu. B²e, converges to singular multivariate Gaussian distribution in R R In (2.1) the two matrices/ ->= py) p.) and E V symmetric with kernel being the linear span 11a
of 1 K+1, and similarly, if a solution exists, we may replace B by B( K+1 which is also solution in (1.5). In this case, all matrix-matrix and matrix-vector multiplications, matrix square root and pseudo-inverse, in (2.T) happen with row-space and column space contained in the orthogonal component of 1 K++ so that the limiting Gaussian distribution in K+ is also supported on this K ---iimnnional subspace.
Since the distribution of the left-hand side of (2.1 is asymptotically pivotal for all null covariates j E [p]. Theorem 2.1 opens the door of statistical inference for multinomial logistic regression in high-dimensional settings. By construction, the multinomial logistic estimate Å e R in (I.9 ensures (Á, Op) is minimizer of (1.5). Therefore, we can deduce the following theorem from Theorem 2.1
Theorem 2.2. Define the matrix R (IK, 0x)² E (K))) using block matrix notation. Let
Assumptions 2.1 10 2.4 be fulfilled. For Å in (1.9) and any j E [p] such that Ho in (1.10) holds,
1 IK VK - R
1/2
( EVV.n
V.R AT 4 NN(OK, Ix) (2.2)
B
matrix 12 x(K+1)
soalar
matric BiK ())CK))))
R(K+ixx
R
where 8i is defined in 1.7 and V, is defined in Theorem
Furthermore, for the same j E [v],
(X.) de n
I( E 1/2 (H RAT e,|
satisfies T2 d xxxx- (2.3)
Note that Equations 2.2 and 2..) is stated using Gz5 e When E is unknown, the quantity Qij above results can be repluced by its consistent estimate Pizj defined in 2.5). and the convergence in distribution results still hold.
Theorem 2.2 is proved Supplementary Section
the best of our knowledge, Theorem 2.2 is the first result that characterizes the distribution of null MLE coordinate Å" e; in high-dimensional multinomial logistic regression with 3 or more classes. It is worth mentioning that the quantities (g:, V., A) used in Theorem 2.2 can be readily computed from the data (X, Y). Therefore, Theo- rem 2.2 lets us test the significance of specific feature: for testing Hg, this theorem suggests the test
statistic (X, Y) (2.3) and the rejection region E def {(X, (2(X,Y) 2 (.)). Under the null hypothesis Ho in (1.10). Theorem 2.2 guarantees ((X,Y) E & - a. In other words, the test that rejects Ho X)Y) E & has type error converging to a. The p-value of this test is
p-value _ fx2 (t)dt,
(2.4)
where fx (+) is the density of the chi-square distribution with K degrees of freedom.
Unknown 22 = TT-lee If E unknown, we describe a consistent estimate of the quantity 2 appcaring in (2.1) (2.2 and 2.3) Under the Gaussian Assumption 2.1 the quantity Uzz is the
7

---

12
15.0
Sample quamtles
10
Sample Quanties
Sample Quantler
10.0
8
7.5
5.0
classical modem
classical modem
2.5
classical
modem
2.5 5.0 7.5 10.0
12.5
Theoretical Quantiles
Theoretical Quantiles
Theoretical Quantiles
(a) (n, K) (4000, 2)
(b) (n, K) = (5000,3) 3)
(c) (n, K) (6000, 4)
Figure 2: Q-Q plots of the test statistic in the left-hand side of (113) (in orange) and in the left-hand
side of (2.3) (in blue) for different (n, K) and p 1000.
reciprocal of the conditional variance (ar(zzy/2i..)) ) which is also the noise variance in the linear model of regressing Xe; onto X-, (the submatrix of X excluding the j-th column). According to standard results in linear models, we have R - xx, -x=,
x²-p+11
Since x (n
p + 1) + almost surely by the strong law of large numbers,
G36 (n - + 1)/[l(", - x-)x,xx))-xxx
(2.5)
is a consistent estimator of 2zzz Therefore, the previous asymptotic results in Theorems 2.1 and 2.2 still hold by Slutsky's 's theorem if we replace Ωzj by the estimate 235 in 2.5
3 Numerical experiments
This section presents simulations and real data analysis to examine finite sample properties of the above results and methods. The source code for generating all of the experimental results in this paper can be found in the supplementary material.
Simulation settings. We set = 1000 and consider different combinations of (n, K). The covariance matrix Σ is specified to be the correlation matrix of AR(1) model with parameter p 0.5, that is, == (0.5 pxxp We generate the regression coefficients A E Ro*k once and for all as follows: sample Ao E R with first [p/4] rows being i.i.d. N (0, Ik), and set the remaining rows to O. We then scale the coefficients by defining A" Ao AJ AA))=1/2 so that A = Ik. With this construction, the p-th variable always a null covariate and we use this null coordinate j = p to demonstrate the effectiveness of our theoretical results presented in Theorem 2.2 and the suggested test for testing Ho as described in [.10]. Using the above settings, we generate the design matrix X E _x from N (0, 2), and then simulate the labels from multinomial logistic model as given in (1.8 using the coefficients A E IRPX For each simulation setting, we perform 5,000 repetitions.
A Assessment of X2 approximations. To assess the x² approximation (2.3) from this paper and that of the classical theory (1.13). we compute the two x test statistics for each sample (2,, 31) Figure shows the empirical quantiles of the two statistics versus the X distribution quantiles. The results demonstrate that the quantiles (in blue) from our high-dimensional theory closely match the 45-dege line (in red), whereas the quantiles (in orange) from the classical theory significantly deviate from the 45-degree line. These findings highlight the accuracy of our proposed x2 approximation (2.3 over the classical result (3.13) when p is not sufficiently small compared to
Uniformity of null p-values. Recall that the p-value from the classical test (.13) is given by (1.4)) while the p-value from this paper taking into account high-dimensionality is given by (.4). Figure 3 displays the histograms of these two sets of p-values out of 5000 repetitions. The results in Figure 3 show that the p-values obtained from the classical test deviate significantly from the uniform distribution, with a severe inflation in the lower tail. This indicates that the classical test tends to produce large type errors due to the excess of p-values close to 0. In contrast, the p-values proposed in this paper exhibit uniform distribution, further confirming the effectiveness and applicability of the theory in Theorem 2.2 for controlling type error when testing for null covariates with (I.10).
8

---

classical
3000
classical
classical
moder em
modeer
mod em
2500
frequency
2000
frequency
2000
frequency
1500
1500
1500
1000
1000
1000
500
500
500
0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
p-values
p-values
p-values
(a) (n, K) (4000, 2)
(b) (n2, K) (5000,
(c) (n, K) (6000, 4)
Figure 3: Histogram for p-values of the classical test (.14) (in orange) and of the proposed test (.4
(in blue) under Ho in simulated data with different (n, K) and p 1000.
Unknown Bzj. In the situation where the covariance matrix E is unknown, we can estimate the diagonal element f2yj = =e,-- by Uz defined in (2.5). To evaluate the accuracy of the normal and chi-square approximations and the associated test with $2,J replaced by Ωzz, we conduct simulations similar to those in Figures and but we replace 22 with its estimate Zzz. The results are presented in Figure ST] The plots are visually Ily indistinguissable from the plots using Пz5- These confirm that the chi-square approximation and the associated test using U2 are accurate.
Mon-Gaussian covariates and unknown 2zj- Although our theory assumes Gaussian covariates, we expect that the same results hold for other distributions with sufficiently light tails. To illustrate covariates:
(i) The design matrix X has i.i.d. Rademacher entries, i.e., P( = = = 2 (ii) Each Zij takes on values 0, 1 and 2 with respectively probabilities a². 2u, - ay), and (1 - -aj)², where aj varies in [0.25, 0.75]. Each columns of X are then centered and normalized to have mean and unit variance. This generation of non-Gaussian covariates is adopted from single-nucleotide polymmorpissms (SNPs) example in Sur and Candès [019]. For these two types of non-Gaussian covariates, we further rescale the feature vectors to ensure that x; has the same covariance as in the Gaussian case the beginning of Section 3 that is Σ =
this point, we consider the following two types of non-Gaussee
We present the Q-Q plots in Figure S2 using the same settings as Figure ST with the only difference being that the covariates in Figure S2 are non-Gaussian distributed. The Q-Q plots of T (X, Y') in (233) plotted in Figure SS still closely match the diagonal
line. These empirical successes suggest that the normal and X approximations 2. 2.3) apply to
a wider range of covariate distributions beyond normally distributed data.
Real data example. We conduct a real data analysis by applying the proposed test to heart disease data from the UCI Machine Learning Repository (link: http://arrchiveicsuu wwwp.e//ww/mmcineeeeeaaaea data). After standard data-cleaning processes, the dataset has 297 instances with 13 features, including age, sex, and other attributes. The response variable was transformed into 3 classes (0, 1, and 2) after converting the labels 3 and 4 to 2. To demonstrate the validity of the proposed significance test, we generate noise variable from standard normal distribution, resulting in a dataset with 297 instances and 14 variables. We test the significance of this noise variable using both the proposed chi-square test and the classical test, repeating the experiment 10,000 times. The type I error of our proposed test is 0.0508, aligning well with the desired type I error of 0.05. In contrast, the classical test exhibits type I error of 0.0734, significantly exceeding the desired rate of 0.05. These results confirm the validity of our proposed test on real data corrupted with fake covariate, while the classical test concludes that the fake covariate is significant in more than 7% of experiments, leading to false discoveries exceeding the desired 0.05 type error.
4 Discussion and future work
Multinomial logistic regression estimates and their p-values are ubiquitous throughout the sciences for analyzing the significance of explanatory variables on multiclass responses. Following the seminal work of (Sur and Candès (2019) in binary logistic regression, this paper develops the first valid tests and p- p-values for multinomial logistic estimates when p and are of the same order. For 3 more

---

classes, this methodology and the corresponding asymptotic normality results in Theorems 2.1 and 2.2 are novel and provide new understanding of multinomial logistic estimates (also known as cross- entropy minimizers) in high-dimensions. We expect similar asymptotic normality and chi-square results to within reach for loss functions different than the cross-entropy or different model for the response yı; for instance Section ST provides an extension to the q-repeated measurements model, where q responses are observed for each feature vector Z;.
Let us point a few follow-up research directions that we leave open for future work. A first open problem regards extensions of our methodology confidence sets for ef B* when Hu in (I.10) is violated for the j-th covariate. This would require more stringent assumptions on the generative model than Assumption 2.1 B* there is not identifiable (e.g., modification of both B and ( in Assumption 2.1 is possible without changing y,). A second open problem is to relate this paper's theory to the fixed-point equations and limiting Gaussian model obtained multiclass models, e.g.. Loureiro et al. [2021]. While it may be straightforward to obtain the limit of 1 -i= -igg and of the empiriaal distribution of the rows of Å in this context (e.g., using Corollary in Loureiro et al.) (2021) D. the relationship between the fixed-point equations and the matrix
" >> V, appearing in (2. IS unclear and not explained by typical results from this literature. A third open problem is to characterize the exact phase transition below which the multinomial logistic MLE exists and is
bounded with highh-robabilltyy (Assumption 2.4); while this is settled for two classes [Candès and
Sur 2020] and preliminary results available for 3 or more classes [Loureiro al. | [202]] Kini and Thrampoulidis) 2021]. a complete understanding of this phase transition is currently lacking. A last interesting open problem is to prove that theory extend to non- nonnGaassian data, observed in simulations. This challenging problem is often referred to as "universality" and has received intense
attention recently (Gerace al. 2022 Han and Shen 2022 Montanari and Saeed 2022 Pesce et al., 2023 Dandi et al. 2023]. showing that in several settings of interest (although none exactly matching the one considered in the present paper), the asymptotic behavior of the minimizers is unchanged if the distribution of the covariates is modified from normal to another distribution with the same covariance.
Acknowledgments and Disclosure of Funding
P.C. Bellec's research was partially supported by the NSF Grant DMS-1945428. The authors acknowledge the Office of Advanced Research Computing (OARC) at Rutgers, The State University of New Jersey for providing access to the Amarel cluster and associated research computing resources that have contributed to the results reported here. URL: htttps://oarc.ruterrss
References
Pierre C Bellec. Observable adjustments in single-index models for regularized m-estimators. arXiv preprint arXiv:2204.06949 2022.
Pierre C. Bellec and Cun-Hui Zhang, Second-order Stein: SURE for SURE and other applications in nigh-dimensional inference. Ann. Statist,, 1964 2021. ISSN 0090-5364.
Raphael Berthier, Andrea Montanari, and Phan-Minh Nguyen. State evolution for approximate
message passing with non-separable functions. Information and Inference: Journal of the IMA, 9(1::33-79, 2020.
Emmanuel J. Candès and Pragya Sur. The phase transition for the existence of the maximum likelihood estimate in high-dimensional logistic regression. The Annals of Statistics, 48(1):27 42, 2020.
Elisabetta Comaccchi,, Francesca Mignacco, Rodrigo Veiga, Cédric Gerbelot, Bruno Loureiro, and Lenka Zdeborová. Learning curves for the multi-class teacher-student perceptron. arXiv preprint arXiv:2203. /20424 2022.
Jan Salomon Cramer. The origins of logistic regression. Tinbergen Institute Working Paper, 2002.
Yatin Dandi, Ludovic Stephan, Florent Krzakala, Bruno Loureiro, and Lenka Zdeborová. Universality laws for gaussian mixtures in generalized lincar models. arXiv preprint arXiv:2302.0893. 2023.
10

---

Kenneth R Davidson and Stanislaw J Szarek. Local operator theory, random matrices and banach spaces. Handbook of the geometry of Banach spaces, 1(317-36):131, 2001.
Oliver Y Feng, Ramji Venkataramanan, Cynthia Rush, Richard Samworth, al. A unifying tutorial on approximate message passing. Foundations and Trendso in Machine Learning, 15(4):335-55 2022.
Federica Gerace, Florent Krzakala, Bruno Loureiro, Ludovic Stephan, and Lenka Zdeborová. Gaus- sian universality of linear classifiers with random labels in high-dimension. arXiv preprint arXiv:2205. 13303, 2022.
Cédric Gerbelot and Raphaël Berthier. Graph-based approximate message passing iterations. arXiv preprint arXiv:2109 109. /905, 2021.
Qiyang Han and Yandi Shen. Universality of regularized regression estimators in high dimensions. arXiv preprint arXiv:2206.079000 07936, 2022.
Ganesh Ramachandra Kini and Christos Thrampoulidis. Phase transitions for one-vs-one and one-vs- all linear separability in multiclass gaussian mixtures. ICASSP 2021-2021 IEEE International Conference Acoustics, Speech and Signal Processing (ICASSP), pages 4020-4024. IEEF, 2021.
Bruno Loureiro, Gabriele Sicuro, Cédric Gerbelot, Alessandro Pacco, Florent Krzakala, and Lenka Zdeborová. Learning gaussian mixtures with generalized linear models: Precise asymptotics in high-dimensions. Advances Neural Information Processing Systems, 34:10144-10157, 2021.
Xiaoyi Mai, Zhenyu Liao, and Romain Couillet. A large scale analysis of logistic regression: Asymptotic performance and new insights. In ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3357-3361. IEEF, 2019.
Andrea Montanari and Basil N Saeed. Universality of empirical risk minimization. In Conference on Learning Theory, pages 4310-4312. PMLR, 2022.
Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Pretten- hofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning Python. Journal of Machine Learning Research,
12:2825-288 2830, 2011.
Luca Pesce, Florent Krzakala, Bruno Loureiro, and Ludovic Stephan. Are gaussian data all you need? extents and limits of universality in high-dimensional generalized linear estimation. arXiv preprint wwXiv22302008923 2023.
Fariborz Salehi, Ehsan Abbasi, and Babak Hassibi. The impact of regularization on high-dimensional logistic regression. Advances in Neural Information Processing Systems, 32, 2019.
Pragya Sur and Emmanuel J. Candès. A modern maximum-likelihoodd theory for high- dinhnsininn logistic regression. Proc. Natl. Acad. Sci. USA,
2019. ISSN 0027-8424.
Kai Tan, Gabriel Romon, and Pierre C Bellec. Noise covariance estimation in multi-task high- dimensional linear models. arXiv preprint wrXiv:2206.07255: 2022.
Christos Thrampoulidis, Samet Oymak, and Mahdi Soltanolkotabi. Theoretical insights into mul-
ticlass classification: A high-dimensional asymptotic view. Advances in Neural Information Processing Systems, 33:8907-8920, 2020.
Aad W Van der Vaart. Asymptotic statistics, volume Cambridge university press, 1998.
J Leo van Hemmen and Tsuneya Ando. An inequality for trace ideals. Communications in Mathe- matical Physics, 76: 143-1,84, 1980.
Steve Yadlowsky, Taedong Yun, Cory Y McLeen, and Alexander D' Amouu. Sloe: A faster method for statistical inference in high-dimensional logistic regression. Advances in Neural Information Processing Systems, 34:29517-29527 2021.
Qian Zhao, Pragya Sur, and Emmanuel Candes. The asymptotic distribution of the mle in high- dimensional logistic models: Arbitrary covariance. Bernoulli, (8()):1835-1811 2022.
Zhu and Trevor Hastie. Classification of gene microarrays by penalized logistic regression. Biostatistics, 5(3):427-443 2004.

---

Supplementary Material of "Multinomial Logistic Regression:
Asymptotic Normality on Null Covariates in High-Dimensions"
Let us define some standard notation that will be used in the rest of this supplement. For vector V E IR", let || maX |og] denote the infinity norm of vector v. is symmetric, we define Amin (A) and || as the minimal and maximal cigenvalues of A, respectively. For two symmetric matrices A, B of the same size, write A 3 B if and only if B A is positive semi-definite.
Diagram: Organization of the proofs
The following diagram summarizes the different theorems and lemmas, and the relationships between them.
Theorem S3.1
Theorem 2.2
Asymptotic normality for 8 c, on null covariates, where B R multinomial logistic MLE using the parameter space from Section [S3.]
Theorem 2.1
Asymptotic normality c null covariates, where E xa* is the multinomial logistic MLE with class K + 1 fixed as the reference class (see
Asymptotic normality for B e, null covariates, where B p=x((+1) is multinomial logistic MI.E (1.5).
The proof conditions in The- orem [S5.]] on the loss function are satisfied by the cross-entropy.
that
K--dimensional orthogo- nal parametrization de- fined by the matrix the
Control of g: and H, for the cross-entropy loss
Section S3.I
defines the matrix B x
Theorem S5.1
discusses a conve- nient parametrization of the model isometric to the subspace orthogonal 1x+1-
Lemmas S6.1 and S6.2 give deterministic ar- guments control the gradients and Hessians of the cross-entropy loss. Lemma [S6.4] con- trols the Hessian of the cross-entropy loss the minimizer, in a specific high-probabilityy event. Lemma [S6.3] defines this high-probability event.
Asymptotic normality null covariates for general loss functions, * Ip. Deduced from Theo- rem S5.2 by rotational invariance.
Lemma $5.5
Lemma $5.5 computes the derivatives of the min- imizer with respect to X, used the proof of Theo- rem S5.2
Theorem
Asymptotic normality on null covariates for general loss functions, Σ I
Lemma $5.3
Normal and x approxi- mations for random vari- ables defined as differen- tiable function of standard normal vectors.
12

---

S1 Extension: q repeated measurements
Let integer 2 constant independent of n. p. Our results readily extend if labels are observed for each observed feature vector 2;, and the corresponding q one-hot encoded vectors are averaged into Yi E {0. q--g" ..., 1}*+ Concretely, for each observation E [n], i.id. labels (Y,") (X;"))<([]) observed with each y, {0,1)K+1 one-hot encoded and Yık q E ik for instance a repeated multinomial regression model with P(Y [ equal to right-hand side of (T4).
) ig] are i.i.d., Assumption 2.3 is satisfied by the law of large numbers if [K+1) P(Y = 1) > 0 since is constant. For this q repeated measurements model, the negative log-likelihood function of parameter B E [Rx((+1)) is
In this case where (>;"")
K+1
K+1 exp(x; Bew)
Ym ik
U.
Bey
log
k=1
k
K+1
Yik
2
T
Bes + log
Σ exp(2 ( Be
k'=1
K+1
4
Yikk Be, + log
exp(2 Bex
L:(BTT:))
q
where the first equality uses yuk ạ Σn= Yo, the second equality uses Σk=1 Yıh = under the following Assumption SI.I and the last equality uses the definition of L, after (.5).
Assumption S1.1. For all [ni], the response y, is in {0, 2/. 2/g, 1, K] with K+1 kk= Yık 1.
In such repeated measurements model, we replace Assumption which the following Theorem ST.T
2.2 with Assumption SI.1 under
Theorem S1.1 Let q 2 be constant. Let Assumptions S1.7
2.1
2.3 and 2.4 be fulfilled. For any
j E (p) such that Ho in (1.10) holds, we have the convergence in distribution (2.1). (..2 and (2.3)
Proof of Theorem $1.1
Under the assumptions in Theorem S1.1
the MLE B for this q repeated
measurements model is the minimizer of the optimization problem
arg min ò
E
LL.(BT x:)
BER+x(X+1)
1=1
as in 1.5 Similar to the non-repeated model, the MLE for the identifiable log-odds model can be expressed as
Å = arg min
((A(,,,))'..))
AEEPXR
The only difference between this q repeated measurements model and the non-repeated model considered in the main text is that the response yık for this q repeated measurements model is now valued in (0, 1/g. 2/q,..,, 1}. Because the proofs of Theorems 2.1 a.2 do not require the value of
Yıu be {0. 1 valued.. Theorem S1.1 can be proved by the same arguments used in the proof of
Theorems 2.1 and 2.2
S2 Implementation details and additional figures
The pivotal quantities in our main results Theorems 2.1 and 2.2 involve only observable quantities that can be computed from the data (I;, ys setn)- this section we provide an efficient way of computing the matrix V, appearing in Theorems 2.1 and 2.2
13

---

Fast computation of V,, Recall the definition of V in Theorem 2.1
Σ He @ (mraf) ( (H, @)).
V, H. (H; 8 x55
The majority of computational cost in calculating V; lies in the step of calculating its second term
(H, @ a
Σ Hi @ (maa)) (H: 8 xx
Here we provide an efficient way to compute this term using the Woodbury matrix identity. Since H; k+1 0K+ we have ker(H, @ (xxaf)) the span of [1x+1 @ eg : [y]], where 1x=1 the all-ones vector in RR+1 Therefore, the second term in V; can be rewritten as
(H: 2 x
Σ
Hi @ (maa)) (H, xz,,
E=1
- [10e)) (H, (H,
= (Ha @ x
Hi & (mx Σ
j=1
We now apply the Woodbury matrix identity to compute the matrix inversion in the above display. Recall He = diag(p;)
p.p?". we have
K+1
Σ H. 00 (I.
Σ (exet)
Pa
Σ (p. 00 (p:
Let A E (exe() @ (Ei=1 Puxzz)?) and U E (R(K+))x((tt)) with the first n columns being (p, (iz)) and the following columns (1 @ez)jEle)- Then the term we want to invert is A-U where a block diagonal matrix and can be inverted by inverting each block separately. By the Woodbury matrix identity, we have
(A- -------- ",===""""""""""",=====""""""""""
The gain of using the above formula is significant for large K: instead of inverting the p(K + x P(K + 1) matrix Σi= Hj @ (zy²7) the left-hand side, the right-hand side only requires to invert block diagonal matrix A and (n + p) (n + p) matrix -lapp U A
14

---

12
15.0
Sample quamtles
10
Sample Quantler
10.0
Sample
7.5
classical modem
2.5
modem
modem
2.5
5.0
12.5
Theoretical Quantiles
Theoretical Quantiles
Theoretical Quantiles
classical modeern
3000
classical
3000
classical
2500
frequency
2000
frequency
2000
frequency
1500
1500
1000
500
500
500
0
0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
p-values
p-values
p-values
(a) (n, K] (4000, 2)
(n, K) (5000,
(c) (n, (6000, 4)
Figure $1: The upper row: Q-Q plots of the test statistics from (..3 (in blue) and (115) (in orange)
for different (n, K) and 1000 using وzzz The lower row: histograms of p-values from classical test and our for different (n,, K) and p 1000 using Uzz-
16
12
14
Sample Quantles
Sample Quantiles
12
Sample Quantiles
12.5
10
8
7.5
6
5.0
classical
classcall
2.5
classical
0
0.0
0.0
2.5
5.0
7.5
10.0
12.5
10
15
5
10
15
Theoretical Quantiles
Theoretical Quantiles
Theoretical Quantiles
12
15.0
Sample Quantlies
10
Sample Quantlies
10.0
8
Sample
7.5
classical modem
2.5
classical
modem
modem
2.5 5.0 7.5 10.0
12.5 Theoretical Quantiles
15
5 20
Theoretical
Theoretical Quantiles
(a) (n, K) = (1000, 2)
(b) (rt, K) = (5000, 3)
(c) (n1, K) (6000, 4)
Figure S2: Q-Q plots of the test statistics from 2.. (in blue) and (1.13) (in orange) for different (n, K) and p 1000 using Ωzy- The upper row: covariates are sampled from Rademacher distribution. The lower row: covariates are sampled from distribution of SNPs.
15

---

desr covengeddd.. P% modam covergeecll.7 7%
corets cwmmgee
modem: civengepii35 2%
(a) 9===1
2
(c)
Figure S3: Seatter plot of pairs (a) except using different
SiTAA,I.VAAp with the same data generating process as in Figure
S3 Proof of Theorem 2.1
Before proving Theorem 2.1 present another parametrization of the multinomial logistic regres- sion model. The asymptotic theory of MLE for this new parametrized multinomial logistic model will be used to prove Theorem 2.1
S3.1 Another parametrization of multinomial logistic regression
Recall the symbol "E" is used in (1.5 to emphasize that the minimizer B in (..) is not unique: if B minimizer of (1.5) then B - 61 K+1 is also minimizer of (1.5). for any b E Rp and the all-ones vector 1K+1 k Besides the log-odds model (T.8). here consider another identifiable
parametrization of multinomial logistic regression, whose unknown parameter, denoted by B", is in Rex
Orthogonal complement. To obtain an identifiable multinomial logistic regression model from we consider the symmetric constraint B = K+1 B = 0 as in Zhu and Hastie 2004 K+i ), where the all-ones vector in k=+ Let Q [[K+] be any
thus B B* (1K+1
11"
matrix such that
IK+1
K+1 11 QQ". Q'Q -I-
(S3.1)
We fix one choice of Q satisfying S3.1 throughout this supplement. Let B = B*Q, then B B* R*t and the model (1.4) can be parameterized using 8 as
Plyx = 1|2;)) = exp(x. QeQu)
KE (1, K +1..
(S3.2)
Σ 1 exp(x B*Qew)
The multinomial logistic MLE of B" in S3.22 is given by
δ arg min L. x2),
(S3.3)
where L; R + R is defined by (u) = (u)) for all u E R By this construction, we have Ô = BQ for any minimizer of (1.5). Furthermore, by the chain rule using the expressions (.7) the gradient and Hessian of L, evaluated at DT x;; are
9: = VL, xx) = Q"'g.
H; =BB")
(S3.4)
Throughout, we use serif upright letters to denote quantities defined on the unidentifiable parameter space +1)
B B E (+)( L; L +R R.
Pr. 8i K+I
H, (K++))((++)
and the normal italic font to denote analogous quantities for the identifiable parameter space BoxK,
B B& E R
L 4:
R,
9i E K
H Rx
16

---

Theorem S3.1 provides the asymptotic normality and the chi-square approximation of null MLE coordinates sigh-dimnssions where n. p. > oo with the ratio n/ converging finite limit.
Theorem $3.1 (Proof is given on pagel 27 Let Assumptions 2.7 2.3 and 2.4 be fulfilled. Assume that
either Assumption 2.2
or Assumption
S7.7 holds. Then for any J E pl such that Ho in 1.10 holds,
v3 1/2
9:9
B eg - N0,,1K
(S3.5)
1=1
where V, = H, (H, @ |E Hi @ (ata) " 2 24).
A direct consequence of S3.5 is the x2 result,
V 1/2
E
1/2
Σ
Vi
B"
egl
XK
(S3.6)
The proof of Theorem S3.1 is deferred to Section SS5 and Section S6 In the next subsection, we prove
Theorem
usiin Theorem S3.1 S3.1
S3.2 Proof of Theorem 2.1
We restate Theorem 2.T for convenience.
Theorem 2.1. Let Assumptions 2. 2.4be fulfilled. Then for any E (p) such that Ho in (.10))
holds, and any minimizer B of 1.5
we have
(2
ΣVV, V. BTe, SN IK+1 R+1 )
11d
gig
(2.1)
12
scolar
R(X
R(K ()((+11) Ri++
RIK ()
where g: -v + p; as 1.7 and V; H, (H, x E Hi @ (m @ xz))
The proof of Theorem 2.1 is a consequence of Theorem S3.I To begin with, we state the following useful lemma.
Lemma S3.2. For V; and V, defined in Theorems 2.1 and S3.1
we have V, QVV.Q.
Proof of Lemma S3.2 Since H1 = C2H H2Q, we have
V, = H1 (H,
H: (xt)?) (H, xx)
-H.-
(HHH.))) o (maa))
(,) ((
-H, - (H, 8 (Q 1.)
@ I H. ((a "("0""000000 (I) @ x1
(Q"
=Q² - Q² (H, 3 of
Σ (H, s n )] (H, x
=Q'V.Q.
where the penultimate equality is proved as follows.
Let A = Q & Ip and = Ei= H₂ @ (x,x7) only in the remaining of this proof. remains to prove
(S3.7) Since H,1 = O, we have D = 0. Since Q² 1 = 0 by definition of Q, we have (2 0. If we write the eigen-decomposiiion of D as D = DPK Nunny then u? (1 @ I,) = 0. Hence, with AT ATu..
A(A T D = Dt.
R
A =
Σ
17

---

Since oz =ww uff AA ux - = [u](( 11
2 Ip]ur = uff uv = [[i i'), we have
K+I
pk
pk
A[AT A
u. A
Σ
X
D1,
1=1
where the second equality uses Ao; = AA"u; u; u4 The proof of S3.7 is complete.
Now are ready to prove that Theorem 2.1 is consequence of Theorem S3.1
Proof of Theorem 2.7 By definition of gi and V,, we have 1 gi and 1 1 0". Thus, we have QQ' = and QQ V₁. Therefore, we can rewrite the left-hand side of 2.. (without Vnl2 )j
1/2
1
1/2
- B ej
Σ V.
gig 1=1
1=1
QQ"
g:g
QQ"))
1/2
QQ"ViQQ")B"O
1
n
1/2
Σ
Q9.9
Q²)
Quo")) B"e,
(
9:9
1/2
Q² - viTT)
B"
= o
Σ
969
1/2
Σ
"
Vg
B BTez
where the first equality uses QQ gi and QQTV, Vi, the second equality g: Q" g,
and VÚ = V,Q from Lemma S3.2. the third equality follows from the same argument of (S3.7) and the last equality uses QTQ I and B = BQ
Therefore, Theorem $3.1 implies that the limiting covariance for the left-hand side of 2. is
11" QQ = Ix K+I This completes the proof.
S4 Proof of Theorem 2.2
We restate Theorem 2.2 for convenience.
Theorem 2.2. Define the matrix R (Ix, (x) T E [K++)) using block matrix notation. Let
Assumptions 2.1 to 2.4 be fulfilled. For Å (I.9) and any j E [p] such that Ho in (I.10) holds,
1 [k+ VK + R
11
1/2
Σ V.R. ATe; ee d N(OK (2.2)
g.g.
2z
matris 1e x(K++1
soalar
matrix R((++))x(++)
R
RK
where gi defined in L.7 and V; defined in Theorem 2.1
Furthermore, for the same j E [p],
del n (x( Rij
1/2
M v. RÀT
satisfies 4xx (2.3)
8:g7
ej
Note that Equations
2.2 and 2.3 is stated using 23 ef When E, is unknown, the
quantity Q4z in above results can be replaced by its consistent estimate Giy defined in (..5.. and the convergence distribution results still hold.
The proof is direct consequence of Theorem 2.1
18

---

Proof of Theorem 2.2 By definition of Á in 1.9). we have A= = B(Ix, - -1k)² and
B(Ik: Tikk:lim 11)) (1x: 0x) B(IK+ - -ex- - B (Bex+1)1", which is of the form B h1 with b = +1x+1. Therefore, (Ixx00 0x) is also a solution of 1.5
À (IK, 0k)
Taking B in Theorem 2.1 to be A (IK, 0x) = AR gives the desired x² result (..3) and
Vn 1/2 ((
1/2 1
11 Σ v. RÂ ej - S,N (K - K+1
71
gig
(S4.1)
n
Multiplying (R" (IK+1
1/2 to the left of the above display gives the desired normality K-T
11
result 2.2 by observing (R" (Ik+1
112s Rt) 11/2 (IK +
1g This completes the proof.
<K+141)
S5 Preliminary results for proving Theorem S3.1
S5.1 Results for general loss functions
In this subsection, we will work under the following assumptions with general convex loss function. Later in Section [S6] we will apply the general results of this subsection to the multinomial logistic loss discussed in Section [S3.]
Assumption S5.1. Suppose we have data (Y.X), where Y E (*+((+1) with rows (y1...., yn). and X E Raxx has i.i.d. rows (I) (fi1...,, Fn)) with i; (00, Σ) and invertible Σ. The observations (i)e(() are i.i.d, and y; has the form = f(U;, B*) for some deterministic function J. deterministic B" E IR and latent random variable U, independent of r;. Assume p/n S 8-1
Assumption S5.2. Given data X), consider twice continuously differentiable and strictly convex loss functions (L.)ie()) with each L, : R >> R depending on y; but not on r,-
Provided that the following minimization problem admits a solution, define
B(Y,X) = arg mii Σ La(B'a,).. x.).
BERPX
1=1
Define for each i E [n].
9 X)- VL,(B(Y,X)*'; x.), =,(Y,))= -2LX))))))000,
so that gi (Y, X) k and H,(YY E kx
Define
G(Y,)) = Σ xY,"
i=1
V(Y,X) = Σ (H, & af
Σ H,(Y, @ (maf)) (H, @ x.)). 1=1
so that G(Y,X) E B and V(Y, X) E g If the dependence on data (Y,X) is clear from context, will simply write B, Di. H, G, and V.
C,, Theorem $5.1. Let Assumptions S5.7 and $5.2 be fulfilled. Let C.,, m., m², K be positive constants
independent of n.p. P. Let U* C PRXA
Raxp be an open set satisfying
(1) IF X) E U*), then the minimizer B(Y, X) in Assumption S5.2
exists. H, 1g for each
iE m), 1 Σ H, & (xx²?) 2 (IK E) and m. K
m" IK.
(2) For any ((,)), X)) c U |G X) G(V, 8)11 S L|I(X 1)2-1/211r holds for some positive constant L.
19

---

Then for any j E [p] such that e, B" = 022, there exists random variable & E RK such that
[[((XX)EEUH (G" (G'))- -Ell
C R-R
and P([EE1 > (a)) S for all E (0,1), C is positive constant depending (co, m. m. K, L) only. If additionally X)(( E U*) + then in the previous display satisfies& & +N(0, Ix) and
G) B a (,(
The proof of Theorem S5.1
given in next subsection.
S5.2 Proof of Theorem $5.1
In this subsection and next subsection, will slightly abuse the notations A and Á, which have different definitions than the definitions in the main text.
Let E' = - Sttte be the singular value decomposition of E" where u1.... uK are the left singular vectors and ,,,..., Uk the right singular vectors. If E is of rank strictly less than K, allow some 8 to be equal to that /1/2 3" > kk= Stuuo still holds with orthonormal (un, u and orthonormal (V1, UK). We consider an orthogonal matrix P RRxx such that
Σ 1/1/ ppr = = P
c1, Puu Cp-K+k
Ak E [K].
(S5.1)
Since BB = 0 implies ez ef 1-12 un now on we fix this matrix and consider the following change of variable,
O, we can always find matrix P satisfying S5. From
2 = XX 1/21/22p
A = P5
(S5.2)
immediately follows that Z has i.i.d. N(O,, 1) entries and the first K rows of A" all zeros. Since the response has the expression y; = [(U,, x! ).), Y is unchanged by the change of variable (S5.2) from ZA* XB. We now work on the multinomial logistic estimation with data (Y, Z) and
the underlying coefficient matrix A" in (S5.2). Parallel to the estimate B of B" in Assumption $5.2
we define the estimate of A" using data (Y, 2)
A(Y, Z) arg min
EL(AA' 2.),
AERPXK
where = 2 e; is the i-th row of Z. By construction, we have Å /////2 B, hence ZÁ = x and e. e/ // Furthermore, the quantities depending XB remain unchanged after the change variable. In particular, the gradient and Hessian
(D. 2)) VL,(A² 2))
V'L( x2) V L.
unchanged. It follows that the matrix G and V are unchanged. Therefore, have
of BV G)
o -11
V2i
In conclusion, with the change of variables S5 we only need to prove Theorem S5.1 in the special case, where the design matrix X i.i.d. N entries and the response Y independent of the first p K columns of X. To this end, we introduce the following Theorem S5.2 and the proof of Theorem $5.1 is consequence of Theorem 55.2 as it proves the desired result for ef ÂV -1/2
Theorem S5.2. Let Ca, M. mi, K be constants independent of n.p. Let Z E Rexp have i.i.d. rows (21, 2n) with 2, (0, Ip). Let yı, Yn RKK+1) such that (y1,...,ye Yn) independent of the first p K columns of Z. Consider twice continuously differentiable and strictly convex loss
20

---

functions (L...= with each L; : R - R depending on y; but not on and define, provided that the minimizer admits a solution,
A(Y,Z) = arg min Σ 2,) AERPX i=1
2)) 9e(Y,Z) VL,(AXX,))
V
G(Y,2) Ei= e. 2) E Raxk, and V(Y = Ei - (H, 8 क > H
(xz7) (H, 2,)) R where we dropped the dependence of H; on (Y, Z) for simplicity.
Let C [xx((+1) x 1aax an open set satisfying
If Z) E O, then the minimizer A(Y, Z) exists, H; 3 I for each E [n], c. c.Ip I
Ei= H,(Y, 2) @ (2+2g") ), and m,lk I 1 Ei= G(Y, G(Y, 2) 3 m Ig-
With the notation G(Y, = Σi= e.g. (Y, z),, we have if two Z, Z E Ruxp satisfy {(Y, 2), (Y, 2)} C O then |G(Y,2)-GOX G(Y, 2 S LI|Z - Z||.
For El RP the first canonical basis vector, there exists a random variable E B such that
B[T((Y,)) OJII eTei 1 -Ell
P-R
and P([E512 > (a) s for all E (0,1), C is positive constant depending
(Ca m. m K, L) only. If additioonall P((Y, Z) E 0) - 1, then & in the previous display satisfies (O, Tk) and
of ÁV G) 1/2 IN0,)k
The proof of Theorem S5.2 is presented in Section S5.3
S5.3 Proof of Theorem S5.2
We first present a few useful lemmas, whose proofs given at the end of this subsection.
Lemma $5.3 (Proof is given on page 24 Let NO, o² and F R² - gaak be weakly
differentiable with [+]]( < cc. Let z be an independent copy of z. Then
def Fe) 27 =<(3)
OF(z)
F())) 0²
30 E
Ozi
0zi
Lemma S5.4 (Proof is given on page 25 IfG, G E R satisfy m.lk
G m² IR and
n
m.lk = GGG m² I for some positive constants m, and m². Then
-
I Lin L - G
S L2n Lgn-112||(G Glle:
where L1, L2 are positive constants depending on (K, m,, m²) only.
Lemma $5.5 (Proof given on page 26). Let the assumptions in Theorem S5.2 be fulfilled. Let
Y E [=x((+1) be fixed. If a minimizer Á(Y, Z) exists at Z, then Z + A(Y, 2) exists and is differentiable neighborhood of with derivative
i vee(A) M(g 2 eg + (H,ATez 8 =i),
Ozij
0g. - (Hi 2 M & e, + (H. eg @ 2z)] - 1(1 i)MlA""
Ozij
where M
En= H, @ (2,27 -| immediately follows that
0g:
[, - (H, s 2 MH, @ =))ATe, e, (H, 0 (H, M (g 8.))
aziz
21

---

Corollary S5.6 (Proof is given on page 26 Ei=1 egg we have for each i [n].] pl.
Under the same conditions of Lemma $5.5 for G
de G(GTG)-1/2
Ōzky
G(G" G)-1/-
= oj Av? G) +
(of @ (ef) @ =)) + GfG
Ozij
Now we are ready to prove Theorem S5.2
Let h O - Raxk K be h(Y = G(Y, () G(Y, 2))- In most of this proof, we will omit the dependence (Y, Z) on h, 9,, Hi G, V to lighten notation. By Lemma S5.4 we know this h is 1Layn=1/2-Lipshhiie in the sense that [|A(Y, Z) h(Y, Z) S T.I.n 1/2 Z|| for all {(Y, Z), (Y, 2)) C O. By Kirszbraun theorem, there exists func- R (an extension of h from O to (= *x*)) such that (Y, Z) h(Y, for all (Y. Z) O, [H 1)op I and || Z) HY, Z) || S
Proof of Theorem $5.2 2
tion H (k+x) x R2xp
LLnn 1/21 F for all {(Y, 2), (Y, 2)} "X gnxp,
For each j E [p]. let 2, = Ze, be the j-th column of Z to distinguish it from the notation 2,, which means the i-th row of Z. Let 2 NO0,1 In) be independent copy of each columns of Z, and
Z (Ip ) + z That is, Z replaces the j-th column of Z by ž. By definition, Z1 1
and 2 1 21
Let ę = [H 2-), E " then |&| < | || since [H 11) < 1. It follows that
P(IE)| > x (o) < P( 1| (a)) =
Note that the first K columns of Z are exchangeable, because they are i..d. and independent of the response Y, we have for each E [p- K].
B[[((Y,) G) - =|1²1
E
zi H(Y (Y,2")) 21
=E((((, A
H(Y 20)13
where the last line holds for any E 1p- K] because (z ef A, 2) (ze, et et A, 29). Therefore,
H(Y
2 =
K
H[, [[(Y, (((())))00|| ÁV(G G) +2 H(Y, z/) /1²
p
K
P- K
T((Y) =
def G)-1/2
+ 2 H(Y, 2タ) - Remz [']
K
Ozi
t=1
where Rem, = 8 (G
Te? T BIG 1/2
from Corollary S5.6
0zs
and M 1> Rh (H, @ ) from Lemma 55.5 Using (a (a+ b)2 < 2a 202, the above display can be bounded by sum of two terms, denoted by (RHS); and (RHS)
For the first term,
(RHS):
2
Σ
K =
2[(1(Y,))
Σ
de deT 7)
+ z HY, 29 17.
D K
2zZt
22

---

Let Let F(zz) (( + = H(Y, 2), then F(2) = H(, 2¹). Apply Lemma S5.3
to F(zz) conditionally on ((1 ete, ) we obtain
E E[T((Y, E 20 =)
Σ
deg h(,,
+ zz F(2) 1"1
Ozit
Des F(zi) - 2 r(a)l"]
=[(((>,)) E =0 F(z)) - Σ
Ozst
E F(z))
F(z))
deg F(z()
2. 1(2)|
Ozit
S 3 Σ R (F()))
Ozit
(Y Z)
12
3
Σ
R
Ozie
F²
where the first equality uses hY, Z) = 0 from the KKT conditions G = 0 and h(Y, = G)G"G)- It follows that
OH(Y, Z)
(RHS) - p-K
K
Ozir
F
t=1 ==1
Note that the integrand in the last display is actually the squared Frobenius norm of the Jacobian of the mapping from Raxp to R= + H(Y, 2). This Jacobian is a matrix with nk rows and up columns, has rank at most nK and operator norm at most LLzn -1/2 because Z + H(Y, Z) is 1/l2n11/22_ppp Lipscchtttz from Lemma $5.4 Using
|| I (ank(A)]AAl |lpp. we obtain
(RHS) 6K (K(LL))²/(> K).
For the second term (RHS)2
2 p-k
Σ [[(((,,) 2) E O}| Rem, [|²]. By definition of Rem
and (a + b)2 < 2a² + 222, we obtain
k
(RHS)2 <
4
E
[[(((Y)))
E cO
Σ
Cff & of MHH, G) 1/2
(S5.3)
p K
8=1
1=1
4
P-K
K
Σ
G) 1/2 2
(S5.4)
K
Ozit
1=1
We next bound (S5.3) and (S5.4) one by one. For (S5.3). we focus on the norm without (GT G) -1/2 which is || EE= (g @ ef MH,, @ 2:)l.. With ||a]| m|xl:|| a" u in mind, let us multiply to the right by unit vector E RR and instead bound
Σ oo Se Se MH = Tr T E(HNOa)s < KII(K llopl EH llop
because the rank of the matrix inside the trace is most K and Tr|.1 S K|| 1lop holds. Then
Gllop S |1 Σ (H.uBeeee De (HHuo)))) llop x8 38 n²) Σ (HHuee))) Gllopp
Next, (H @ cee ) lop (H, & c
(u @ In)llop < 1 because H, - I and ||u||| 1.
In summary, the norm in (S5.3) is bounded from above by
G) 1op M op 2 GG
To bound (S5.3). since in the event (Y, Z) O, m.lk 3 L 3 m and 1 E (H, & 2127T) ž c. C.I we have ||G|lop S v Vm"n. n. hence || S c-l Thus, the above display can be bounded by
n
K(m.n) -1/2, |1 lop V (m" 1/2c 1lop.
23

---

Since Z E RTxP has i..d. [0, 1) entries,
Davidson and Szarek 2001 Theorem 11.13] implies that
Z op S vn + /P= - n/n. Therefore.
S5.3 S C(e,, K,C,,KLU L)m
Now we bound S5.44 Since
K
IG G) 1/2
Ozit
==1
0(G G) 1/2
Ozij
1=1
K
0G ((GTG))-12 G1-11/2
G
ex!
Ōzij
K
OGT G)-
)
Ge
eN
Ozij
n
G)-1/2
|
Ozyj
Using | <
nK, and the mapping Z (G* is LLan-1 Lipschitz O using we conclude that
Lemma S5.4
4K LL1 K). S5.4
Combining the above bounds on (RHS)1 and (RHS)2, we have
[[(( 2) ( G)) - =|1]]
C(c,, K, m. K, m², L, L L2)
(S5.5)
p-
where the constant depends on (c., K, m. m², n L)) only because In and L2 are constants depending on (K, m., m²) only.
If additionally 2) E 0) + L, we have P((Y, E O) 1 using (Y,Z) d (Y, 21). Therefore,
|I(YY, (1) 0) - [M 24) 21) 2') -) (S5.6)
By S5.5 know (G' -)-- 1
- = when P((Y E 0) 1. Hence, we conclude
(/22) A d 40, and 1/ A e, 1| +x
We next prove Lemmas S5.3 to $5.5 and corollary S5.6
Proof of Lemma $5.3
Let 20 (z TT 22n, then N(0, (I For each l [K], let
f(k) 22n
+ IR 2m
be
(z0) ([((2) - F(E)Jex
On
so that 2/ (20) 2T)F(2) F(k) and div ((k)
-=
8e F(g)es
Applying the
second order Stein formula [Bellec and Zhang 202 (see also
Tan et al. 2022 Lemma F.1] for
24

---

collection of such formulas) to f(() gives, with Jac denoting the Jacobian.
T F(z)e 0² de F(z)ex F())(() '1
Oz
[ (th - o² div ((8) 91
E
o'BIS (zn) )| + +o ((+ (20))
(FE) P(8)]ek| o EE (J)c[F(z)ex]
( - -Jac[F(3)ex)
'1
Onxn
Onxn
20² - BF()))l + (TTTa(/))))T((((((()) I 30'E| Jac(F(z)
where the last incquality uses the Gaussian Poincaré incquality, and the Tr(A²) S |- Summing over [K] gives the desired inequality.
incquality
Proof of Lemma S5.4 We first prove G GT GTG Lipschitz by noting
[[G - O (G-
s |G - Gllool + ||G|lopp
2v Wm Gllop.
Then we show GTG + GGT is Lipschitz. Let A = GT and À ÕT we have
||- - À - 1lop
||A - A)
< || - Al|op - |lopl Ā --11
op
S (m.n) 3|1 Ållop.
We next prove (GT () 1/2 is Lipsehitz. Let S =
and if
with ||u||| = 1 is the eigenvector vS V with eigenvalue d, then
o u=(S - (5 = n" (V5 V55 of VT
du- VSu du Su
du duss V Ju.
As d can be chosen as +IIV V op inequality (van Hemmen and Ando 1980))
(this argument is a special case of the Hemmen-Ando this implies
IVSS - V5l = |u S)u)
S 1S-511 1lop
S $1
5
u (T(5 v) Amina ( v/$+V5
2/ n
Combining the above Lipschitz results, we have
m²
- ( -G llop I (m"n) (2(mn))))))))) (mm GGlopp m
n -IG Gllop
immediately follows that
-16 S VK m² m I Gl
That is, the mapping G + (GT 1)- is L Lipscchitz where L1 V 2m"
25

---

For the second statement, the result follows by
- lop
< || +|0
- G)-1/2 -111) op
s |G G (m)
Ln Gllop
Hence,
VK (mmm 1/2 + 1/3|G - GlIe;
where L2 VV (m= 1/2 + (m1))
Proof of Lemma $3.5
Recall the KKT conditions Ei= "9/ = 0p OpxK. We look for the derivative
with respect to Zij. Denoting derivatives with dot, find by the chain rule and product rule
Zz
& kzij
I(1 = i)ej,
0gi
0gi
OA
Hi[AT= 2 + I(I= i)A e,].
91
kzsj
TAT Ozij
Thus, differentiating the KKT conditions w.r.t. Zij by the product rule gives
Σ 1=1
o) + + 21 + 1(1 = i)ATe) e,) " Hi] = O.
That is,
e,gi
Σ
पकर AH + zee AH = 0.
1=1
We then move the term involving À to one side, and vectorize both sides,
9: 2 + (H, (=))= Σ H1 8 2T²7 vec(A).
With M = = (H) 2 -1 obtain
vec(A) M e, (H, e, @ 2)).
Hence, using veo( H 21) = vec(= AHH = (H) @ of veol A) gives
ÿ1 = (Hi @ 2 veel A) + i( ej
- (1l of M 99 2 8,, (H1,AT e, 21) + 1(1= i)HliA"e
Thus,
02) {,
0: = -H,, 8 xf
-(H, @ of (
(H, 2 &c,
[H, (H, s 2 M &
V.A"' - (H, 8 M 0c,)
where V; = -H; - (H, & 2! MHH 8 2.).
26

---

Proof of Corollary 55.6 For each i E [n].j E [p], we have by the product rule
def
G(G G)
1/2 11/2
Ozij
Ouf (G' G) G) e T G)
dzij
Ozi
-Vv'A'e (H, @ o MM e,)l (G" e C
0G" G) 0GG
Ozij
3 () -)1
=v [- e & G)
1/2 + T G
0zig
With V =EE V₁, we further have
a deL G)G²
Ozig
G) Ozij
of AV. (G G)- -1/2 Σ (of @
o
S6 Proof of Theorem S3.1
Recall that Theorem S5.1 holds for general loss function L; : R + R provided that conditions (1) hold. In this section, consider the multinomial logistic loss function L; defined in Section S3.1 To be specific,
and (2) in Theorem S5.1
K+1
K+1
[,(u)= = Σ
Yuee Qu + log
Σ k==1
exp(ee
Vu E IR..
(S6.1)
In order apply Theorem S5.1 we need verify that, when L; in (S6.1) is used, the two conditions (1) and (2) in Theorem [5.1 hold. To this end, we present a few lemmas in the following two subsections, which will be useful for asserting the conditions (1) and (2) when apply Theorem SS.1 to prove Theorem S
S6.1 Control the singular values the gradients and Hessians
Before stating the lemmas that assert the conditions in Theorem S5.1 define
U = XY, E ( x Raxp, - exists, [X
11 K M1
Ug {r e [2)((+1)
[yya = 1) 2 ym for all [K [K+1]]
==1
Lemma S6.1 (deterministic result on gradient). Let L; be defined as in (6.11. Assume that either Assumption 2.2 or Assumption S1.7 holds. IfY Uyr any M BaaK such that |MQ S nt.
we have
m. m. 3 n Σ (L oT) KII.
where m. is a positive constant depending on (K, (K,Y,T) only
Proof of Lemma 1 Without loss of generality, let's assume that yn is an integer. Otherwise, we can replace it with the greatest integer less than equal yn, denoted [yn].
Y E Uy there exists at least yn many disjoint index sets (SI..... Syn such that the following hold for each I hyn].
(i) St [n]: (ii) |S|| = K + (1M1
Σ Yik 1, iES:
Ak E [K 1].
27

---

Since S are disjoint and U St C [n2], we have
M
IQM eill S Σ @M
11Q < nT.
1=1 FE
It follows that at most an many of / E {1, yn s.t. Eies, IIMMT eill > r/a, otherwise the previous display can not hold. In other words, there exists a subset L c {1, 2,,,,.., .... yn} with | 2 (7 a)n s.t. Eies, [|QMTe. T/a for all / E L Define the index set
UJ. then |4] 2 (K + 1) a)., and IQMee:lla VT/a for all E I. Let us take 7/2, then [ 2 n and T 2 <( + 1)n/2. Recall that [u(u) = L.(Qu). we have VI,(u) = Q" (Quu).
Thus,
VL = Q'VL(QMM = (-yx +p))
where Pt E RR++ and its k-th entry satisfying
Pik
Σ K+I exple QMTe)
E 1e.1-d. c.
(S6.2)
for some constant c E (0, 1) depending on (T, a K) only. Therefore,
n
( e4 VI. e2)
n
(ys p:) TO
m n
(y; p:)(y; P.
1=1 IES:
2 Q" Σ
(y; pi)ly; Q.)
TeL' VES)
n Q2" Σ Af A.Q.
where A₁ [[(*+))((++) has K + rows (y: - i E Si). We further note that A, is of the form (Ix+1 Pt) up rearrangement of the columns, where P, E ((++))((++) is stochastic matrix with entries of the form
exple exple[QMTe) QMT K+I Σ exp(e
i [K(++]
By (6.22), for each / E L", the stochastic matrix P; is irreducible and aperiodic and ker(Ix41 Pt) the span of the all-ones vector 1. Therefore.
11 (IK+1 - (IK+1 Pi)(T41- R (JK+1 Pr).
follows that
= -P)) raaakk( P()QQT) rank( - P.)Q) s K.
We conclude that the rank of (1K+1 P K.
28

---

If P denotes the set of matrices {P ((++))((++)) stochastic with entries in [c, 1 - d],, and
Sk-1= {a E RR ||a|| 1}. By compactness of P and s we obtain
VI,(M e))
Amin (Q" A(Q)
LEL"
Σ let.
min a Q AZ AQQa aESK-
at Q² (IK41 - (Ix+1 - P)Qa
min aESK-1,PE
n
7 0, Q (Ix+1 PP.) (IK+1 PP.)Qa,
0² (Ix+1 -P)) (IK+1 P.)Qa.
2
:=
where a, E SSK1,, E P, and m. a positive constant depending on (K, 7) only. The first inequality above uses the property Amin (A B) Amin(A) + Amin (B), where and B two positive semi-definite matrices.
VL, = m. m. For the upper bound, since | S 1 by definition of Q and all the entries of (y, -p.) p.)" are between -1 and 1 if Assumption 2.2 or Assumption SI.1 holds, we have
In other words, 1 E
c (.VL(M&&., e 112 > (y: p.) P.l S nK.
Lemma S6.2 (deterministic result on Hessian). Let L; be defined as (S.. For all i E [m], we have L2L,(w) IK for any u E R and
V2Li(u) c.lk:
min
WERK
Q
Sr
where c, is positive constant depending on (K, r) only.
Proof of Lemma S6.2 Recall that L,(u) = LQu),, we have
where V2L,(Qu) diag(p:) - P.Pp" and the k-th entry of p. E k+ is defined as Pik exple Qu)
EKE exple, Qu) for all i [m]. [K 1]. Thus, Pik 1 for all [n], [K 1].
For any vector E 5k-1 we have
T a =a T² [diag(p:) - piPp'1Q
-a (pag(p;)QQ
< = diag (p, p1lppp
$1,
where the last inequality uses || < and Pa - for any k. Hence, V2L;(u) = Ig for any i [n].
Now we prove the lower bound. For any u E R such that |Qu|loo S r, we have
exple (qp) E Ec.11-0 1
Σ
K+1 Σ == exp( c Qu)
29

---

for some constant depending on (K, r) only. For any vector SK-1 let 17 = Qa E RK+1, then
1 I r 0 and
Li a² Q" (diag(p))) (p))
1a
p.P
_rd(diaa(p:) - P:P?
K+1
K+1
Pixnn
Pix:TT))
k=1
2 Puenh
Σ Σ
7
Pike
Σ
Pu-n (1
Σ Pak)
=0,
where the last equality uses EKET Pik 1, and the inequality follows by Cx & pam)) (Ex vPa pPt-m)?) S Ex Pik Σx Pa using the Cauchy-Schwand inequality, and here
holds if and only if /Pik /Pisill: for each k, which is not true since Pik E [c, c] and 1 7 0.
Let H (Q" (diag(p) pp')) PE
}. then H is compact and min aESR HEH
min (ERK Sr
Amin (V2((u)) 2
a²Haa a._ H.a, >
Sr
for some a, E SK-1 and H. E H. Therefore.
min HERK Qullse
V Li(1)) 2 c.lk
where c, is positive constant depending on (K,r) only.
S6.2 Lipschitz conditions
We first restate the definitions of following sets,
U = {(Y E (R=(((-1) x R : exists, | [B
KFI )1b < <nr)
11T
U, [
Σ T(ya = 1) 2 7m 7m for k E [K + 11}
Lemma $6.3. Assume p/n < 8-1 < -1- a for some E (0, and 8
Let
[n]
[ = [n(1 } and P = EvEl c.e Define
U, = (X RRPP, min Amin( IEI
|XE 1/2 3 12X 2 62 Vn
< (²}
(S6.3)
12
for some positive constants Φ., 8², which depend on (8, a) only. Let U* = ((, :
Uy X U.}. Then under Assumptions 2.1 2.3 and 2.4 and if either Assumption 2.2 or Assump- tion SI.1 Thols,, we have
() -
(ii) Let G be defined in Assumption S5.2
If (Y, C U*,, we have
X|G(Y,X) - -(Y S LII(X -X2E-1/111,
where I. is positive constant depending (K, T,) only.
Proof of Lemma S6.3 We first prove statement (i). Under Assumption 2.1
Bellec 2022 Lemma
7.7] implies
((min Amin ( 1/2xx PXXE
)?)
30

---

for some positive constant o, depending on (d) only. Furthermore, Davidson and Szarek 2001 Theorem II.13] implies
P( ||XX-1/2le ~Φ*) +1 vn
for some positive constant 6² depending on only. Therefore, P(X E U2) + 1. Under Assump-
tion 2.3 we have P(Y E Un) 1. Under Assumption 2.4 we have X( U) + 1. In
conclusion, under Assumptions 2.1 2.3 and 2.4 we have P[(,, X) E U*) + as a,
Now we prove the statement (ii). For a fixed Y. let
X) U B, B be their corresponding minimizers of (S3.3), and G, G be their corresponding gradient matrices. We first provide some useful results derived from the KKT conditions. From the KKT conditions ==G= = 0, we have
=
- - -B) -------------------------------- -B,
=
Since [v)llp S 1 for any E RK from Lemma S6.2 (L,(-) is I-Lipschitz. Thus,
(XB - G 6) = Σ
1 x BT 24 VL, x.) - VL,(B"I)))
2
VL((',)))))). -V((()))V
==1
-1G-
If (X,X) e U, we have
IX P/ + |X BQT - 2nt. That is,
Σ 11QB x4 12 + 1Q < 2nT.
Define the index set
I= fie [n| |Q n, 1 IQB م 2
then we have |/| 2 (1
a)n by Markov's inequality. Thus, for all i E I, we have QQQ x;11 V
i: I
V
2r
a
Applying Lemma S6.2 with r
2r gives min
L'L 2 c.lk:
[|QuI>>S 2x
where c, is a constant depending on (K, a). Therefore,
c. (XB
=c Σ 1 iTā,l²
iE!
B x,- B" n,
IET
<
(BT xi BT n., -L
i=1
=XX -XB-XXB XB,G-
+ ((X - X)B.G- ((-)))-G( = -(X xx) B))
We next bound the first line from below by expanding the squares,
[P
8
x B3) Px x) B1
2 ( - B>I (2 (B)
2 n + -2( -x))B,
31

---

where in the last inequality we use the constant 6- in (S..). Therefore, we obtain
c, B)1?
((- ((--) -2. - PXX 8)B).
Together with the inequality that | 1 (XB - XB, G), we obtain
c032 / - B)1 + |G 3
S 5-2((X - x) - B) + 2((X - x)B.C G G) - -20, - B), P((X 8)
5(4+ |I(X - XX)--11 op (1 - B)llF V || G(lF) - (11 B|1 ||
[G]]lom) vn. vn
Vn
where we bound (X(B PrX X)B) by definition of d", (R(B- -
- P/X-
-(EV)(B-BB) 137 B)B3
(B>) ())B--))))()PX-)))"" op ||x - XX-=-111 1/2 |opp|>> Blr
SV B)) - 1/2=-1121lop ||:11/2 Bl|l.
S / To this end, since Φ |E s ||P:XB||=/vr S || F V |X ||//
Now derive a bound of the form
S VT.
Applying Lemma S6.1 to M xB, we have - gi = m.lk. Therefore,
Tr(g.g 2 Km,. |GG| = 1 Σ 1 =
n
This implies that
3 B } STS
|
Km. (I)
In conclusion, if ((Y,X), C U then
( || G - Cn 1/2 1 1/= S CK][(X X)E 1/2
where C constant depending (K, 7: T, a) only. Note that |G]|e S Vnk since all entries of G are in [-1
Lemma S6.4. If p/n S 5-1
< <1 a) for some a E 10, and 5 1. If (Y, E and X Un.
where U2 is defined in Lemma $6.3
we have
Σ H. & (I, = 01 (Tx @ E),
where ei is a positive constant depending on (K, at )b) only.
Proof of Lemma $6.4 If (Y, E U, we have |X I 1 S Nt. Define the index set I= {i E [m] : 1Q + 5}
then we have I (1 an by Markov's 's inequality. Therefore, for any E I, l[Q8 1//000 S T
Applying Lemma S6.2 with u = 8T I, and r = 2, we have for any i E 1, H; = V2L,( ax)) >
T, a) only. Therefore, (Y,X) E and c. K for some positive constant c, depending on (K,
X E Uz, we have
H. @ (I 2 - Σ
H. 8 (x
n iEl
m c. n iET
1k (Z,47)
CL @ X
- cI 8 2
cI ② E),
32

---

where P1 = ExEl ege? and is positive constant depending on (K, T, a, 64) only.
S6.3 Proof of Theorem S3.1
The proof of Theorem S3.
a direct consequence Theorem S5.1 by noting
Vn 1/2
9:9
-1/2
B ej = G)
which is a consequence of the identities G = Ei=1 e. and V = - V.
It thus remains to verify the conditions (1) and (2) in Theorem S5.1 from the assumptions in Theorem S3.]
Applying Lemma S6.3 with a chosen as 1 - 8-1/2, we have for (V,)) CU.
XG(Y, G(Y, X)((= - LI(X x)--1/111,
where L is positive constant depending on ?K 6) only.
Apply Lemma $6.4 with the same 5-1/2, we have for (Y,X)= E U*,
H4 (I ž c+ @ 2), Σ 1=1
where c, is positive constant depending on (K, T, 5) only.
Applying Lemma S with M =X,, we have for (Y, E U*
m. I n VLMM & = KIK.
where m, is positive constant depending on 7, only. Therefore, the conditions (1) and (2) in Theorem S5.1 hold when the multinomial logistic loss is used. This completes the proof of Theorem S3.1
S7
Other proof
S7.1
Proof of Equation 1.12
(Classical asymptotic theory with fixed p)
Here we provide derivation of the asymptotic distribution of MLE under classical setting, where fixed and n tends infinity.
We first calculate the Fisher information matrix of the multinomial logistic log-odds model [18
with covariate x ~ No0, Σ) and response y E (0, 1) one-hot encoded satisfying KK+
Note that the model 1.8) can be rewritten as
exp(x A
Plys 1/z) =
Vk K}
1 ΣE=1 exp(x" Atee))
P(y*+1 1|x) + ER=1 exp(xT Ateyr)
The likelihood function of parameter A BRXK is
K exp(x" Aex)
9K-1
L(A)
k=1 W exp(x* Aew)
Σ K -- exp(x* Aew)
K
(exp(x T Aex)]
1+CK
1 E exp(x Aex)
33

---

where we used K+1 = 1. Thus, the log- leelikkllioodd function is
((A) yua" Ae log 1 exp(x Aew)
It is more convenient to calculate the Fisher information matrix on the vector space Rek instead of the matrix space IRP* To this end, let vec(A"), ), then x" Aee ex A x (x" e vec(A"")
(x of )0 and the log-likelihood function parameterized by 0 is
K
K
{(O)
Σ yx k=1
& Ee
log[1 +
exp((x" 2 el
By multivariate calculus, obtain the Fisher information matrix evaluated at 0" = vec(A*
I(0") = -E á ae(o)
50 00
0=0-
= BB((x' " & (diag(T") x π x ")},
where n² E R with k-th entry क
-+ exp(2 Ategr) for each k E [K].
explx " A
From classical maximum likelihood theory, for instance (an der Vaari 1998 Chapter 5], have
vn(A - 8*) (0,
where i = voc(AT) and À is the MLE of A Furthermore, if the j-th covariate is independent of the response, we know e. A" = O7, then
VnA ej = vn(A"e, e, T,) vn(ez @ [k) 8*) (N((,,))
where
5 (o) I I (e, Ix)
=e] cov(z) e, [(diag(TT) - -
holds by the independence between the j-th covariate and the response under Hg in (L10). This completes the proof of Equation (1.12)
34


Reviews:
[{"review_id": "LPpU0MqiUS", "reviewer": "Reviewer_FsXG", "summary": "This paper studies the asymptotic distribution of the MLE of the multinomial logistic regression model when the sample size and the number of parameters are of the same order.  The validity of the asymptotic theories is evaluated through extensive simulation studies. The paper is overall very well written and the presentation is clear.", "questions": "1. Related to my question in the weakness section, one could use the proposed tests to exclude some noise covariates in the multinomial regression and see how much improvements can be achieved in terms of classification errors. Preferably, one should compare its performance to some existing methods. \n\n2. I would strongly suggests applying the proposed methodology to some real-world benchmark data set. The statistical tests may provide additional insights and the classification errors can be evaluated through cross-validation.\n\n3. The paper's theoretical framework appears to heavily rely on the seminal work of Sur and Candes (2019), because of which I still have some reservation on the novelty of the theory. To address concerns about the novelty of the theory, it would be wise to explicitly highlight the unique theoretical challenges posed by the multinomial logistic regression model compared to the more commonly studied binary logistic regression model. ", "limitations": "NA", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is very well written, and the problem under consideration is of great interest. Theoretical results are sound and the numerical experiments are sufficient.", "weaknesses": "One of the motivating example of the paper is to study classification with 3 or more classes. And one of the most important goal in classification is the classification error or AUC. In my opinion, the impact of the proposed method on classification errors should be evaluated through simulation studies and some real data examples. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the insightful and positive feedback. \n\n> **Q1**: Related to my question in the weakness section, one could use the proposed tests to exclude some noise covariates in the multinomial regression and see how much improvements can be achieved in terms of classification errors. Preferably, one should compare its performance to some existing methods.\n\n**A1:** \nThe tests we propose provide strong guarantees for controlling Type I error: rejecting a null covariate (i.e., wrongly making a discovery regarding an insignificant covariate) happens with probability approximately $\\alpha$ for any fixed $\\alpha\\in(0,1)$. However, this does not guarantee that the features that are not rejected should be safely removed from the model (which would require Type II error control). There could be instances where $B^*$ has many small coefficients on all features, so small that the test would not reject $H_0$ for all features; but overall the MLE $\\hat B$ may still be more accurate (say, in prediction error) than random guess in such cases. In such situations, removing the features for which $H_0$ was not rejected would lead to a model with dimension 0 (with prediction probabilities $1/K$ for all classes) and would perform worse than the MLE for estimation of $B^*$ or prediction on a test set. In short, while the method suggested by the referee could be implemented, there are situations where using the MLE on the full model performs better.\n\nTo our knowledge, there is no existing testing method with the same goal as ours for classification problems with 3 or more classes. \n\n> **Q2**: I would strongly suggests applying the proposed methodology to some real-world benchmark data set. The statistical tests may provide additional insights and the classification errors can be evaluated through cross-validation.\n\n**A2:** \nThanks for the suggestion. We conducted a real data analysis by applying the proposed test to heart disease data from the UCI Machine Learning Repository. After standard data-cleaning processes, the dataset has 297 instances with 13 features, including age, sex, and other attributes. The response variable was transformed into 3 classes (0, 1, and 2) after converting 3 and 4 to 2. To demonstrate the validity of the proposed significance test, we generated a noise variable from a standard normal distribution, resulting in a dataset with 297 instances and 14 variables. We applied our test to this noise variable using both the proposed test and the classical test, repeating the experiment 10,000 times. The type I error of our proposed test is 0.0508, aligning well with the desired type I error of 0.05. In contrast, the classical test exhibits a type I error of 0.0734, significantly exceeding the desired rate of 0.05. These results confirm that the classical test tends to include noise variables, leading to false discoveries. On the other hand, our proposed test is more reliable and achieves the desired type I error.\n\nRegarding\n> The statistical tests may provide additional insights and the classification errors can be evaluated through cross-validation \n\nwe refer to our response to the previous question, where we explained that the proposed method cannot be readily used to safely drop features from the model. Our method provides a strong guarantee for Type I error control, but this does not readily lead to improved classification error on unseen data.\n\n> **Q3**: The paper's theoretical framework appears to heavily rely on the seminal work of Sur and Candes (2019), because of which I still have some reservation on the novelty of the theory. To address concerns about the novelty of the theory, it would be wise to explicitly highlight the unique theoretical challenges posed by the multinomial logistic regression model compared to the more commonly studied binary logistic regression model.\n\n**A3:** \nThe proposed theory is novel because it is different from the seminal work of [SC19] in at least three aspects: \n\n(i) The results in [SC19] only apply to binary classification problems, while our proposed theory also applies to classification problems with three or more classes. \n\n(ii) The techniques used in the proof are significantly different from those of [SC19]. Case-in-point: even for $K=2$, our argument does not rely on or involve the nonlinear system with 3 equations and 3 unknowns which is central to [SC19] for binary classification (or its generalization to $K\\ge 3$ classes). The Approximate Message Passing proofs used by [SC19] are also not employed here. Our result for asymptotic multivariate normal and chi-square distributions is based on applying Lemma S5.3 of the supplement to a carefully chosen function $F(\\tilde z)$ valued in the set of matrices with all singular values equal to 1, which is a new result of independent interest to exhibit asymptotic chi-square distributions in such problems. This further poses unique challenges to ensure that the derivatives are well-behaved despite the normalization requiring all singular values to be equal to 1, and we solve some of these challenges specific to the cross-entropy in section S6 of the supplement.\n\n(iii) From the computational side, all the quantities involved in our results can be readily computed from the data, while the relevant quantities in [SC19] require the estimation of the solutions to the nonlinear system mentioned in (ii). For $K=2$, this is not an issue, as [SC19] proposed the ProbeFrontier method to estimate the solutions to the nonlinear system. \nAny extension of their work to 3 or more classes (which, as far as we are aware, does not currently exist for the multinomial logistic model) would also require solving solutions for a new nonlinear system. We expect the ProbeFrontier method to break down since the new nonlinear system would now include several scalar parameters due to the matrix nature of $B^*$ (as opposed to the single unknown scalar $\\gamma=Var[x_i^T\\beta^*]$ in [SC19]).\n\n[SC19]: Sur and Candès (2019), PNAS"}, {"review_id": "tqpPu6UK76", "reviewer": "Reviewer_7JUy", "summary": "In this paper, the authors provide asymptotic characterization of the behavior of the maximum likelihood estimator (MLE) of multinomial logistic model (with more than two classes), in the high-dimensional regime where the dimension and the sample size of data go to infinity at the same rate.\n\nUnder some technical assumptions (that may need some further elaborations), this paper develops asymptotic normality and asymptotic chi-square results for the multinomial logistic MLE on null covariates (see Theorem 2.1 and 2.2). The proposed results can be used for statistical inference and significance test of some specific features. Numerical experiments on synthetic data are provided in Section 3 to validate the proposed theory.", "questions": "1. While I am almost fine with Assumption 2.1 and 2.2 (just being curious, is some upper bound on the spectral norm of the covariance $\\Sigma$ needed? or it is just a matter of scaling with respect to $n,p$?), I am a bit confused by Assumption 2.3 and 2.4: are they something intrinsic or for the ease of technical analysis? What happens if, say Assumption 2.3 is violated? Can we have something similar but just more involved or the MLE is totally different? Also, Assumption 2.4 is a bit misleading, in the sense that the assumption is not instinct, and should perhaps be reduced into some assumption on the dimension ratio $p/n$ and/or statistics of the data? I believe it makes more sense to assume something like \"the dimension ratio $p/n$, covariance $\\Sigma$ and xxx satisfy that xxx\". I am also confused by the paragraph after Assumption 2.4 and I am not sure the convergence of some multinomial regression solver can be used as a rigorous theoretical indicator. The algorithm may converge (or believed to converge) due to many reasons. Perhaps some better (numerical) criterion can be proposed by, e.g., checking the gradient and/or Hessian of the point of interest.\n2. for the sake of presentation and use, it be helpful to present Theorem 2.2 and the estimation of $\\Omega_{jj}$ in form of an algorithm.\n3. Almost nothing is mentioned for the proof of the theoretical results (Theorem 2.1 and 2.2): is the proof technically challenging or contains some ingredients and/or intermediate results that may be of independent interest? Could the authors elaborate more on this?", "limitations": "This paper is primarily of a theoretical nature, and I do not see any potential negative societal impact of this work.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "This papers focuses on the fundamental and important problem of MLE of multinomial logistic model in the modern high-dimensional regime.\nThe proposed theory improves prior art in characterizing the asymptotic normality and asymptotic chi-square results, both of significance to statistics and ML. \nThe paper is in general well written and easy to follow.\n", "weaknesses": "I do not have strong concerns to raise for this paper. \nSee below for some detailed comments and/or questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable and encouraging feedback to help us improve our paper. Here we provide answers to your concerns point by point.\n\n> **Q1:** While I am almost fine with Assumption 2.1 and 2.2 (just being curious, is some upper bound on the spectral norm of the covariance $\\Sigma$ needed? or it is just a matter of scaling with respect to $n,p$?), I am a bit confused by Assumption 2.3 and 2.4: are they something intrinsic or for the ease of technical analysis? What happens if, say Assumption 2.3 is violated? Can we have something similar but just more involved or the MLE is totally different? Also, Assumption 2.4 is a bit misleading, in the sense that the assumption is not instinct, and should perhaps be reduced into some assumption on the dimension ratio $p/n$ and/or statistics of the data? I believe it makes more sense to assume something like \"the dimension ratio $p/n$, covariance $\\Sigma$ and xxx satisfy that xxx\". \n\n**A1:**\n\nThe upper bound on the spectral norm of $\\Sigma$ is not needed in our theory thanks to the rotational invariance: a key step in our proof involves rotating and scaling the covariate with covariance $\\Sigma$ to the identity matrix $I_p$ (see page 20 in the supplement).\n\nAssumption 2.3 is indeed necessary for our technical analysis. One intuition that this is a fairly mild assumption is that if for some label $k_0$ we observe $o(n)$ observations, the dataset obtained by dropping all observations with label $k_0$ still has sample size $\\tilde n = (1-o(1))n$ because we dropped $o(n)$ observations. We end up with a dataset with one fewer label, and almost the same sample size. \nFurthermore, if $y_i$ for $i=1,...,n$ are iid from a fixed distribution independent of $n,p$, if $P(y_{ik}=1)>0$ then by the weak law of large numbers, the probability to observe $\\gamma*n$ responses $y_i$ with label $k$ converge to 1 for any $\\gamma\\in(0,P(y_{ik}=1))$. On the other hand, if $P(y_{ik}=1)=0$ then this label $k$ may be removed as it will not appear at all. \n\nFor Assumption 2.4, we agree with you that imposing conditions on the MLE $\\mathsf{\\hat B}$ is not intrinsic. It is a reasonable assumption to us because (a) the event in Assumption 2.4 can be observed directly, and (b) because the goal of this paper is to study the property of $\\mathsf{\\hat B}$ provided it exists (it would not make sense to study such tests based on the multinomial MLE in settings where this MLE does not exist with high or positive provability; nothing can be done with MLE in the event that it does not exist). \nWe suspect Assumption 2.4 could be reduced to \"the dimension ratio p/n, covariance $\\Sigma$ and xxx satisfy that xxx\" as suggested by the referee, and as Candès and Sur (2021) did for binary logistic regression.\nHowever, such analysis for 3 or more classes does not exist yet (only partial answers are known, cf. lines 326-328 page 9)\nand a complete answer to this phase transition for multinomial logistic regression is beyond the scope of our paper.  \n\n> I am not sure the convergence of some multinomial regression solver can be used as a rigorous theoretical indicator. The algorithm may converge (or be believed to converge) due to many reasons. Perhaps some better (numerical) criterion can be proposed by, e.g., checking the gradient and/or Hessian of the point of interest.\n\nWe agree that checking the gradients/Hessian at the current iterate is a good criterion. We will clarify this paragraph: Our stance is to leverage (and not reinvent the wheel on) the well-studied convex solvers from optimization which often have their own built-in methods to assess convergence, for instance as returned in the ``success`` and ``message`` attributes of ``scipy.optimize.OptimizeResult`` for L-BFGS which is the default solver in scikit-learn. In practice, a definitive sign that convergence fails and the MLE does not exist is when the solver returns a $\\hat B$ that perfectly separates the data in the sense that the corresponding softmax$(x_i^T\\hat B)$ perfectly interpolates $y_i$ for all $i$. This holds for $K=1$: gradient descent is known to converge to such interpolator (max-margin classifier) on separable data (Soudry, Hoffer, Nacson, Gunasekar and Srebro 2018)., although we are not aware of similar works for $K>1$.\n\n> **Q2:** for the sake of presentation and use, it be helpful to present Theorem 2.2 and the estimation of $\\Omega_{jj}$ in form of an algorithm.\n\n**A2:**\nThanks. \nFor the camera-ready version, we will include another equation in Theorem 2.2 with $\\Omega_{jj}$ replaced by its estimator (2.5).\n\n> **Q3:** Almost nothing is mentioned for the proof of the theoretical results (Theorem 2.1 and 2.2): is the proof technically challenging or contains some ingredients and/or intermediate results that may be of independent interest?\n\n**A3:**\nThe rigorous proof is non-trivial and, due to page limits, we decided to present it in the supplement. The first few pages of the supplement\nhas a diagram explaining the relationship between the different theorems/lemmas, and how the lemmas in the supplement are combined to obtain the theorems stated in the main text.\n\nThe main idea of the method used in the supplement consists of applying Lemma S5.3 to a carefully chosen function $F(\\tilde z)$ taking values in the set of matrices with all singular values equal to 1, which is a new result of independent interest to exhibit asymptotic chi-square distributions in such problems.\nThis poses unique challenges to ensure that the derivatives are well-behaved despite the normalization requiring all singular values to be equal to 1, and we solve some of these challenges specific to the cross-entropy loss in section S6 of the supplement. We expect some of these tools to be of independent interest and useful more broadly (beyond the cross-entropy loss) when one wishes to prove asymptotic chi-square distributions in related problems.\n"}, {"review_id": "poWJlk25B1", "reviewer": "Reviewer_uvEK", "summary": "This paper studies the asymptotic distribution of multinomial logistic regression in the high dimensions. Under the assumption of a Gaussian design (and a few other assumptions), it establishes and characterizes the asymptotic normality of the coefficient of a null feature. This extends previous results for high-dimensional binary logistic regression to $K \\geq 3$ number of classes. This result enables proper significant testing in high-dimensions, for which the classical fixed-$p$ asymptotic fails to control the type-I error. ", "questions": "1. Per my first comment on the \"Weaknesses\", can authors comment on the extent to which $H_0$ can be weakened? \n2. To what extent is the \"universality\" expected to hold for this problem? Does the asymptotic normality fail under a random, heavy-tailed design? What about fixed designs?", "limitations": "I do not foresee any potential negative societal impact of this work. In terms of limitations, I think obtaining results for the non-null features would make this work a much stronger paper. ", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper extends previous results for high-dimensional logistic regression to the multi-class setting. \n2. The paper is well-written and the technical content is presented rigorously. \n3. The simulations results seem to match the asymptotic theory very well. ", "weaknesses": "1. The hypothesis $H_0$ seems stronger than the hypothesis $H_0'$: the population coefficient (in terms of KL projection) for feature $j$ is zero. Supposedly $H_0'$ reduces to $H_0$ when the model is well-specified. In terms of the data generating mechanism given by Assumption 2.1, it seems possible that when $y_i = f(U_i, x_i^{T} B^{\\ast})$ does not hold, a feature that does not satisfy $H_0$ may still have zero as its population \"true\" coefficient. I am wondering whether $H_0$ can be suitably weakened.\n2. The result relies on a random Gaussian design, although it is suggested that the result probably extends to other random designs with light tails. \n3. The paper does not characterize the distribution for non-null features, so the study of parameter inference for high-dimensional multinomial logistic regression is still incomplete. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful feedback to help us improve our paper. Here, we respond to your comments point by point.\n\n> **Q1**: Per my first comment on the \"Weaknesses\", can authors comment on the extent to which can be weakened?\n\n**A1**: The question regarding the possibility to weaken $H_0$ to $H_0'$ is interesting. We do not have an answer right away. The proof is written with a generative model of the form $y_i = f(U_i, x_i^T B^*)$ in mind, and it is not clear at this point that all steps in the proof where this model is used could be generalized to a model-free setting. (Note that many choices for $f$ are allowed which provides significant model flexibility; the stronger assumption here is indeed that the response depends on $x_i$ only through $x_i^TB^*$).\n\n> **Q2**: To what extent is the \"universality\" expected to hold for this problem? Does the asymptotic normality fail under a random, heavy-tailed design? What about fixed designs?\n\n**A2**: Thanks for the great question. We expect our asymptotic normality results (Theorems 2.1 and 2.2) to hold when the covariate $x_i$ is iid sub-Gaussian or has bounded moments of sufficiently large order. This is expected because \n\n1. This is empirically confirmed by the simulation results using the Radermacher distribution and SNP example in Section 3. However, the asymptotical normality fails to hold under random, heavy-tailed design; we conducted additional simulation results using heavy-tailed covariates (e.g., Pareto and Weibull); the corresponding quantile-quantile lines deviate significantly from the 45-degree line. \n2. In the same regime (proportional asymptotics for p/n), several recent works show that universality holds under a set of assumptions that typically includes iid entries with bounded moments of sufficiently large order:\n    - \"Universality of empirical risk minimization\" by Montanari and Saheed (arxiv:2202.08832), cf. section 3.3 there for iid sub-gaussian entries,\n    - \"Universality laws for Gaussian mixtures in generalized linear models\" by Dandi, Stephan, Krzakala, Loureiro, Zdeborová (arxiv:2302.08933) \n    - \"Universality of regularized regression estimators in high dimensions\" by Han and Shen (arxiv:2206.07936) which handles iid entries with $6+\\epsilon$ bounded moments.\n\n\nThese three works provide hope for extending our results to iid entries (or other design distributions treated in these works). However, all these works assume some form of \"delocalization\" of $\\hat B$ in the form of an $L_\\infty$-norm constraint in the minimization problem (see, e.g., equation (29) in Montanari and Saheed). This is fine if $\\hat B$ is already known to satisfy such $L_\\infty$-norm bound, however in logistic regression (even in two classes) such $L_\\infty$-norm bound is not known, and these universality results cannot yet be applied. Given the intricacies involved in proving the delocalization and universality results, we prefer to leave this problem open for future research as these challenging universality questions warrant their own papers. \n\n> I think obtaining results for the non-null features would make this work a much stronger paper.\n\nWe agree that an analysis for non-null features would be very interesting. However, \na major obstacle to obtaining such results is that we operate under the model\n$y_i=f(U_i,x_i^TB^*)$, allowing for a general function $f$ ($f$ must output a\none-hot encoded vector, but is otherwise unrestricted).  In our setting, $B^*$\nis thus not identifiable because $y_i=\\tilde f(U_i, x_i^T\\tilde{B}^*)$, where\n$\\tilde{B}^*=B^* M$ and $\\tilde{f}(u,v^T)=f(u,v^TM^{-1})$ for any invertible matrix\n$M$.  Consequently, turning to confidence intervals for entries of $B^*$ or\nconfidence sets for rows of $B^*$ would require either a model more specific\nthan the current model of the paper (in order to address the identifiability\nissue), or considering that the parameter of interest is the column space of\n$B^*$. Both of these avenues appear interesting for future work. The rationale for\nleaving such extensions for future work is that they would require mathematics\nsignificantly different from the present work (mathematics relying on\nadditional assumptions compared to our model $y_i=f(U_i,x_i^TB^*)$) and we are\nnot aware of known techniques that could advance the analysis for non-null\nfeatures."}, {"review_id": "0MFFS5O1eC", "reviewer": "Reviewer_U9Cr", "summary": "The paper aims to study the statistics of the multinomial logistic MLE on null covariates in a multiclass classification setting in the asymptotic proportional regime, i.e., in the limit in which the size of the dataset $n$ and of the dimension of the covariates $p$ are both sent to infinity but with finite ratio. The work extends previous results by Sur and Candès in which the normality of the MLE on null covariates was proven in the case of binary classification and introduces a proper statistic for a $\\chi^2$ test on the relevance of a feature.", "questions": "Out of curiosity, the authors state that they expect their results to hold for covariates with distribution having \"sufficiently light\" tails. Could they comment on this? Do they expect, for example, sub-gaussianity to be sufficient for their results to hold? \n\nAs a (very) minor remark, Theorem 2.1 and Theorem 2.2 exhibit the quantities $\\mathsf y_i-\\hat{\\mathsf p}_i\\equiv -\\mathsf g_i$, but a different notation is used. The analogy of the two results might appear more evident at first sight using the same notation in both theorems (eg replacing $-\\mathsf y_i+\\hat{\\mathsf p}_i$ with $\\mathsf g_i$ in the first theorem).", "limitations": "The work is of theoretical nature and all limitations are within the well specified hypotheses of the reported theorems.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is clearly and carefully written, the validity conditions of the results are clearly stated and the outcome of the theoretical analysis is supported by robust numerical evidence. The characterization of the statistics of the MLE on null covariates in the case of multiclass classification in the high dimensional limit is timely. The authors also compare the high-d and the classical theory, showing the remarkable difference in the provided predictions. Finally, they also test the validity of their results beyond the Gaussian hypothesis adopted for the derivation of their theorem.", "weaknesses": "A possible theoretical weakness of the paper is the fact that the existence of the MLE is assumed (Assumption 2.4) and not fully characterized. This \"weakness\" is acknowledged by the authors themselves in the final section of the main text.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the insightful and encouraging feedback. Here, we respond to your comments point by point.\n\n> **Q1:** Out of curiosity, the authors state that they expect their results to hold for covariates with distribution having \"sufficiently light\" tails. Could they comment on this? Do they expect, for example, sub-gaussianity to be sufficient for their results to hold?\n\n**A1:**\nThis is a great question. Regarding this question of universality with respect to the distribution of $x_i$, the associated challenges and some recent references, we refer to our answer to question **Q2** from Reviewer uvEK.\n\n> **Q2:** As a (very) minor remark, Theorem 2.1 and Theorem 2.2 exhibit the quantities $\\mathsf{y}_i - \\mathsf{\\hat p}_i = -\\mathsf{g}_i$\n, but a different notation is used. The analogy of the two results might appear more evident at first sight using the same notation in both theorems (eg replacing $-\\mathsf{y}_i + \\mathsf{\\hat p}_i$ with $\\mathsf{g}_i$ in the first theorem).\n\n**A2:**\nThank you very much for your attention to detail. \nWe adopt your suggestion and revised the manuscript by replacing \n$-\\mathsf{y}_i + \\mathsf{\\hat p}_i$ with $\\mathsf{g}_i$ in Theorem 2.1 to make the analogy between the two results more evident."}]

Author rebuttal:

Meta-review: This paper studies the asymptotic distribution of the maximum likelihood estimator for the multinomial logistic regression under the proportional growth regime (where the sample size and the number of unknowns are of the same order). The results developed in this paper generalized the previous theory by Sur et al. to accommodate classification with more than two classes. Asymptotic normality and asymptotic chi-square distributions have been established, which in turn suggest a new hypothesis test to test the significance of any given feature. The paper is well-written, and all reviewers recommend acceptance after the discussion phase. While the paper is mainly of theoretical interest, it would have broader impacts if the authors could demonstrate the efficacy of the proposed solution on more practical data.
Decision: Accept (poster)