The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
90,613 characters
Moment-Based Inference for Regression with Latent Dirichlet Covariates
\maketitle
\begin{abstract}
Topic models are often used as first-stage dimension-reduction tools before
regression, with estimated document-level topic shares treated as observed
covariates. This plug-in workflow creates two inferential difficulties:
valid inference requires a regular first-stage-to-second-stage expansion
that propagates topic-estimation uncertainty, and, at fixed document length,
a document's topic mixture is not consistently recoverable from its own
words even when the population topic matrix is known. Corrected spectral
moment methods for LDA provide a natural starting point: when the total
Dirichlet concentration parameter is known, low-order word moments can be
corrected to yield operators diagonal in the latent topic basis. We extend
this idea to downstream regression. Under a finite latent Dirichlet
allocation model with response residuals orthogonal to the low-order token
moments used for identification, response-weighted word moments admit the
same correction, and the resulting supervised operator identifies the
regression coefficient \(\beta\) directly, without estimating document-level
topic shares. The main theoretical obstacle is that the spectral correction
depends on the unknown total concentration \(\alpha_0\). We show that, for
\(k\ge3\) topics and under a generic finite-probe condition, \(\alpha_0\) is
identifiable by commutativity: at the true value, a family of corrected
word-moment operators commute, whereas away from the truth they generically
do not. This yields a feasible estimator and allows uncertainty in
\(\hat\alpha_0\) to be propagated into inference for \(\beta\). The estimator
is asymptotically linear as the number of documents grows with fixed document
length, with sandwich standard errors based on document-level moment
contributions. Simulations show near-nominal coverage where plug-in
topic-share regressions can undercover, and an application to top economics
journals illustrates contrast inference for latent topic effects.
\end{abstract}
\section{Introduction and Literature Review}\label{sec:intro}
Topic models are widely used as preprocessing tools for empirical
regression analysis. A corpus of high-dimensional text is first mapped
into a low-dimensional representation, typically estimated document-level
topic shares, and these topic-derived variables are then used as
regressors, controls, or predictors in a downstream model. Examples
include newspaper topics in regressions for political violence
\citep{MuellerRauh2018}, news-topic predictors in macroeconomic and
financial forecasting regressions \citep{LarsenThorsrud2019,Thorsrud2020},
topic-based measures of business news in vector autoregressions and
stock-market timing models \citep{BybeeKellyManelaXiu2024}, and topic
measures constructed from central-bank transcripts in econometric studies
of transparency and deliberation \citep{HansenMcMahonPrat2018}. In this
role, latent Dirichlet allocation and related topic models are best
viewed as dimension-reduction steps: they convert sparse word-count data
into a small number of latent coordinates intended for subsequent
statistical use.
This perspective connects topic-model preprocessing to the use of
estimated factors in factor-augmented regressions, where a high-dimensional
panel is compressed into a lower-dimensional latent representation before
a second-stage analysis \citep{StockWatson2002,Bai2003,BaiNg2006}. The
analogy is useful because it highlights the inferential issue. Once an
estimated latent representation is used as a regressor, uncertainty from
the first step must be reconciled with inference in the second step. In
the empirical topic-modeling workflow, this reconciliation is often
absent: researchers estimate topics, construct document-level topic
shares, and then condition on those estimated shares as if they were
ordinary observed covariates.
We consider a canonical finite-LDA version of this problem. Let document
\(i\) contain \(N\) word tokens \(x_{i1},\ldots,x_{iN}\), let \(k\) denote
the number of topics, let \(h_i\) denote its latent topic mixture, and let
\(O\) denote the topic matrix. Conditional on \(h_i\), the word tokens are
generated from the mixture \(Oh_i\). Suppose a scalar response can be written
as
\[
Y_i=\beta^\top h_i+\varepsilon_i,
\]
where the residual is orthogonal to the low-order token moments used for
identification; a sufficient model-level condition is
\[
\mathbb{E}(\varepsilon_i\mid h_i,x_{i1},\ldots,x_{iN})=0 .
\]
The object of inference is the downstream coefficient \(\beta\). We treat
the topic matrix, Dirichlet parameters and regression coefficient as fixed
unknown population parameters; the Dirichlet law for \(h_i\) is a
random-effects distribution for document heterogeneity, not a prior for a
posterior analysis over these unknown parameters. The usual plug-in workflow
estimates topic shares \(\hat h_i\) and then regresses \(Y_i\) on
\(\hat h_i\). This creates two difficulties.
The first is a first-stage regularity problem. Let \(n\) denote the number
of independent documents. For classical downstream inference, the relevant
object is not only the first-stage topic estimate but the entire map from
the observed corpus to the final regression coefficient. A standard
delta-method argument requires this map to be regular in the sense that
\[
\sqrt n(\hat\beta-\beta)
=
n^{-1/2}\sum_{i=1}^n \phi_i+o_p(1),
\]
for an influence function \(\phi_i\) with consistently estimable variance.
Many topic-modeling procedures used in applied work, including Bayesian
and approximate-Bayesian implementations of LDA, produce posterior or
approximate-posterior summaries under chosen prior, hyperparameter and
computational specifications. Such summaries may be useful for exploration,
prediction or posterior inference under the specified model, but they do
not by themselves provide frequentist confidence intervals for the
downstream coefficient \(\beta\). Frequentist validity of the resulting
two-stage map must be established separately. In standard plug-in workflows,
the second-stage regression instead conditions on estimated topic shares as
if they were observed covariates. The resulting standard errors then use
only the second-stage regression variation and omit the first-stage
component of the influence function; unless that omitted component is
asymptotically negligible, the confidence intervals are not valid for the
sampling distribution of the full two-step procedure.
The second difficulty is more fundamental. Even if the population topic
matrix \(O\) were known, a document-specific topic mixture \(h_i\) cannot
be consistently recovered from the words in document \(i\) when document
length is fixed. The number of documents may grow, but the information
about any particular document's latent mixture remains bounded. Hence the
plug-in regressor \(\hat h_i\) contains non-vanishing document-level
error. This is the generated-regressor problem emphasized by
\citet{BattagliaChristensenHansenSacher2024} for regressions using
variables extracted from unstructured data. In the present setting it
arises even under perfect recovery of the topic matrix: increasing the
number of documents improves estimation of population moments, but it
does not make each individual \(h_i\) observed. A regression on \(\hat h_i\) can therefore have a wrong probability
limit or a persistent finite-document bias, not merely standard errors that
condition on a noisy generated regressor.
Integrated latent-variable likelihood procedures or fully Bayesian joint
models could also target \(\beta\) under a more complete specification by
integrating over the latent document mixtures. The contribution here is
complementary: we give a low-order moment route to frequentist confidence
intervals for \(\beta\), based on the supervised moment restrictions stated
below, without specifying a full likelihood for the joint distribution of
responses and words or a prior over the unknown population parameters. The
comparison with plug-in regressions is therefore aimed at the common two-step
topic-share regression workflow rather than at exhausting all possible
latent-variable estimators.
The spectral approach to LDA provides a useful starting point for
resolving both difficulties. \citet{AnandkumarGeHsuKakadeLiu2012} show
that, when the total Dirichlet concentration parameter is known, the
topic matrix can be recovered from low-order cross-token word moments
after applying Dirichlet-specific corrections. These corrections subtract
lower-order terms induced by the common document-level topic mixture,
leaving a population object that is diagonal in latent topic coordinates.
\citet{WangZhu2014} extend this corrected-moment logic to supervised
LDA under a related latent-mixture regression structure, introducing
response-weighted corrected moments that recover the regression weights.
\citet{RenWangZhu2018} develop this supervised spectral approach further,
giving both a two-stage method that recovers LDA parameters before
recovering the regression model and a single-phase method that jointly
recovers topic and regression parameters.
We take this supervised spectral moment structure as the starting point. Our
supervised identity uses the same response-weighted moment logic, but we
formulate it as a conceptually close and operationally different observed-space
operator. The response-weighted correction is formed in observed word space,
and after projection onto the recovered topic basis the representation
\(O^+H^y_{\alpha_0}O\) is diagonal with entries proportional to the components
of \(\beta\). The distinction is therefore not the underlying supervised-LDA
moment model, but the inferential target and estimation problem: we study
frequentist inference for a downstream regression coefficient at fixed document
length, without estimating document-level topic shares.
The remaining roadblock, common to the spectral LDA and supervised spectral
LDA procedures just discussed, is that the corrected moments are indexed by
the total Dirichlet concentration parameter. Write this concentration as
\(\alpha_0\). In these procedures, \(\alpha_0\) is supplied as an input when
forming the corrected moments. In empirical work this quantity is rarely known,
and treating it as a tuning parameter is unsatisfactory for two reasons. First,
the downstream procedure then rests on a choice of \(\alpha_0\) rather than
on an identification argument for \(\alpha_0\) from the observable word
distribution. Practical selection rules are therefore difficult to
interpret as estimators of a population parameter. Second, even if a
particular rule performs well computationally, its uncertainty has no
automatic route into the final confidence interval for \(\beta\). We show
that, when the number of topics satisfies \(k\ge3\) and a generic finite-probe
condition holds, \(\alpha_0\) is itself identified from observable moments:
at the true concentration value, a family of corrected word-moment operators
commute, whereas away from the truth they generically fail to commute. This
gives a feasible estimator of \(\alpha_0\) and allows its uncertainty to be
propagated through the estimator of the downstream coefficient.
We work in a fixed-document-length asymptotic regime: the number of
independent documents grows, while the number of words in each document
may remain bounded. In this regime, the target is inference on the
population coefficient \(\beta\), not recovery of every document-specific
topic mixture. Within this setting, the paper makes three contributions.
First, for known total concentration, we give an observed-space operator
formulation of the supervised spectral moment identity. This identity starts
from the same response-weighted corrected moments as the supervised spectral
LDA literature, but casts them in a form tailored to downstream inference:
the supervised operator is an observed-space matrix whose topic-basis
representation \(O^+H^y_{\alpha_0}O\) is diagonal, and those diagonal entries
are proportional to the components of \(\beta\). This identifies the downstream
coefficient directly from observable word--response moments, without
constructing document-level topic-share regressors.
Second, we identify and estimate the total Dirichlet concentration parameter.
The spectral LDA and supervised spectral LDA procedures discussed above form
their corrected moments conditional on \(\alpha_0\). We show that, when the
number of topics is at least three and a generic probe condition holds,
\(\alpha_0\) is the unique value for which a finite family of corrected
word-moment operators commute. This gives a feasible estimator of \(\alpha_0\)
from the same low-order word moments used for topic recovery.
Third, relative to the supervised spectral recovery results in
\citet{WangZhu2014} and \citet{RenWangZhu2018}, we establish frequentist
inference for the downstream coefficient in a fixed-dimensional regime. The
primitive empirical quantities are distinct-token word moments and
response-weighted word moments, averaged across independent documents. With
fixed vocabulary dimension \(d\), fixed topic dimension \(k\), and fixed document
length \(N\), or after conditioning on a fixed admissible compression dimension,
these moments satisfy an ordinary root-\(n\) central limit theorem over
documents. The concentration-parameter, spectral, and supervised coefficient
maps are smooth under rank and separation conditions. The resulting estimator
of \(\beta\) is asymptotically linear with a feasible sandwich variance
estimator, so the contribution is a first-order sampling expansion and feasible
standard errors rather than only consistency or sample-complexity recovery.
The Monte Carlo experiments are designed to match this inferential logic.
They evaluate the finite-sample performance of the proposed confidence
intervals and compare them with plug-in regressions using estimated
document-topic shares. The proposed intervals have coverage close to
nominal in the designs considered. Plug-in regressions undercover
substantially, including a known-\(O\) plug-in reconstruction procedure in
which the topic matrix is supplied but document-level mixtures are still
reconstructed from finite document counts. This isolates the fixed-document
generated-regressor problem: increasing the number of documents improves
population moment estimation, but it does not make the document-level topic
mixtures observed. The simulations are intended to isolate this mechanism
under favorable rank, separation and probe conditions, not to exhaust all
finite-sample stress tests of the feasible workflow.
The rest of the paper is organized as follows. Section~\ref{sec:model}
introduces the finite-LDA model, the downstream regression target, and
the observable word and word--response moments. Section~\ref{sec:identification}
derives the corrected moment operators, proves population identification
of the topic matrix and the downstream coefficient, and establishes
identification of \(\alpha_0\) through commutativity. Section~\ref{sec:estimation}
defines the proposed moment-based estimators and gives the asymptotic distribution and
sandwich variance estimator for the downstream coefficient. Section~\ref{sec:simulation}
reports Monte Carlo evidence on finite-sample performance and on the
failure of plug-in topic regressions at fixed document length. Section~\ref{sec:application}
gives a real-data illustration using articles from the five general-interest
economics journals. Section~\ref{sec:discussion} discusses limitations and
extensions. Proofs and extensions, including fixed-dimensional observed
controls and split-sample linear compression for high-dimensional preprocessing,
are collected in the Supplement.
\section{Model and target parameter}\label{sec:model}
This section specifies the finite-LDA sampling experiment and the downstream regression target. The point of the formulation is to separate the latent document-level topic mixture, which is useful for defining the model, from the population regression coefficient, which is the object of inference. The corrected spectral moments used for identification are introduced in Section~\ref{sec:identification}.
\subsection{Observed data and latent topic structure}
We observe independent documents
\[
(Y_i,x_{i1},\ldots,x_{iN}),\qquad i=1,\ldots,n,
\]
where \(Y_i\in\mathbb{R}\) is a scalar document-level response. The vocabulary has
\(d\) terms, and each word token is represented by a one-hot vector
\(x_{ij}\in\{e_1,\ldots,e_d\}\subset\mathbb{R}^d\). The main theory is stated for
a common document length \(N\ge3\), fixed as \(n\to\infty\). The condition
\(N\ge3\) is used because the spectral construction relies on third-order
cross-token moments. Variable document lengths \(N_i\ge3\) can be handled by
normalizing the within-document moment averages document by document, as
recorded in Appendix~\ref{app:variable-lengths}; we maintain the same
sampled-document LDA moment structure throughout.
Let \(k\) denote the number of topics, and let
\[
O=[O_{:1},\ldots,O_{:k}]\in\mathbb{R}^{d\times k}
\]
denote the topic matrix. Each column \(O_{:\ell}\) lies in the vocabulary
simplex \(\Delta^{d-1}\) and gives the word distribution for topic
\(\ell\). The matrix \(O\) and the Dirichlet parameter below are treated
as fixed unknown population parameters.
For document \(i\), the latent topic mixture is
\[
h_i\sim\operatorname{Dirichlet}(\alpha),
\qquad
\alpha=(\alpha_1,\ldots,\alpha_k)^\top\in(0,\infty)^k,
\]
with total concentration
\begin{equation}\label{eq:alpha0-def}
\alpha_0:=\sum_{\ell=1}^k\alpha_\ell .
\end{equation}
The usual LDA sampling scheme draws latent topic labels
\[
z_{ij}\mid h_i \sim \operatorname{categorical}(h_i),
\qquad j=1,\ldots,N,
\]
independently across token positions, and then draws words according to
\[
x_{ij}\mid z_{ij}=\ell
\sim \operatorname{categorical}(O_{:\ell}) .
\]
The topic labels \(z_{ij}\) are not observed and will not be used
directly. Integrating them out gives the equivalent conditional word
distribution
\begin{equation}\label{eq:cond-prob}
\Pr(x_{ij}=e_v\mid h_i)=(Oh_i)_v,
\qquad v=1,\ldots,d .
\end{equation}
Thus, conditional on \(h_i\), the word tokens are independent and
identically distributed with conditional mean
\begin{equation}\label{eq:cond-mean}
\mathbb{E}[x_{ij}\mid h_i]=Oh_i .
\end{equation}
Marginally, tokens in the same document are dependent through the common
latent mixture \(h_i\), which is the dependence exploited by the
cross-token moment identities below.
The document-specific vector \(h_i\) is not observed. Because \(N\) is
fixed, the analysis does not require, and does not assume, that \(h_i\)
can be consistently recovered for each document. The inferential target
below is instead a finite-dimensional population parameter identified
from the joint law of \((Y_i,x_{i1},\ldots,x_{iN})\).
\subsection{Downstream response model}
The response is linked to the latent topic mixture by
\begin{equation}\label{eq:supervised-linear-model}
Y_i=\beta^\top h_i+\varepsilon_i,
\end{equation}
for some \(\beta\in\mathbb{R}^k\), with \(\mathbb{E}[Y_i^2]<\infty\). The primitive
restriction used by the supervised moment identities is low-order
response-token orthogonality over the token positions entering the averaged moments:
\begin{equation}\label{eq:response-token-orthogonality}
\mathbb{E}(\varepsilon_i)=0,\qquad
\mathbb{E}(x_{ia}\varepsilon_i)=0\quad\text{for all }a,\qquad
\mathbb{E}(x_{ia}x_{ib}^\top\varepsilon_i)=0\quad\text{for all }a\ne b .
\end{equation}
A sufficient model-level condition is
\begin{equation}\label{eq:strong-response-condition}
\mathbb{E}(\varepsilon_i\mid h_i,x_{i1},\ldots,x_{iN})=0,
\end{equation}
or equivalently
\(\mathbb{E}(Y_i\mid h_i,x_{i1},\ldots,x_{iN})=\beta^\top h_i\). The weaker
condition \(\mathbb{E}(Y_i\mid h_i)=\beta^\top h_i\) alone is not enough for the
response-weighted word moments used below if residual variation is
systematically related to realized token choices beyond \(h_i\).
It should be noted that, since \(\mathbf{1}^\top h_i=1\), a separate intercept is not identified from a
common shift of all topic coefficients. If a model is written as
\(c+\tilde\beta^\top h_i\), then it is observationally equivalent to
\((\tilde\beta+c\mathbf{1})^\top h_i\). We therefore adopt the normalization that absorbs the intercept into
\(\beta\). This is the same normalization issue that arises in a regression
with a full set of category indicators: the common level is a convention, while
contrasts such as \(\beta_a-\beta_b\) are invariant and often the substantively
meaningful objects. Topic labels are arbitrary, so \(\beta\) is always
interpreted in the same ordering as the columns of \(O\). Appendix~\ref{app:observed-controls}
gives the corresponding fixed-dimensional extension with observed controls.
\subsection{Observable low-order moments}
All population moments below refer to a generic document, and we suppress the document index. By exchangeability, any distinct token positions have the same joint distribution, so we use positions \(1,2,3\) for notation. Write
\[
\pi:=\frac{\alpha}{\alpha_0},
\qquad
D:=\operatorname{diag}(\alpha).
\]
The first two Dirichlet moments are
\begin{equation}\label{eq:dir-first}
\mathbb{E}[h]=\pi,
\end{equation}
and
\begin{equation}\label{eq:dir-second-matrix}
\mathbb{E}[hh^\top]
=\frac{D+\alpha\alpha^\top}{\alpha_0(\alpha_0+1)} .
\end{equation}
The corresponding first and second cross-token word moments are
\begin{equation}\label{eq:first-moment}
\mu:=\mathbb{E}[x_1]=O\pi,
\end{equation}
\begin{equation}\label{eq:second-moment}
\begin{aligned}
M_2:=\mathbb{E}[x_1x_2^\top]
&=O\mathbb{E}[hh^\top]O^\top \\
&=\frac{1}{\alpha_0(\alpha_0+1)}ODO^\top
+\frac{\alpha_0}{\alpha_0+1}\mu\mu^\top .
\end{aligned}
\end{equation}
The use of distinct positions in \(M_2\) is important: \(M_2\) is a cross-token moment, not the second moment of a single multinomial draw.
For a contraction direction \(\eta\in\mathbb{R}^d\), define the flattened third cross-token moment
\begin{equation}\label{eq:T-eta-def}
T(\eta):=\mathbb{E}\{x_1x_2^\top\langle x_3,\eta\rangle\}.
\end{equation}
Let \(\gamma:=O^\top\eta\). Conditional independence and \eqref{eq:cond-mean} give
\begin{equation}\label{eq:T-eta-pullout}
T(\eta)=O\,\mathbb{E}\{hh^\top(\gamma^\top h)\}\,O^\top .
\end{equation}
The contracted third Dirichlet moment is
\begin{equation}\label{eq:dir-contracted-third}
\begin{aligned}
\mathbb{E}\{hh^\top(\gamma^\top h)\}
=\frac{1}{\alpha_0(\alpha_0+1)(\alpha_0+2)}
\Big
(& (\alpha^\top\gamma)\alpha\alpha^\top
+\alpha\alpha^\top\operatorname{diag}(\gamma)
+\operatorname{diag}(\gamma)\alpha\alpha^\top \\
& + (\alpha^\top\gamma)D
+2\operatorname{diag}(\alpha\circ\gamma)
\Big).
\end{aligned}
\end{equation}
Thus
\begin{equation}\label{eq:T-eta-explicit}
\begin{aligned}
T(\eta)
=\frac{1}{\alpha_0(\alpha_0+1)(\alpha_0+2)}
O\Big
(& (\alpha^\top\gamma)\alpha\alpha^\top
+\alpha\alpha^\top\operatorname{diag}(\gamma)
+\operatorname{diag}(\gamma)\alpha\alpha^\top \\
& + (\alpha^\top\gamma)D
+2\operatorname{diag}(\alpha\circ\gamma)
\Big)O^\top .
\end{aligned}
\end{equation}
The corrected moments in the next section subtract lower-order terms induced by the common document mixture; at \(\tau=\alpha_0\), this removes the rank-one term in \eqref{eq:second-moment} and all lower-order third-moment terms in \eqref{eq:T-eta-explicit}, leaving only the diagonal factor proportional to \(\operatorname{diag}(\alpha\circ\gamma)\).
Finally, define the supervised observable moments
\begin{equation}\label{eq:supervised-response-moments}
m_y:=\mathbb{E}[Y],
\qquad
v_y:=\mathbb{E}[x_1Y],
\qquad
T^y:=\mathbb{E}[x_1x_2^\top Y].
\end{equation}
Under \eqref{eq:supervised-linear-model} and
\eqref{eq:response-token-orthogonality},
\begin{equation}\label{eq:supervised-moment-factorizations}
m_y=\beta^\top\pi,
\qquad
v_y=O\mathbb{E}[hh^\top]\beta,
\qquad
T^y=O\mathbb{E}\{hh^\top(\beta^\top h)\}O^\top .
\end{equation}
Consequently, whenever \(O\) has full column rank, there exists \(\eta_\beta\in\mathbb{R}^d\) such that
\begin{equation}\label{eq:eta-beta-bridge}
O^\top\eta_\beta=\beta,
\end{equation}
and the supervised moments satisfy the bridge identities
\begin{equation}\label{eq:supervised-bridge-identities}
m_y=\langle \eta_\beta,\mu\rangle,
\qquad
v_y=M_2\eta_\beta,
\qquad
T^y=T(\eta_\beta).
\end{equation}
These identities are the population reason why the downstream coefficient can be recovered from observable word--response moments without first constructing document-level estimates of \(h_i\).
\section{Population identification through corrected operators}\label{sec:identification}
This section gives the population identities on which the estimators are based. The central point is that the document-level mixture \(h_i\) need not be recovered document by document. Instead, low-order observable moments can be corrected so that the resulting operators are diagonal in the latent topic coordinates. One operator recovers the topic directions; its supervised analogue recovers the downstream coefficient vector. A second consequence of the same diagonalization is an identifying restriction for the unknown concentration mass \(\alpha_0\).
Throughout the section, expectations are under the model of Section~\ref{sec:model}. We write \(\tau>0\) for a candidate value of the total concentration and reserve \(\alpha_0\) for the true value. Let
\[
\mu=\mathbb{E}[x_1],\qquad
M_2=\mathbb{E}[x_1x_2^\top],\qquad
T(\eta)=\mathbb{E}\{x_1x_2^\top\langle x_3,\eta\rangle\},
\]
where \(\eta\in\mathbb{R}^d\). Define the \(\tau\)-corrected second moment
\begin{equation}\label{eq:id-Btau}
B_\tau
:=M_2-\frac{\tau}{\tau+1}\mu\mu^\top,
\end{equation}
and the \(\tau\)-corrected contracted third moment
\begin{equation}\label{eq:id-Atau}
\begin{aligned}
A_\tau(\eta)
:= {}&T(\eta)
-\frac{\tau}{\tau+2}
\{M_2\eta\mu^\top+\mu\eta^\top M_2+\langle \eta,\mu\rangle M_2\} \\
&\quad
+\frac{2\tau^2}{(\tau+1)(\tau+2)}
\langle \eta,\mu\rangle\mu\mu^\top .
\end{aligned}
\end{equation}
The associated observed-space operator is
\begin{equation}\label{eq:id-Htau}
H_\tau(\eta):=A_\tau(\eta)B_\tau^+,
\end{equation}
where \(B_\tau^+\) is the Moore--Penrose inverse. The operator is generally not symmetric, so eigenvectors below are right eigenvectors.
\subsection{Corrected LDA moments}\label{subsec:corrected-moments}
The following lemma is the algebraic core of the paper. It states that the correction in \eqref{eq:id-Btau}--\eqref{eq:id-Atau} removes the off-diagonal Dirichlet terms exactly at \(\tau=\alpha_0\). Let
\[
D:=\operatorname{diag}(\alpha),\qquad
C_2:=\alpha_0(\alpha_0+1),\qquad
C_3:=\alpha_0(\alpha_0+1)(\alpha_0+2).
\]
\begin{lemma}[Corrected moment factorization]\label{lem:id-corrected-factorization}
Assume \(O\in\mathbb{R}^{d\times k}\) has full column rank and \(\alpha\in(0,\infty)^k\). For any \(\tau>0\) and \(\eta\in\mathbb{R}^d\), put \(w=O^\top\eta\). Then
\begin{equation}\label{eq:id-B-factor}
B_\tau=O S_B(\tau)O^\top,
\qquad
S_B(\tau)=\frac{1}{C_2}D+
\frac{\alpha_0-\tau}{\alpha_0^2(\alpha_0+1)(\tau+1)}\alpha\alpha^\top .
\end{equation}
Moreover, \(S_B(\tau)\) is positive definite for every \(\tau>0\), with
\begin{equation}\label{eq:id-SB-inv}
S_B(\tau)^{-1}=C_2D^{-1}-(\alpha_0-\tau)\mathbf{1}\mathbf{1}^\top .
\end{equation}
The third-order correction satisfies
\begin{equation}\label{eq:id-A-factor-general}
A_\tau(\eta)=O S_A(\tau;w)O^\top,
\end{equation}
where \(S_A(\tau;w)\) is linear in \(w\). At the true concentration,
\begin{equation}\label{eq:id-true-diagonal-moments}
B_{\alpha_0}=O\frac{D}{C_2}O^\top,
\qquad
A_{\alpha_0}(\eta)=O\frac{2}{C_3}\operatorname{diag}(\alpha\circ w)O^\top .
\end{equation}
Consequently,
\begin{equation}\label{eq:id-H-true-diagonal}
H_{\alpha_0}(\eta)
=O\left\{\frac{2}{\alpha_0+2}\operatorname{diag}(O^\top\eta)\right\}O^+ .
\end{equation}
\end{lemma}
The proof is a direct calculation from the Dirichlet second and third moments and is given in the Supplement. The number of topics is also identified at this point: since \(\mathbb{E}[hh^\top]\) is positive definite and \(O\) has full column rank, \(\operatorname{rank}(M_2)=k\).
\subsection{Topic directions and supervised coefficients when \texorpdfstring{\(\alpha_0\)}{alpha0} is known}\label{subsec:known-alpha-identification}
Equation~\eqref{eq:id-H-true-diagonal} gives the usual spectral recovery of the topic directions, but in observed coordinates rather than whitened coordinates. If \(\eta\) is chosen so that the coordinates of \(O^\top\eta\) are distinct and nonzero, the right eigenvectors associated with the \(k\) nonzero simple eigenvalues of \(H_{\alpha_0}(\eta)\) are exactly the topic vectors.
\begin{theorem}[Topic identification with known concentration]\label{thm:id-known-alpha-topic}
Assume \(O\in\mathbb{R}^{d\times k}\) has full column rank and \(\alpha\in(0,\infty)^k\). Let \(\eta\) be drawn from any distribution on \(\mathbb{R}^d\) that is absolutely continuous with respect to Lebesgue measure. Then, with probability one over \(\eta\), the matrix \(H_{\alpha_0}(\eta)\) has exactly \(k\) nonzero simple eigenvalues,
\begin{equation}\label{eq:id-topic-eigenvalues}
\lambda_j(\eta)=\frac{2}{\alpha_0+2}(O^\top\eta)_j,
\qquad j=1,\ldots,k,
\end{equation}
and the associated right eigenvectors are \(O_{:1},\ldots,O_{:k}\). Thus \(O\) is identified up to column permutation and scaling; the simplex normalization \(\mathbf{1}^\top O_{:j}=1\) fixes the scaling.
\end{theorem}
The downstream coefficient is identified by a parallel supervised correction. Using the response moments from Section~\ref{sec:model}, define for a candidate \(\tau\)
\begin{equation}\label{eq:id-Aytau}
\begin{aligned}
A_\tau^y
:={}&T^y
-\frac{\tau}{\tau+2}
\{v_y\mu^\top+\mu v_y^\top+m_y M_2\} \\
&\quad
+\frac{2\tau^2}{(\tau+1)(\tau+2)}m_y\mu\mu^\top,
\end{aligned}
\end{equation}
and
\begin{equation}\label{eq:id-Hytau}
H_\tau^y:=A_\tau^y B_\tau^+ .
\end{equation}
The key identity is obtained at \(\tau=\alpha_0\).
\begin{theorem}[Direct identification of the downstream coefficient]\label{thm:id-beta-supervised}
Assume the model of Section~\ref{sec:model}, including the response-token
orthogonality conditions in \eqref{eq:response-token-orthogonality}. If
\(O\) has full column rank, then
\begin{equation}\label{eq:id-Ay-diagonal}
A_{\alpha_0}^y
=O\frac{2}{C_3}\operatorname{diag}(\alpha\circ\beta)O^\top
\end{equation}
and
\begin{equation}\label{eq:id-Hy-diagonal}
H_{\alpha_0}^y
=O\left\{\frac{2}{\alpha_0+2}\operatorname{diag}(\beta)\right\}O^+ .
\end{equation}
Consequently, once the topic columns have been labeled by the unsupervised spectral step,
\begin{equation}\label{eq:id-beta-projection}
\beta_j
=\frac{\alpha_0+2}{2}\{O^+H_{\alpha_0}^yO\}_{jj},
\qquad j=1,\ldots,k .
\end{equation}
This identification does not require the entries of \(\beta\) to be distinct and does not require document-level estimates of \(h_i\).
\end{theorem}
The theorem follows from the bridge identities
\(m_y=\langle\eta_\beta,\mu\rangle\), \(v_y=M_2\eta_\beta\), and
\(T^y=T(\eta_\beta)\) for any \(\eta_\beta\) satisfying
\(O^\top\eta_\beta=\beta\). These identities are implied by the low-order
orthogonality conditions in \eqref{eq:response-token-orthogonality}; the
stronger conditional mean restriction \eqref{eq:strong-response-condition}
is only one sufficient primitive condition. Thus
\(A_{\alpha_0}^y=A_{\alpha_0}(\eta_\beta)\), and
Lemma~\ref{lem:id-corrected-factorization} applies. Equation~\eqref{eq:id-Hy-diagonal}
shows that the supervised operator is diagonal in the topic basis. If the
entries of \(\beta\) are distinct, this relation may be read as an eigenvalue
statement. For inference, however, we use the projection formula
\eqref{eq:id-beta-projection}, which avoids the complications caused by
repeated or zero coefficient values once the topic basis has been fixed.
\subsection{Identification of the concentration mass}\label{subsec:alpha0-identification}
The preceding identities still require the total concentration \(\alpha_0\)
as an input. This is a central practical limitation of corrected spectral LDA
moments: the correction that diagonalizes the word moments depends on a
quantity that is rarely known in applications. We now show that the same
operator representation identifies \(\alpha_0\) from observable word moments.
The identifying restriction is commutativity.
At the true value, \eqref{eq:id-H-true-diagonal} implies that every corrected
operator is diagonal in the same topic basis. Hence any two such operators
commute: for every pair \(\eta_1,\eta_2\),
\begin{equation}\label{eq:id-commute-true}
[H_{\alpha_0}(\eta_1),H_{\alpha_0}(\eta_2)]=0,
\end{equation}
where \([R,S]=RS-SR\). We show that this property is generically unique to
the true concentration: away from \(\alpha_0\), mean-orthogonal contractions
produce nonzero commutators. Thus commutativity supplies an identifying
restriction for the concentration mass. Let
\begin{equation}\label{eq:id-projection-def}
P_\mu:=I_d-\frac{\mu\mu^\top}{\|\mu\|^2}.
\end{equation}
If \(\eta\in\mu^\perp\), then \(w=O^\top\eta\in\alpha^\perp\), because
\[
\alpha^\top O^\top\eta=\alpha_0\mu^\top\eta=0.
\]
\begin{theorem}[Identification of \texorpdfstring{\(\alpha_0\)}{alpha0} by commutativity]\label{thm:id-alpha-commutativity}
Assume \(k\ge3\), \(O\in\mathbb{R}^{d\times k}\) has full column rank, and \(\alpha\in(0,\infty)^k\). For \(\tau>0\), write
\[
H_\tau(\eta)=O T_\tau(O^\top\eta)O^+,
\qquad
T_\tau(w):=S_A(\tau;w)S_B(\tau)^{-1}.
\]
Then the following statements hold.
\begin{enumerate}[label=\textnormal{(\roman*)},leftmargin=2.5em]
\item At \(\tau=\alpha_0\), the family \(\{H_{\alpha_0}(\eta):\eta\in\mathbb{R}^d\}\) is pairwise commuting.
\item If \(\tau\ne\alpha_0\) and \(w_1,w_2\in\alpha^\perp\) are non-collinear, then
\begin{equation}\label{eq:id-latent-commutator}
[T_\tau(w_1),T_\tau(w_2)]
=c(\tau,\alpha_0)
\{(\alpha\circ w_1)w_2^\top-(\alpha\circ w_2)w_1^\top\},
\end{equation}
where
\begin{equation}\label{eq:id-c-commutator}
c(\tau,\alpha_0)
=\frac{4(\alpha_0-\tau)(\alpha_0\tau+\alpha_0+\tau)}
{\alpha_0(\alpha_0+1)(\alpha_0+2)^2(\tau+2)^2} .
\end{equation}
In particular, \([T_\tau(w_1),T_\tau(w_2)]\ne0\), and hence any \(\eta_1,\eta_2\in\mu^\perp\) satisfying \(O^\top\eta_j=w_j\) give
\[
[H_\tau(\eta_1),H_\tau(\eta_2)]\ne0 .
\]
\end{enumerate}
Consequently, \(\alpha_0\) is the unique value of \(\tau>0\) for which the full family \(\{H_\tau(\eta):\eta\in\mu^\perp\}\) is pairwise commuting.
\end{theorem}
The restriction \(k\ge3\) is used only here. It ensures that \(\alpha^\perp\) has dimension at least two, so non-collinear mean-orthogonal contractions exist. The present commutator argument is therefore formulated for topic models with at least three topics.
For estimation, one does not need to check commutativity over all \(\eta\in\mu^\perp\). Let \(s_1,\ldots,s_L\in\mathbb{R}^d\), define
\begin{equation}\label{eq:id-finite-probes}
v_\ell=P_\mu s_\ell,
\qquad \ell=1,\ldots,L,
\end{equation}
and let \(\mathcal I\subset\{(\ell,q):1\le \ell<q\le L\}\) be non-empty. Consider the population finite-probe criterion
\begin{equation}\label{eq:id-pop-criterion}
Q_{\mathcal I}(\tau)
:=\sum_{(\ell,q)\in\mathcal I}
\big\|[H_\tau(v_\ell),H_\tau(v_q)]\big\|_F^2 .
\end{equation}
If there is at least one pair \((\ell,q)\in\mathcal I\) for which \(O^\top v_\ell\) and \(O^\top v_q\) are non-collinear, then
\begin{equation}\label{eq:id-Q-identifies-alpha}
Q_{\mathcal I}(\tau)=0
\quad\Longleftrightarrow\quad
\tau=\alpha_0 .
\end{equation}
If two raw probes \(s_1,s_2\) are drawn independently from distributions that are absolutely continuous on \(\mathbb{R}^d\), this non-collinearity condition holds with probability one. The sample estimator in Section~\ref{sec:estimation} is the plug-in version of \eqref{eq:id-pop-criterion}, using empirical moments and empirical mean-orthogonal projections.
The commutativity result also has a purely unsupervised implication. Combined
with the corrected-moment eigenvector identity in Theorem~\ref{thm:id-known-alpha-topic},
it implies that the finite-dimensional LDA model is identified from low-order
word moments without supplying the total concentration as an input. We record
this consequence explicitly.
\begin{corollary}[Unsupervised finite-LDA identification]\label{cor:unsupervised-lda-identification}
Under the assumptions of \Cref{thm:id-alpha-commutativity}, the word-side parameters of the finite LDA model are identified from the joint law of three distinct tokens, without supplying \(\alpha_0\) as an input. In particular, \(k=\operatorname{rank}(M_2)\), \(\alpha_0\) is identified by the commutativity condition in \Cref{thm:id-alpha-commutativity}, the topic matrix \(O\) is identified up to column permutation by \Cref{thm:id-known-alpha-topic} using any contraction direction with distinct nonzero latent coordinates, and the simplex normalizations \(\mathbf{1}^\top O_{:j}=1\) fix the column scales. The Dirichlet parameter is then identified by
\[
\pi=O^+\mu,\qquad \alpha=\alpha_0\pi .
\]
\end{corollary}
\begin{proof}
The rank statement follows from \(M_2=O\mathbb{E}(hh^\top)O^\top\), since \(\mathbb{E}(hh^\top)\) is positive definite and \(O\) has full column rank. The commutator condition identifies \(\alpha_0\). Given \(\alpha_0\), \Cref{thm:id-known-alpha-topic} identifies the columns of \(O\) up to permutation and scale, and the simplex normalizations fix the scale. Finally, \(\mu=O\alpha/\alpha_0\), so \(\pi=\alpha/\alpha_0=O^+\mu\), and \(\alpha=\alpha_0\pi\). Equivalently, the same correction can be viewed as the supervised bridge with pseudo-response \(Y_\eta=\langle x_3,\eta\rangle\), using \(x_1,x_2\) for the word moments. Conditional independence of distinct tokens gives
\[
\mathbb{E}(Y_\eta)=\langle\eta,\mu\rangle,\qquad
\mathbb{E}(x_1Y_\eta)=M_2\eta,\qquad
\mathbb{E}(x_1x_2^\top Y_\eta)=T(\eta),
\]
so the corrected pseudo-response operator is exactly \(H_{\alpha_0}(\eta)\).
\end{proof}
\section{Estimation and inference}\label{sec:estimation}
This section turns the population identities of Section~\ref{sec:identification} into a feasible estimator. The construction has three steps. First, estimate the observable cross-token moments by within-document averages over distinct token positions. Second, estimate the concentration mass \(\alpha_0\) by minimizing a finite set of empirical commutators. Third, recover the topic matrix from one unsupervised spectral operator and recover the downstream coefficient from the supervised operator. The limit theory treats these steps as one smooth map of empirical moments, after rank selection and topic order have been determined.
\subsection{Empirical cross-token moments}\label{subsec:empirical-moments}
For document \(i\), define the within-document empirical mean
\begin{equation}\label{eq:est-mu-i}
\hat\mu_i:=\frac1N\sum_{a=1}^N x_{ia}.
\end{equation}
The second cross-token moment is estimated by averaging over ordered distinct pairs,
\begin{equation}\label{eq:est-M2-i}
\hat M_{2,i}
:=\frac{1}{N(N-1)}\sum_{a\ne b}x_{ia}x_{ib}^\top .
\end{equation}
Because ordered pairs include both \((a,b)\) and \((b,a)\), \(\hat M_{2,i}\)
is symmetric. Equivalently, if \(c_i=\sum_{a=1}^N x_{ia}\) is the word-count vector, then
\begin{equation}\label{eq:est-M2-counts}
\hat M_{2,i}
=\frac{c_ic_i^\top-\operatorname{diag}(c_i)}{N(N-1)} .
\end{equation}
For a contraction direction \(v\in\mathbb{R}^d\), the third cross-token moment is estimated by
\begin{equation}\label{eq:est-T-i}
\hat T_i(v)
:=\frac{1}{N(N-1)(N-2)}
\sum_{\substack{a,b,c=1\\ a,b,c\;\mathrm{distinct}}}^N
x_{ia}x_{ib}^\top\langle x_{ic},v\rangle .
\end{equation}
The corresponding sample averages are
\begin{equation}\label{eq:est-raw-averages}
\hat\mu:=\frac1n\sum_{i=1}^n\hat\mu_i,
\qquad
\hat M_2:=\frac1n\sum_{i=1}^n\hat M_{2,i},
\qquad
\hat T(v):=\frac1n\sum_{i=1}^n\hat T_i(v).
\end{equation}
The hats in this section refer to the sample size \(n\), which is suppressed when no ambiguity can arise. The estimators in \eqref{eq:est-raw-averages} are unbiased for \(\mu\), \(M_2\), and \(T(v)\), respectively. The use of distinct token positions removes the single-token multinomial contribution that would otherwise appear on the diagonal of the second moment.
The supervised moments are estimated in the same way. Define
\begin{equation}\label{eq:est-supervised-moments}
\hat m_y:=\frac1n\sum_{i=1}^nY_i,
\qquad
\hat v_y:=\frac1n\sum_{i=1}^nY_i\hat\mu_i,
\qquad
\hat T^y:=\frac1n\sum_{i=1}^nY_i\hat M_{2,i}.
\end{equation}
Thus \(\hat T^y\) is the response-weighted analogue of the second cross-token moment. The full third-order tensor need not be stored in computation; only the finite number of contractions \(\hat T(v)\) used by the estimator are required.
\emph{Large-vocabulary preprocessing.} The formulae above are written in the
original vocabulary dimension. In large vocabularies, one can instead compress
tokens linearly before forming the corrected moments. Appendix~\ref{app:pca-compression}
shows that if a fixed-dimensional compression matrix does not lose a topic
direction, in the sense that \(\operatorname{rank}(R^\top O)=k\), then the corrected moment
identities and the coefficient \(\beta\) are unchanged. When observed controls
are included after compression, the control adjustment also requires the latent
topic-share scale. Appendix~\ref{app:pca-compression} gives a compressed
scale-recovery step based on the compressed first and corrected second moments;
this step is used before constructing the control cross-moment \(M_{qh}\). The
split-sample PCA version estimates an \(m\)-dimensional compression, with fixed
\(m\ge k\), from the uncentered cross-token second moment on an independent
split, holds it fixed, and runs the moment estimator on the compressed tokens in
the second split. This is intended as a compatibility result for fixed \(k\) and
fixed \(m\), not a rigorous growing-vocabulary or growing-rank asymptotic theory.
\subsection{Plug-in corrected operators}\label{subsec:plugin-operators}
For each candidate \(\tau>0\), define the empirical corrected moments
\begin{equation}\label{eq:est-Btau}
\hat B_\tau
:=\hat M_2-\frac{\tau}{\tau+1}\hat\mu\hat\mu^\top,
\end{equation}
and
\begin{equation}\label{eq:est-Atau}
\begin{aligned}
\hat A_\tau(v)
:={}&\hat T(v)
-\frac{\tau}{\tau+2}
\{\hat M_2v\hat\mu^\top+\hat\mu v^\top\hat M_2
+\langle v,\hat\mu\rangle\hat M_2\} \\
&\quad
+\frac{2\tau^2}{(\tau+1)(\tau+2)}
\langle v,\hat\mu\rangle\hat\mu\hat\mu^\top .
\end{aligned}
\end{equation}
Similarly, the supervised corrected moment is
\begin{equation}\label{eq:est-Aytau}
\begin{aligned}
\hat A_\tau^y
:={}&\hat T^y
-\frac{\tau}{\tau+2}
\{\hat v_y\hat\mu^\top+\hat\mu\hat v_y^\top+\hat m_y\hat M_2\} \\
&\quad
+\frac{2\tau^2}{(\tau+1)(\tau+2)}
\hat m_y\hat\mu\hat\mu^\top .
\end{aligned}
\end{equation}
These are the sample analogues of \eqref{eq:id-Btau}, \eqref{eq:id-Atau}, and \eqref{eq:id-Aytau}.
The corrected second moment has rank \(k\) in population, so the inverse in the empirical operator is taken on an estimated rank-\(k\) subspace. Let
\[
\hat B_\tau=\sum_{j=1}^d\hat\lambda_j(\tau)\hat u_j(\tau)\hat u_j(\tau)^\top,
\qquad
\hat\lambda_1(\tau)\ge\cdots\ge\hat\lambda_d(\tau),
\]
be a spectral decomposition. For a candidate rank \(m\), set
\begin{equation}\label{eq:est-truncated-inverse}
\hat B_{\tau,m}^+
:=\sum_{j=1}^m\hat\lambda_j(\tau)^{-1}
\hat u_j(\tau)\hat u_j(\tau)^\top,
\end{equation}
whenever \(\hat\lambda_m(\tau)>0\). If \(m=0\), set \(\hat B_{\tau,0}^+=0\). Values of \(\tau\) for which \(\hat\lambda_m(\tau)\le0\) are excluded from the minimization below. The empirical unsupervised and supervised operators are
\begin{equation}\label{eq:est-Htau}
\hat H_{\tau,m}(v):=\hat A_\tau(v)\hat B_{\tau,m}^+,
\qquad
\hat H_{\tau,m}^y:=\hat A_\tau^y\hat B_{\tau,m}^+ .
\end{equation}
If \(k\) is not specified in advance, estimate it from the rank of \(M_2\). Let \(\tilde\lambda_1\ge\cdots\ge\tilde\lambda_d\) be the eigenvalues of \(\hat M_2\), and let \(a_n\) satisfy
\begin{equation}\label{eq:est-threshold-rate}
a_n\to0,
\qquad
\sqrt n\,a_n\to\infty .
\end{equation}
Define
\begin{equation}\label{eq:est-khat}
\hat k:=\sum_{j=1}^d\mathbf 1\{\tilde\lambda_j>a_n\}.
\end{equation}
When the number of topics is chosen by design or by an external model-selection step, the same formulas apply with \(\hat k\) replaced by that chosen value.
In finite samples, implementing the rank rule in \eqref{eq:est-khat} still requires a numerical choice of threshold, and nearby candidate ranks may have similar empirical second-moment spectra. When a split-sample implementation is used, as in the application below, rank choice can be made on the first split before the response-weighted moment estimation step. We use the commutator identity as a practical response-free rule of thumb: for each candidate \(k\), compute the commutator-based estimate \(\hat\alpha_0(k)\) across several admissible probe draws on the first split, and prefer ranks for which this estimate is interior and stable. This is a finite-sample rule of thumb for rank selection, not an additional rank-consistency theorem. The motivation is that, under the maintained finite-LDA model and the correct rank, admissible finite-probe commutator criteria have the same population minimizer \(\alpha_0\). Thus a working rank for which \(\hat\alpha_0(k)\) moves substantially across reasonable probe draws is less attractive than one for which the concentration estimate is interior and stable. Subsequent estimation and standard errors are computed on the independent second split, conditional on the selected rank, compression, and realized probes.
This split-sample design also allows the commutator calculation to be used a second time as a diagnostic step. Since the working rank is chosen on the first split, the commutator profile can be inspected again on the second split as a holdout diagnostic for the selected finite-dimensional specification. If the same sample were used both to choose the rank and to assess commutator stability, this second calculation would mainly repeat the selection criterion. With an independent second split, boundary behavior, flatness, multimodality, or probe sensitivity provides a separate warning about the stability of the selected spectral diagonalization. We discuss these post-selection diagnostics in Section~\ref{subsec:estimate-alpha0}.
\subsection{Estimating the concentration mass}\label{subsec:estimate-alpha0}
Let \(P(u)=I_d-uu^\top/\|u\|^2\) for \(u\ne0\). Choose raw probe vectors \(s_1,\ldots,s_L\in\mathbb{R}^d\), with \(L\ge2\), and define their empirical mean-orthogonal projections
\begin{equation}\label{eq:est-probes}
\hat v_\ell:=P(\hat\mu)s_\ell,
\qquad
\ell=1,\ldots,L .
\end{equation}
Let \(\mathcal I\subset\{(\ell,q):1\le\ell<q\le L\}\) be a non-empty set of probe pairs. For fixed rank \(m\), define the empirical commutator criterion
\begin{equation}\label{eq:est-Q}
\hat Q_{\mathcal I}^{(m)}(\tau)
:=\sum_{(\ell,q)\in\mathcal I}
\left\|
[\hat H_{\tau,m}(\hat v_\ell),\hat H_{\tau,m}(\hat v_q)]
\right\|_F^2,
\end{equation}
where \([R,S]=RS-SR\). If \(\hat\lambda_m(\tau)\le0\), set \(\hat Q_{\mathcal I}^{(m)}(\tau)=+\infty\). The estimator of the total concentration is
\begin{equation}\label{eq:est-alpha0}
\hat\alpha_0
:=\min\arg\min_{\tau\in\Theta}
\hat Q_{\mathcal I}^{(\hat k)}(\tau),
\end{equation}
where \(\Theta=[\underline\tau,\overline\tau]\subset(0,\infty)\) is compact and contains the true \(\alpha_0\) in its interior. For the true rank \(m=k\), the finite-probe identification result in Section~\ref{sec:identification} implies that the population version of \eqref{eq:est-Q} is uniquely minimized at \(\alpha_0\), provided at least one tested pair of projected probes is non-collinear in the latent coordinates. Thus the sample criterion has the correct local target on the event \(\hat k=k\). No uniqueness claim is made here for a misspecified truncation rank.
The probability-one genericity statements are identification results, not
finite-sample guarantees. In applied work one should inspect the profile
of \(\hat Q_{\mathcal I}^{(m)}(\tau)\), whether \(\hat\alpha_0\) lies near
the boundary of \(\Theta\), local curvature around the minimizer, eigenvalue
gaps for the ordering operator, and sensitivity to the realized probe and
ordering directions. A flat, multimodal or boundary-attaining commutator
criterion, or estimates that move materially across admissible probes, should
be treated as a substantive warning. It may indicate weak finite-sample
identification, an unstable rank choice, ill-conditioned topic separation,
an inadmissible compression, or misspecification of the maintained LDA moment
structure.
\subsection{Topic and downstream coefficient estimators}\label{subsec:topic-beta-estimators}
The concentration estimator is then plugged into a single spectral operator to recover the topic columns. Choose an ordering direction \(r\in\mathbb{R}^d\), independent of the data if random, and set
\begin{equation}\label{eq:est-ordering-direction}
\hat\eta:=P(\hat\mu)r .
\end{equation}
Define
\begin{equation}\label{eq:est-topic-operator}
\hat H^{o}:=\hat H_{\hat\alpha_0,\hat k}(\hat\eta).
\end{equation}
On the high-probability event that \(\hat H^{o}\) has \(\hat k\) eigenvalues
separated from the zero cluster, and that these selected eigenvalues are real,
simple and mutually separated, let
\(\hat o_1(r),\ldots,\hat o_{\hat k}(r)\) be the associated right eigenvectors.
The selected nonzero eigenpairs are labeled by decreasing real eigenvalue within this selected cluster, and the columns are normalized by
\begin{equation}\label{eq:est-topic-normalization}
\mathbf{1}^\top\hat o_j(r)=1 .
\end{equation}
Off this event, the estimator may be defined by any fixed deterministic convention; this off-event convention is asymptotically irrelevant.
Set
\begin{equation}\label{eq:est-Ohat}
\hat O(r):=[\hat o_1(r),\ldots,\hat o_{\hat k}(r)].
\end{equation}
The role of \(r\) is only to select a stable ordering of the topic columns. At the population value, the eigenvectors are the columns of \(O\), ordered by the entries of \(O^\top P(\mu)r\).
Finally, form the supervised operator at the estimated concentration,
\begin{equation}\label{eq:est-Hy-final}
\hat H^y:=\hat H^y_{\hat\alpha_0,\hat k}
=\hat A_{\hat\alpha_0}^y\hat B_{\hat\alpha_0,\hat k}^+ .
\end{equation}
The downstream coefficient estimator is
\begin{equation}\label{eq:est-beta}
\hat\beta(r)
:=\frac{\hat\alpha_0+2}{2}
\operatorname{diag}\{\hat O(r)^+\hat H^y\hat O(r)\},
\end{equation}
where \(\operatorname{diag}(\cdot)\) extracts the diagonal as a vector. This is the direct sample analogue of \eqref{eq:id-beta-projection}. It does not estimate document-level topic shares and does not run a second-stage regression on generated regressors.
\subsection{Asymptotic linearity and standard errors}\label{subsec:asymptotics-inference}
The asymptotic argument is fixed-dimensional: \(d\), \(k\), and \(N\) are fixed while \(n\to\infty\). Let \(\hat{\mathcal T}_i\) denote the ordered-distinct triple tensor
\[
\hat{\mathcal T}_i
:=\frac{1}{N(N-1)(N-2)}
\sum_{\substack{a,b,c=1\\ a,b,c\;\mathrm{distinct}}}^N
x_{ia}\otimes x_{ib}\otimes x_{ic},
\]
so that \(\hat T_i(v)=\hat{\mathcal T}_i(I,I,v)\). Collect the document-level moments needed by the estimator in
\begin{equation}\label{eq:est-Zi}
Z_i:=\{\hat\mu_i,\hat M_{2,i},\hat{\mathcal T}_i,
Y_i,Y_i\hat\mu_i,Y_i\hat M_{2,i}\}.
\end{equation}
Let \(Z_0=\mathbb{E} Z_i\) and \(\bar Z=n^{-1}\sum_i Z_i\). Since the word components are bounded, the condition \(\mathbb{E} Y^2<\infty\) is enough for
\begin{equation}\label{eq:est-Z-clt}
\sqrt n(\bar Z-Z_0)\rightsquigarrow N(0,\Omega),
\qquad
\Omega:=\operatorname{var}(Z_i).
\end{equation}
On the event that the rank estimate is correct and the selected nonzero-cluster ordering eigenvalues are real, simple, mutually separated, and separated from the zero cluster, all estimators above are functions of \(\bar Z\). The commutator minimization is differentiable at the population value because the stacked commutator has a nonzero derivative with respect to \(\tau\) at \(\alpha_0\). The eigenvector and pseudoinverse maps are differentiable on the corresponding fixed-rank, separated-cluster stratum.
The following theorem records the resulting first-order expansion. If the probes or ordering direction are random, the statement is conditional on their realized values.
This inferential expansion is one of the key distinctions from the supervised spectral recovery literature:
the goal is not only consistent recovery of supervised-LDA parameters, but a first-order sampling
expansion and feasible variance estimator for the downstream coefficient.
\begin{theorem}[Feasible asymptotic inference]\label{thm:est-main-asymptotics}
Assume the model of Section~\ref{sec:model} with fixed \(d\), fixed \(k\ge3\), fixed \(N\ge3\), \(O\) full column rank, \(\alpha\in(0,\infty)^k\), and \(\mathbb{E} Y^2<\infty\). Let \(\Theta\) be compact with \(\alpha_0\in\mathrm{int}(\Theta)\), and let the rank threshold satisfy \eqref{eq:est-threshold-rate}. Suppose the finite-probe criterion has at least one pair \((\ell,q)\in\mathcal I\) for which \(O^\top P(\mu)s_\ell\) and \(O^\top P(\mu)s_q\) are non-collinear. Suppose also that the ordering direction \(r\) satisfies that the coordinates of \(O^\top P(\mu)r\) are distinct and nonzero.
Then
\begin{equation}\label{eq:est-consistency}
\Pr(\hat k=k)\to1,
\qquad
\hat\alpha_0\xrightarrow{p}\alpha_0,
\qquad
\hat O(r)\xrightarrow{p}O_r,
\qquad
\hat\beta(r)\xrightarrow{p}\beta_r,
\end{equation}
where \(O_r\) is the topic matrix in the ordering induced by \(r\), and \(\beta_r\) is \(\beta\) in the same ordering. Moreover, there exists a mean-zero influence function \(\phi_i(r)\) with finite second moment such that
\begin{equation}\label{eq:est-joint-al}
\sqrt n
\begin{pmatrix}
\hat\alpha_0-\alpha_0\\
\operatorname{vec}\{\hat O(r)-O_r\}\\
\hat\beta(r)-\beta_r
\end{pmatrix}
=\frac1{\sqrt n}\sum_{i=1}^n\phi_i(r)+o_p(1).
\end{equation}
Consequently,
\begin{equation}\label{eq:est-joint-clt}
\sqrt n
\begin{pmatrix}
\hat\alpha_0-\alpha_0\\
\operatorname{vec}\{\hat O(r)-O_r\}\\
\hat\beta(r)-\beta_r
\end{pmatrix}
\rightsquigarrow N\{0,V(r)\},
\qquad
V(r):=\operatorname{var}\{\phi_i(r)\}.
\end{equation}
If the probes and \(r\) are drawn independently from distributions that are absolutely continuous with respect to Lebesgue measure, the finite-probe non-collinearity condition and the distinct nonzero ordering coordinates hold with probability one over those draws. The positive spectral separation of the corrected second-moment cluster follows from the model and compactness of \(\Theta\); the ordering-operator separation from the zero cluster follows from the nonzero ordering coordinates.
\end{theorem}
\begin{corollary}[Conditional inference after split compression and controls]\label{cor:compressed-control-inference}
Under the conditions of Theorem~\ref{thm:est-main-asymptotics} and
Assumption~\ref{ass:split-compressed-controls}, consider the split-sample
compressed estimator with observed controls described in
Appendices~\ref{app:observed-controls}--\ref{app:pca-compression}. Let
\(\mathcal F_1\) denote the sigma-field generated by the first-split compression,
selected working rank, finite probes, and ordering direction, and let
\(n_2=|I_2|\to\infty\). Then, conditionally on \(\mathcal F_1\),
\begin{equation}\label{eq:est-compressed-control-corollary}
\sqrt{n_2}\{\hat\beta^c_{\hat R_1}(r)-\beta_{\hat R_1,r}\}
=\frac{1}{\sqrt{n_2}}\sum_{i\in I_2}\phi^c_{\hat R_1,i}(r)
+o_p(1\mid \mathcal F_1),
\end{equation}
where
\[
\phi^c_{\hat R_1,i}(r)
=D^c_{\beta,r}\{Z^c_{\hat R_1,i}-Z^c_{\hat R_1,0}\},
\qquad
\mathbb{E}\{\phi^c_{\hat R_1,i}(r)\mid\mathcal F_1\}=0 .
\]
The empirical sandwich estimator obtained from the enlarged compressed control
moment vector \(Z^c_{\hat R_1,i}\) is consistent for the conditional covariance
of the right side of \eqref{eq:est-compressed-control-corollary}. The first-split
compression and first-split rank choice add no separate first-order term after
conditioning. If the conditional covariance converges in probability to a
nonrandom limit, the same normal approximation holds unconditionally.
\end{corollary}
The proof of the theorem and corollary is given in the Supplement. Its main steps are standard but useful to separate. The document-level moments satisfy the central limit theorem \eqref{eq:est-Z-clt}. The minimizer \(\hat\alpha_0\) has the usual one-dimensional minimum-distance expansion obtained by differentiating the stacked commutator equations at \(\alpha_0\). The topic estimator then follows from first-order perturbation theory for simple right eigenvectors. Finally, \eqref{eq:est-beta} is a smooth map of \((\hat\alpha_0,\hat B_{\hat\alpha_0,\hat k}^+,\hat A_{\hat\alpha_0}^y,\hat O)\); the first-order perturbation of \(\hat O\) contributes only off-diagonal terms inside \(\hat O^+\hat H^y\hat O\), so the diagonal map recovers the coefficient influence function without requiring distinct entries of \(\beta\).
For implementation, let \(G_{\beta,r}\) denote the population map from \(Z_0\) to \(\beta_r\) induced by \eqref{eq:est-alpha0}--\eqref{eq:est-beta}, with the realized probes and ordering direction held fixed. Let
\begin{equation}\label{eq:est-D-beta}
D_{\beta,r}:=\left.\frac{\partial G_{\beta,r}(z)}{\partial z^\top}\right|_{z=Z_0} .
\end{equation}
Then the lower block of \(V(r)\) is
\begin{equation}\label{eq:est-Vbeta}
V_\beta(r)=D_{\beta,r}\Omega D_{\beta,r}^\top .
\end{equation}
All reported standard errors use the influence-function implementation of this derivative. Write the finite-probe commutator map as \(g(\tau,z)\), and let \(\alpha(z)\) be the local minimum-distance solution. At an interior solution,
\begin{equation}\label{eq:est-D-alpha-impl}
D_\alpha
:=\left.\frac{\partial \alpha(z)}{\partial z^\top}\right|_{z=Z_0}
=-\{G_\tau^\top G_\tau\}^{-1}G_\tau^\top G_z,
\qquad
G_\tau=\left.\frac{\partial g(\tau,z)}{\partial \tau}\right|_{(\alpha_0,Z_0)},\quad
G_z=\left.\frac{\partial g(\tau,z)}{\partial z^\top}\right|_{(\alpha_0,Z_0)} .
\end{equation}
For fixed \(\alpha_0\), the derivative of the coefficient map is obtained by differentiating the corrected moment operators, the truncated inverse of the corrected second moment, the normalized right eigenvectors used to form \(\hat O(r)\), and the final diagonal map in \eqref{eq:est-beta}. The full derivative has the chain-rule form
\begin{equation}\label{eq:est-D-beta-chain}
D_{\beta,r}=D_{\beta,r}^{\rm fix}+D_{\beta,\alpha}D_\alpha,
\end{equation}
where \(D_{\beta,r}^{\rm fix}\) is the derivative with \(\alpha_0\) held fixed and \(D_{\beta,\alpha}\) is the derivative of the coefficient map with respect to \(\alpha_0\). The empirical influence values are therefore
\begin{equation}\label{eq:est-emp-if-beta}
\hat\phi_{\beta,i}(r)=\hat D_{\beta,r}\{Z_i-\bar Z\} .
\end{equation}
Equivalently, with
\begin{equation}\label{eq:est-sandwich}
\hat\Omega:=\frac1n\sum_{i=1}^n\{Z_i-\bar Z\}\{Z_i-\bar Z\}^\top,
\qquad
\hat V_\beta(r)=\hat D_{\beta,r}\hat\Omega\hat D_{\beta,r}^\top
=\frac1n\sum_{i=1}^n\hat\phi_{\beta,i}(r)\hat\phi_{\beta,i}(r)^\top,
\end{equation}
where \(\hat D_{\beta,r}\) is the sample analogue of \eqref{eq:est-D-beta-chain}. Equations~\eqref{eq:est-emp-if-beta}--\eqref{eq:est-sandwich} describe the no-control moment vector. If observed controls are included, the same formulas are applied after replacing \(Z_i\), \(\bar Z\), \(D_{\beta,r}\), and \(\Omega\) by the enlarged control moment vector \(Z_i^c\), its mean, the derivative of the control-adjusted map, and \(\operatorname{var}(Z_i^c)\). If the data are first compressed, \(Z_i^c\) is formed from the compressed tokens and includes the compressed scale-recovery moments used by the control adjustment. A pointwise \((1-a)\) confidence interval for the \(j\)-th ordered coefficient is
\begin{equation}\label{eq:est-ci}
\hat\beta_j(r)
\pm
z_{1-a/2}\left\{\frac{\hat V_{\beta,jj}(r)}{n}\right\}^{1/2}.
\end{equation}
The same covariance estimate is used for Wald tests and linear contrasts of the ordered coefficient vector. For pointwise confidence intervals we assume \(V_{\beta,jj}(r)>0\). For Wald tests of \(C\beta=c\), we assume \(CV_\beta(r)C^\top\) is nonsingular.
\section{Monte Carlo experiments}\label{sec:simulation}
This section evaluates the finite-sample behavior of the proposed estimator in the
setting for which the theory is designed: downstream regression with a fixed number
of words per document. The experiments are organized around two questions. First,
does the direct moment-based estimator deliver calibrated inference for the
coefficient vector \(\beta\)? Second, what goes wrong if the analyst instead
constructs document-level topic estimates and treats them as observed regressors?
The number of topics is supplied to all methods, so the exercise isolates the
inferential problem rather than the separate rank-selection step.
To make the fixed-document plug-in issue explicit, consider any deterministic reconstruction rule
\(\tilde h_i=g_N(c_i;O)\) based on the document counts and the known population topic matrix.
If a no-intercept population least-squares regression of \(Y_i\) on \(\tilde h_i\) is nonsingular and
\(\mathbb{E}(\tilde h_i\varepsilon_i)=0\), its probability limit is
\begin{equation}\label{eq:plugin-plim}
\beta_g
=\{\mathbb{E}(\tilde h_i\tilde h_i^\top)\}^{-1}\mathbb{E}(\tilde h_i h_i^\top)\beta .
\end{equation}
Thus the plug-in coefficient equals \(\beta\) only under additional restrictions, such as
\(\mathbb{E}(\tilde h_i h_i^\top)=\mathbb{E}(\tilde h_i\tilde h_i^\top)\). The correctly specified posterior mean
\(\tilde h_i=\mathbb{E}(h_i\mid c_i)\) is a special case satisfying this identity by iterated expectations, but a generic finite-\(N\) reconstruction rule, including the constrained least-squares rule below, need not satisfy it. The simulations therefore use the known-\(O\) plug-in regression as an analytic and numerical benchmark for the persistent generated-regressor problem at fixed document length.
\subsection{Design and methods}
The vocabulary size is \(d=100\), the number of topics is \(k=10\), and each
document contains \(N=100\) words. The topic matrix \(O\) is generated once and
then held fixed across Monte Carlo replications. Its columns are independently
generated from a Dirichlet distribution with concentration parameter \(0.3\),
subject to screening restrictions that remove nearly duplicate topics and nearly
anchor-word designs. For the realized matrix, the minimum pairwise Hellinger
distance between topic columns is \(0.588\), the maximum row dominance is
\(0.879\), the minimum Hellinger distance from any row-normalized word profile to
a simplex vertex is \(0.250\), and there are no exact zeros.
For document \(i\),
\[
h_i\sim\operatorname{Dirichlet}(\alpha),\qquad
c_i\mid h_i\sim\operatorname{Multinomial}(N,Oh_i),
\]
and the response is
\[
Y_i=\beta^\top h_i+\varepsilon_i,
\qquad
\varepsilon_i\sim N(0,\sigma^2),
\qquad
\beta=(1.0,0.9,\ldots,0.1)^\top .
\]
The noise variance is calibrated to give population \(R^2=0.35\). We consider a
symmetric Dirichlet design, \(\alpha_j=0.5\) for \(j=1,\ldots,10\), so that
\(\alpha_0=5\), and an asymmetric design, \(\alpha_j=2j/55\), so that
\(\alpha_0=2\). The main downstream comparison uses \(n\in\{3000,5000,15000\}\); each design cell uses \(1000\) Monte Carlo replications. The concentration-parameter table below also reports the original \(n=30000\) design point.
We compare four procedures. The first is an infeasible no-intercept latent-oracle regression
that uses the true \(h_i\)'s. The second is the proposed direct spectral estimator,
which estimates \(\alpha_0\), recovers the topic matrix, and estimates \(\beta\)
from the supervised corrected moments. Its confidence intervals use the
influence-function implementation of the first-order variance formula developed
above: the derivative of \(\hat\alpha_0\) is obtained from the implicit
minimum-distance expansion of the commutator criterion, and the derivative of
\(\hat\beta\) is the fixed-\(\alpha_0\) spectral derivative plus the chain-rule
contribution from estimating \(\alpha_0\). The third
procedure is a known-\(O\) plug-in regression: it uses the true topic matrix \(O\), but
estimates document-level topic mixtures by
\[
\tilde h_i
=\arg\min_{h\in\Delta^{k-1}}
\left\|\frac{c_i}{N}-Oh\right\|_2^2,
\]
and then applies no-intercept ordinary least squares treating \(\tilde h_i\) as
observed. The no-intercept specification is natural because
\(\mathbf 1^\top h_i=1\). The fourth procedure is the analogous spectral plug-in regression, replacing \(O\) by
the unsupervised spectral estimate \(\hat O\). Thus the final two procedures
mimic the modular workflow used in many empirical applications. Topic labels are
aligned to the true labels by Hungarian matching only when computing simulation
errors and coverage.
For each method, Table~\ref{tab:sim-beta} reports the root mean squared
error
\[
\operatorname{RMSE}_\beta
=
\left\{
\frac{1}{1000}\sum_{r=1}^{1000}
\frac{1}{k}\|\hat\beta^{(r)}-\beta\|_2^2
\right\}^{1/2},
\]
the average coordinate-wise coverage of nominal \(95\%\) confidence intervals, the
minimum coordinate-wise coverage, average interval length, and average Hellinger
distance between matched topic columns.
\begin{table}[!htbp]
\centering
\small
\setlength{\tabcolsep}{4.5pt}
\begin{threeparttable}
\caption{Downstream inference at fixed document length}
\label{tab:sim-beta}
\begin{tabular}{@{}llrrrrr@{}}
\toprule
\(n\) & Method
& \(\operatorname{RMSE}_\beta\)
& Cov.
& Min cov.
& Length
& Topic \(H\) \\
\midrule
\multicolumn{7}{@{}l}{\textit{Panel A: symmetric Dirichlet, \(\alpha_0=5\)}} \\
3,000 & Latent oracle & 0.022 & 0.949 & 0.935 & 0.085 & -- \\
3,000 & Direct spectral & 0.073 & 0.978 & 0.957 & 0.452 & 0.205 \\
3,000 & Plug-in, true \(O\) & 0.046 & 0.607 & 0.126 & 0.086 & 0.000 \\
3,000 & Plug-in, \(\hat O\) & 0.068 & 0.529 & 0.270 & 0.088 & 0.205 \\
\addlinespace
5,000 & Latent oracle & 0.017 & 0.954 & 0.943 & 0.066 & -- \\
5,000 & Direct spectral & 0.058 & 0.969 & 0.944 & 0.331 & 0.178 \\
5,000 & Plug-in, true \(O\) & 0.044 & 0.492 & 0.016 & 0.067 & 0.000 \\
5,000 & Plug-in, \(\hat O\) & 0.063 & 0.461 & 0.226 & 0.068 & 0.178 \\
\addlinespace
15,000 & Latent oracle & 0.010 & 0.951 & 0.934 & 0.038 & -- \\
15,000 & Direct spectral & 0.027 & 0.962 & 0.937 & 0.113 & 0.125 \\
15,000 & Plug-in, true \(O\) & 0.042 & 0.272 & 0.000 & 0.039 & 0.000 \\
15,000 & Plug-in, \(\hat O\) & 0.052 & 0.288 & 0.028 & 0.038 & 0.125 \\
\addlinespace[1mm]
\multicolumn{7}{@{}l}{\textit{Panel B: asymmetric Dirichlet, \(\alpha_0=2\)}} \\
3,000 & Latent oracle & 0.024 & 0.950 & 0.940 & 0.087 & -- \\
3,000 & Direct spectral & 0.072 & 0.966 & 0.948 & 0.447 & 0.168 \\
3,000 & Plug-in, true \(O\) & 0.045 & 0.818 & 0.422 & 0.091 & 0.000 \\
3,000 & Plug-in, \(\hat O\) & 0.079 & 0.647 & 0.328 & 0.091 & 0.168 \\
\addlinespace
5,000 & Latent oracle & 0.018 & 0.954 & 0.938 & 0.067 & -- \\
5,000 & Direct spectral & 0.056 & 0.967 & 0.947 & 0.300 & 0.145 \\
5,000 & Plug-in, true \(O\) & 0.042 & 0.746 & 0.168 & 0.070 & 0.000 \\
5,000 & Plug-in, \(\hat O\) & 0.069 & 0.613 & 0.252 & 0.070 & 0.145 \\
\addlinespace
15,000 & Latent oracle & 0.011 & 0.952 & 0.936 & 0.039 & -- \\
15,000 & Direct spectral & 0.029 & 0.962 & 0.943 & 0.141 & 0.103 \\
15,000 & Plug-in, true \(O\) & 0.039 & 0.506 & 0.000 & 0.041 & 0.000 \\
15,000 & Plug-in, \(\hat O\) & 0.051 & 0.460 & 0.059 & 0.040 & 0.103 \\
\bottomrule
\end{tabular}
\begin{tablenotes}
\footnotesize
\item[] \textit{Notes:} Cov. is the average coordinate-wise coverage of nominal
\(95\%\) confidence intervals; Min cov. is the minimum over the ten coordinates.
Length is the average coordinate-wise interval length. Topic \(H\) is the
average Hellinger distance between matched topic columns; it is not applicable for
the latent-oracle regression and equals zero for methods using the true topic
matrix.
\end{tablenotes}
\end{threeparttable}
\end{table}
The contrast between the direct estimator and the plug-in regressions is the main
finding. In the symmetric design, the direct spectral estimator has average
coverage \(0.978\), \(0.969\), and \(0.962\) as \(n\) increases from \(3000\) to
\(15000\). Over the same sample sizes, the spectral plug-in regression has
coverage \(0.529\), \(0.461\), and \(0.288\). The known-\(O\) plug-in regression,
which is given the true topic matrix, also undercovers severely: its coverage
falls from \(0.607\) to \(0.272\). Thus the failure is not primarily a topic
recovery error. It is the fixed-document generated-regressor error induced by
estimating \(h_i\) from only \(N=100\) words and then treating \(\tilde h_i\) as
observed.
The asymmetric design gives the same qualitative conclusion. The direct spectral
estimator remains close to nominal coverage, whereas the true-topic plug-in
coverage falls from \(0.818\) to \(0.506\), and the spectral plug-in coverage
falls from \(0.647\) to \(0.460\). The plug-in intervals have lengths comparable
to the latent-oracle intervals, so they shrink at the usual \(n^{-1/2}\) rate even
though their centers retain finite-document bias. With \(1000\) replications, the
Monte Carlo standard error of a single coordinate's nominal \(95\%\) coverage
estimate is about \(\{0.95(0.05)/1000\}^{1/2}=0.0069\), far smaller than the
reported plug-in undercoverage.
As a further diagnostic, we reran the spectral plug-in regression while giving it
the true value of \(\alpha_0\). Across the six design cells in
Table~\ref{tab:sim-beta}, this changed average coverage by at most
\(0.039\). The plug-in failure is therefore not explained by uncertainty in the
concentration parameter; it comes from using finite-document topic estimates as
regressors.
\subsection{Estimation of the concentration parameter}
Table~\ref{tab:sim-alpha} reports the behavior of the commutator estimator of
\(\alpha_0\). Bias and RMSE decrease steadily with \(n\) in both Dirichlet
designs. At the two largest sample sizes, the empirical standard deviation and
the average estimated standard error are close, supporting the first-order
linearization used in the downstream variance formula. Together with
Table~\ref{tab:sim-beta}, this indicates that the proposed method learns
\(\alpha_0\) from word moments and delivers calibrated downstream inference
without plug-in document-topic regressors.
\begin{table}[!htbp]
\centering
\scriptsize
\setlength{\tabcolsep}{5pt}
\begin{threeparttable}
\caption{Estimation of \(\alpha_0\)}
\label{tab:sim-alpha}
\begin{tabular}{@{}lrrrrrrr@{}}
\toprule
\(n\)
& Mean
& Bias
& RMSE
& Cov.
& Emp. sd
& Mean se
& se/sd \\
\midrule
\multicolumn{8}{@{}l}{\textit{Panel A: symmetric Dirichlet, \(\alpha_0=5\)}} \\
3,000 & 5.984 & 0.984 & 1.364 & 0.961 & 0.944 & 1.109 & 1.174 \\
5,000 & 5.634 & 0.634 & 1.000 & 0.949 & 0.774 & 0.812 & 1.050 \\
15,000 & 5.209 & 0.209 & 0.472 & 0.951 & 0.423 & 0.429 & 1.012 \\
30,000 & 5.112 & 0.112 & 0.322 & 0.941 & 0.302 & 0.297 & 0.984 \\
\addlinespace
\multicolumn{8}{@{}l}{\textit{Panel B: asymmetric Dirichlet, \(\alpha_0=2\)}} \\
3,000 & 2.284 & 0.284 & 0.413 & 0.928 & 0.300 & 0.324 & 1.080 \\
5,000 & 2.167 & 0.167 & 0.273 & 0.940 & 0.216 & 0.236 & 1.095 \\
15,000 & 2.054 & 0.054 & 0.134 & 0.951 & 0.123 & 0.127 & 1.034 \\
30,000 & 2.028 & 0.028 & 0.091 & 0.948 & 0.087 & 0.088 & 1.018 \\
\bottomrule
\end{tabular}
\begin{tablenotes}
\footnotesize
\item[] \textit{Notes:} Results are for the proposed direct spectral estimator.
Cov. is the empirical coverage of the nominal \(95\%\) confidence interval for
\(\alpha_0\). The last column is the ratio of the average estimated standard
error to the Monte Carlo standard deviation of \(\hat\alpha_0\).
\end{tablenotes}
\end{threeparttable}
\end{table}
The simulation designs deliberately hold fixed the topic dimension, topic
separation, document length, concentration search interval and probe
construction. They are therefore best read as tests of the fixed-document
generated-regressor mechanism and of the first-order variance formula in
well-separated designs. They do not replace application-specific diagnostics
for the commutator criterion, ordering eigengaps, rank selection or
compression stability.
\section{Application: abstract topics and citations in economics}\label{sec:application}
We illustrate the method on articles published in the five general-interest economics journals
\emph{American Economic Review}, \emph{Econometrica}, \emph{Journal of Political Economy},
\emph{Quarterly Journal of Economics}, and \emph{Review of Economic Studies}. Article metadata and
citation counts are taken from OpenAlex \citep{PriemPiwowarOrr2022,OpenAlexData2026}. The response is
\[
Y_i=\log\{1+\text{citation count}_i\},
\]
and the document consists of the article title and abstract. The analysis is descriptive: a topic
coefficient is interpreted as an association between latent abstract content and citation impact,
not as a causal effect of writing on that topic. The exercise should also be read under the maintained
low-order response-token orthogonality restrictions, which cannot be verified from the observed text alone.
The maintained response equation in this application is
\[
Y_i = \beta^\top h_i + \delta^\top q_i + \varepsilon_i,
\]
where \(h_i\) is the latent topic mixture and \(q_i\) contains publication-year indicators. Since
\(\mathbf{1}^\top h_i=1\), a common shift in all topic coefficients is observationally equivalent to an
intercept shift. Thus the absolute level of an individual coefficient \(\beta_j\), and a test of
\(\beta_j=0\), depends on the normalization implicit in the intercept. The invariant objects are
topic contrasts,
\[
\beta_a-\beta_b,
\]
which compare the conditional mean response when latent topic mass is shifted from topic \(b\) to
topic \(a\). A reallocation of 10 percentage points of topic mass from \(b\) to \(a\) corresponds to
\(0.1(\beta_a-\beta_b)\) log-points in \(Y_i\), holding the year controls fixed. We therefore report
the topic coefficient estimates to orient the reader, but interpret the empirical results through
joint and pairwise contrast tests.
The raw corpus contains records published between 2000 and 2018. We exclude records whose titles or
metadata indicate comments, replies, corrections, editorial material, or very short documents. We
lowercase tokens, remove standard English and scholarly stopwords, and keep terms whose document
frequency is at least 40. The main specification does not stem words. The resulting analysis sample
has 5409 documents and 1771 terms. Because the vocabulary dimension is too large for direct dense
third-order moment calculations, we use the split-sample compression described in
Appendix~\ref{app:pca-compression}. A 10\% split of documents is used to estimate a 100-dimensional
orthonormal compression from the uncentered cross-token second moment, and the remaining 90\% is
used for estimation and inference. All reported confidence intervals and standard errors are computed
on the second split, conditional on the first-split compression. This follows the split-compression
construction in Appendix~\ref{app:pca-compression}, under which the first-stage compression is held
fixed for second-split inference. The empirical moments use document-specific length normalizations
as in Appendix~\ref{app:variable-lengths}.
Following the split-sample finite-rank diagnostic described in Section~\ref{subsec:plugin-operators}, the working topic dimension is selected using only first-split word moments, without using citation outcomes. For each candidate \(k\in\{10,\ldots,20\}\), we estimate \(\alpha_0\) from the commutator criterion across 20 independent sets of probe directions, using five mean-orthogonal commutator partner directions for each run. Candidates with boundary or failed concentration estimates are discarded. Among the remaining candidates, we choose the value of \(k\) with the smallest standard deviation of \(\log \hat\alpha_0\). This first-split commutator-stability diagnostic selects \(k=15\).
After fixing this rank and the first-split compression, we run the estimator on the independent second split and inspect the second-split commutator profile, boundary behavior, compressed scale recovery, and control-adjustment conditioning as diagnostics for the selected specification. In a second-split probe diagnostic using five admissible probe draws, \(\hat\alpha_0\) ranges from 5.09 to 7.59, and none of the corresponding commutator profiles has a boundary minimum. The second-split estimator uses year fixed effects, with year 2000 as the omitted category. The coefficient estimates and standard errors are computed in compressed coordinates after the moment-based topic-scale recovery used in the observed-control adjustment; lifted word-space directions are used only to label topics. The selected specification has
\[
\hat\alpha_0=5.09,
\qquad
n_{\rm est}=4868,
\qquad
\operatorname{cond}(\hat\Gamma_q)=6.56,
\]
where \(\hat\Gamma_q\) is the matrix in the observed-control adjustment. The minimum absolute recovered scale factor is 1.78, and the value of \(\operatorname{cond}(\hat\Gamma_q)\) is moderate in this application, so the year-control adjustment is not close to singular by this diagnostic.
Table~\ref{tab:top5-main} reports the topic labels and centered coefficient estimates in the selected ordering. The labels are assigned from the leading positive words of the lifted word-space directions. The centered estimates subtract the average topic coefficient and are shown only to orient the reader; the invariant inferential objects are the joint equality test and pairwise contrasts reported below.
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{Topic labels and centered coefficients for top-five economics articles}
\label{tab:top5-main}
\scriptsize
\setlength{\tabcolsep}{3.6pt}
\begin{tabular}{r>{\raggedright\arraybackslash}p{0.24\linewidth}>{\raggedright\arraybackslash}p{0.38\linewidth}rr}
\toprule
Topic & Label & Leading words & Centered est. & S.e. \\
\midrule
1 & Risk / insurance / asset preferences
& risk, insurance, financial, market, sharing, utility
& 0.317 & 0.364 \\
2 & Capital / growth / inequality
& capital, growth, tax, model, market, income
& 0.035 & 0.289 \\
3 & Tax / income / health / public finance
& tax, income, taxes, consumption, health, reform
& -0.016 & 0.327 \\
4 & Macro / monetary / shocks
& model, rate, policy, monetary, shocks, inflation
& 1.096 & 0.216 \\
5 & Markets / policy / welfare
& market, policy, price, markets, optimal, welfare
& -0.124 & 0.244 \\
6 & Monetary policy / inflation
& policy, monetary, rate, inflation, growth, shocks
& 0.442 & 0.273 \\
7 & Labor search / wages
& market, model, price, labor, workers, search
& 0.138 & 0.384 \\
8 & Consumer choice / prices / demand
& price, prices, choice, preferences, utility, consumers
& 0.287 & 0.422 \\
9 & Auctions / prices / products
& price, prices, auction, auctions, firms, demand
& -0.381 & 0.189 \\
10 & Trade / international
& trade, countries, free, country, international, growth
& 0.563 & 2.203 \\
11 & Game theory / auctions / equilibrium
& equilibrium, games, trade, game, players, auctions
& -1.334 & 2.381 \\
12 & Mechanism design / principal--agent
& optimal, agent, principal, mechanisms, contracts, risk
& -0.750 & 0.904 \\
13 & Firms / trade / productivity
& firms, firm, model, trade, workers, productivity
& 0.962 & 0.274 \\
14 & Private information / mechanism design
& information, private, agents, agent, models, equilibrium
& -0.460 & 0.302 \\
15 & Econometrics / identification / asymptotics
& models, estimator, estimators, estimation, asymptotic, identification
& -0.773 & 0.230 \\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\scriptsize
\item Notes: The response is \(\log\{1+\text{citation count}\}\). Documents are title and abstract text.
The sample contains top-five economics articles from 2000--2018. The estimator uses a 10\%/90\%
split-sample PCA compression with \(m=100\), first-split commutator-stability selection of \(k\)
using word moments only, and year fixed effects. The coefficient estimates and standard errors are computed in compressed coordinates after moment-based topic-scale recovery. The displayed estimates are \(\hat\beta_j-\bar{\hat\beta}\), with standard errors computed from the full analytic sandwich covariance. Centering removes the common level that absorbs the intercept normalization because topic shares sum to one; pairwise contrasts are unchanged by this centering.
\end{tablenotes}
\end{threeparttable}
\end{table}
The topics align with recognizable fields and methods in economics. The recovered topic labeled
``Econometrics / identification / asymptotics'', for example, has leading words
\emph{models}, \emph{estimator}, \emph{estimators}, \emph{estimation}, \emph{asymptotic},
and \emph{identification}. Other topics correspond to macroeconomic shocks, monetary policy,
labor search, public finance, firms and productivity, auctions, mechanism design, and private
information. This is a qualitative check that the spectral directions are not merely recovering common
article boilerplate.
The natural omnibus null in this setting is not \(\beta=0\), but equality of all topic coefficients,
\[
H_0:\beta_1=\cdots=\beta_{15}.
\]
This null says that, after controlling for publication year, latent topic composition has no association
with citation impact. Using the full influence-function covariance matrix, the Wald statistic for this
14-dimensional restriction is 56.26, with \(p=5.25\times10^{-7}\). Thus, within the maintained
latent-topic model, abstract topic composition is associated with citation impact.
Table~\ref{tab:top5-contrasts} reports the pairwise contrasts that remain significant after Holm
adjustment over all 105 pairwise comparisons. The largest contrast is between the macro/monetary
shocks topic and the econometrics/identification topic. A ten-percentage-point shift of latent topic
mass from the latter to the former corresponds to 0.187 log-points in the citation outcome, holding
publication year fixed. The firms/trade/productivity topic also has a larger citation association than
the econometrics/identification topic and the auctions/prices/products topic. These are descriptive
contrasts in latent topic coordinates, not causal comparisons across research fields.
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{Holm-significant topic contrasts}
\label{tab:top5-contrasts}
\scriptsize
\setlength{\tabcolsep}{3.0pt}
\begin{tabular}{>{\raggedright\arraybackslash}p{0.255\linewidth}>{\raggedright\arraybackslash}p{0.285\linewidth}rrrr}
\toprule
Higher topic & Lower topic & Difference & S.e. & \(p\) & Holm \(p\) \\
\midrule
Macro / monetary / shocks
& Econometrics / identification / asymptotics
& 1.868 & 0.323 & \(7.6\times10^{-9}\) & \(8.0\times10^{-7}\) \\
Macro / monetary / shocks
& Auctions / prices / products
& 1.477 & 0.292 & \(4.2\times10^{-7}\) & \(4.36\times10^{-5}\) \\
Firms / trade / productivity
& Econometrics / identification / asymptotics
& 1.734 & 0.364 & \(2.0\times10^{-6}\) & \(2.01\times10^{-4}\) \\
Macro / monetary / shocks
& Private information / mechanism design
& 1.556 & 0.394 & \(8.0\times10^{-5}\) & 0.0082 \\
Macro / monetary / shocks
& Markets / policy / welfare
& 1.220 & 0.310 & \(8.2\times10^{-5}\) & 0.0083 \\
Firms / trade / productivity
& Auctions / prices / products
& 1.343 & 0.356 & \(1.6\times10^{-4}\) & 0.016 \\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\scriptsize
\item Notes: Contrasts are differences between the corresponding topic coefficients in
Table~\ref{tab:top5-main}. Holm \(p\)-values adjust over all \(15\times14/2=105\) pairwise topic
comparisons. Standard errors use the full analytic sandwich covariance, including the contribution
from estimating \(\alpha_0\) and the compressed scale-recovery step used for the observed-control adjustment.
\end{tablenotes}
\end{threeparttable}
\end{table}
The direction of these contrasts is consistent with existing evidence on citation heterogeneity across
fields in economics. \citet{CardDellaVigna2013} study the same five general-interest journals and
document substantial differences in citation patterns across fields, including relatively high citation
performance for more recent Development and International Economics papers and lower relative
citation performance for more recent Econometrics and Theory papers. \citet{AnauatiGalianiGalvez2016}
classify top-five articles into applied, applied theory, econometric methods, and theory, and find
that citation life cycles differ markedly across these categories: applied and applied-theory papers
have more favorable citation profiles, while theory and econometric-method papers have lower
typical citation profiles, with econometric methods displaying a more heterogeneous upper tail. Our
topic-level estimates should not be read as a replication of those field classifications, but the broad
pattern is similar: topic content is associated with citation outcomes, and the more applied macro and
firm/trade/productivity directions have higher citation associations than some econometric and
theory-oriented directions.
The application is intended as a numerical illustration of the moment-based inferential workflow. It
uses a single first-split selection rule, fixed before estimating the second-split coefficients, and all
reported uncertainty is computed on the independent second split, conditional on the first-split
compression and selected working dimension. Appendix~\ref{app:application-robustness} reports second-split diagnostics, adjacent-rank checks, and a preprocessing robustness check. The main text is kept to one specification to avoid turning the illustration into a model-selection exercise.
\section{Discussion}\label{sec:discussion}
This paper develops direct inference for regression with latent Dirichlet
topic covariates. The main distinction from the usual plug-in workflow is
that document-level topic shares are not estimated and then treated as
observed regressors. At fixed document length, those document-level
quantities remain noisy even when the topic matrix is known. The proposed
response-weighted corrected moments instead identify the downstream
coefficient \(\beta\) directly from observable word--response moments under
explicit low-order response-token orthogonality conditions. The same
corrected moment structure also yields an estimating equation for the
unknown total concentration parameter \(\alpha_0\), through the
commutativity of corrected word-moment operators.
The results are model-based and should be interpreted within the stated
asymptotic regime. The formal theory keeps the vocabulary dimension, number
of topics and document length fixed while the number of documents increases.
The commutativity identification of \(\alpha_0\) is formulated for
\(k\ge3\) and relies on rank, separation and generic probe conditions.
These genericity conditions should not be read as a substitute for numerical
diagnostics in applied use. Boundary solutions, flat commutator profiles,
small ordering eigengaps, unstable rank choices or substantial sensitivity to
admissible probes are substantive warnings. They may reflect weak finite-sample
identification, ill-conditioned topic separation, an unstable compression, or
misspecification of the maintained LDA moment structure. The split-compression
results in the Supplement identify \(\beta\) conditional on an admissible
projection; compressed directions are not, by themselves, simplex-normalized
topics.
Integrated latent-variable likelihood or Bayesian joint-model procedures
could also be used to target \(\beta\) under fuller specifications by
integrating over the latent document mixtures. The contribution here is a
different tradeoff: low-order moment-based frequentist inference for
\(\beta\), without document-level topic-share regression and without a full
likelihood or prior specification for all unknown population parameters.
Several extensions remain. A theory allowing the vocabulary dimension or
the compression dimension to grow with the number of documents would better
match large-scale text applications. It would also be useful to study
robustness under departures from the Dirichlet topic-mixture assumption,
including non-Dirichlet latent mixtures, correlated topics and misspecified
topic dimension. Finally, the same moment logic may be useful beyond topic
models, in other latent simplex models where the inferential target is a
downstream regression coefficient rather than recovery of each unit's latent
coordinate.
\newpage