Latent Factor Point Processes for Patient Representation in Electronic Health Records

TL;DR

This paper introduces a latent factor point process model for EHR data that captures complex temporal structures and underlying disease processes, improving patient representation for classification and clustering tasks.

Contribution

It proposes a novel low-rank point process model and Fourier-Eigen embedding for efficient, interpretable patient representations from high-dimensional EHR data.

Findings

Effective in capturing subgroup-specific temporal patterns

Improves classification and clustering of patient data

Demonstrated on Alzheimer's disease cohort

Abstract

Electronic health records (EHR) contain valuable longitudinal patient-level information, yet most statistical methods reduce the irregular timing of EHR codes into simple counts, thereby discarding rich temporal structure. Existing temporal models often impose restrictive parametric assumptions or are tailored to code level rather than patient-level tasks. We propose the latent factor point process model, which represents code occurrences as a high-dimensional point process whose conditional intensity is driven by a low dimensional latent Poisson process. This low-rank structure reflects the clinical reality that thousands of codes are governed by a small number of underlying disease processes, while enabling statistically efficient estimation in high dimensions. Building on this model, we introduce the Fourier-Eigen embedding, a patient representation constructed from the spectral density matrix of the observed process. We establish theoretical guarantees showing that these embeddings efficiently capture subgroup-specific temporal patterns for downstream classification and clustering. Simulations and an application to an Alzheimer's disease EHR cohort demonstrate the practical advantages of our approach in uncovering clinically meaningful heterogeneity.

Tables

Table 1. Table 1 : Average AUC and adjusted Rand index over 100 replications for varying specifications of the transfer kernel function β ​ ( t ) \beta(t) . The score of the best performing method for each β ​ ( t ) \beta(t) is rendered in bold. The final column gives the average score across all four specifications of β ​ ( t ) \beta(t) .

Metric	Embedding					Average
AUC	Fourier-Eigen	$0.8868$	$0.7406$	$0.6645$	0.7256	$0.7544$
	PMI	0.8763	0.6654	0.6252	$0.8126$	0.7449
	Counts	0.8866	0.5547	0.5650	0.6836	0.6725
Adj. Rand	Fourier-Eigen	$0.9606$	0.8595	0.9364	$0.8911$	$0.9119$
	PMI	0.9467	0.7935	0.7764	0.6640	0.7952
	Counts	0.8348	$0.9389$	$0.9432$	0.6777	0.8487

Table 2. Table 2 : Results of Cox regression analysis with cluster membership indicators.

Embedding	Estimated HR	Lower 95% CI	Upper 95% CI	p-value
Fourier-Eigen	1.344	1.17	1.545	3.04e-05∗∗∗
PMI	0.9586	0.8805	1.044	0.33
BERT	0.8237	0.6741	1.006	0.07

Equations

$${\lambda }^N(t)dt=\mathbb{P}(dN(t)=1|{\mathcal{H}}_t),$$

$$V_{ij}(\tau )=\frac{\mathbb{E}[d{N}_i(t)d{N}_j(t-\tau )]}{dtd(t-\tau )}-{\Lambda }_i{\Lambda }_j-{\delta }_{ij}(\tau ){\Lambda }_i,i,j=1,\dots ,d,$$

$$G_i\sim \pi (\mathcal{G}),$$

$${\mathbf{M}}_i|G_i=g\sim \text{PoisProc}({\mu }^{(g)},T_i),$$

$$\displaystyle\lambda^{N_{ij}}(t)=\nu_{j}(\bm{X}_{i})+\sum_{\ell=1}^{k}\big{(}\omega_{j\ell}*dM_{i\ell}\big{)}(t).$$

$${\lambda }^{N_{ij}}(t)$$

$$={\nu }_j({\mathbf{X}}_i)+\sum _{\ell =1}^k{\int }_0^t{\omega }_{j\ell }(t-u)d{M}_{i\ell }(u)$$

$$={\nu }_j({\mathbf{X}}_i)+\sum _{\ell =1}^k\sum _{\text{Miℓ occurs at time u<t}}{\omega }_{j\ell }(t-u).$$

$$V_{j{j}^{\prime }}^{{\mathbf{N}}_i}(\tau )=\sum _{\ell =1}^k{\mu }_{i\ell }\cdot ({\omega }_{j\ell }\ast {\omega }_{j^{\prime }\ell })(\tau ).$$

$$V^{\bm{N}_{i}}(\tau)=\big{(}\mathbf{\omega}\star\mathbf{D}_{i}\star\mathbf{\omega}^{\top}\big{)}(\tau).$$

$$\mathsf{F}\big{\{}V^{\bm{N}_{i}}_{jj^{\prime}}\big{\}}(\xi)=\sum_{\ell=1}^{k}\mu_{i\ell}\cdot\bar{\omega}_{j\ell}(\xi)\cdot\bar{\omega}_{j^{\prime}\ell}(\xi)$$

$$\mathsf{F}\big{\{}V^{\bm{N}_{i}}\big{\}}(\xi)=\mathbf{W}{(\xi)}\mathbf{D}_{i}(\mathbf{W}{(\xi)})^{\top}.$$

$$\mathsf{F}\big{\{}V^{\bm{N}_{i}}\big{\}}(\xi_{0})=\mathbf{U}_{i}\Lambda_{i}\mathbf{U}_{i}^{-1},$$

$$f_i=({\lambda }_{i1},\dots ,{\lambda }_{ik}{)}^{\top }$$

$$\mathsf{F}\big{\{}\widehat{V}^{\bm{N}_{i}}\big{\}}(\xi_{0})=\widehat{\mathbf{U}}_{i}\widehat{\Lambda}_{i}\widehat{\mathbf{U}}_{i}^{-1},$$

$${\hat{f}}_i=({\hat{\lambda }}_{i1},\dots ,{\hat{\lambda }}_{ik}{)}^{\top },$$

$$\widehat{h}=\mathop{\mathrm{argmin}}_{h\in\mathcal{H}}\frac{1}{n}\sum_{i=1}^{n}L\big{(}G_{i},h(\widehat{f}_{i})\big{)},$$

$$\widehat{h}=\mathop{\mathrm{argmin}}_{h\in\mathcal{H}}\frac{1}{n}\sum_{i=1}^{n}L\big{(}G_{i},h(\widehat{f}_{i})\big{)}.$$

$$\widehat{V}^{\bm{N}_{i}}_{jj^{\prime}}(\tau)=\begin{cases}\frac{1}{Th}\int\int_{[0,T]^{2}}K\big{(}\frac{(t^{\prime}-t)+\tau}{h}\big{)}dN_{j^{\prime}}(t^{\prime})dN_{j}(t)-\frac{1}{T^{2}}N_{j^{\prime}}([0,T])\frac{1}{T^{2}}N_{j}([0,T])&j\neq j^{\prime}\\ \frac{1}{Th}\int\int_{[0,T]^{2}/\{t=t^{\prime}\}}K\big{(}\frac{(t^{\prime}-t)+\tau}{h}\big{)}dN_{j^{\prime}}(t^{\prime})dN_{j}(t)-\frac{1}{T^{2}}\big{(}N_{j}([0,T])\big{)}^{2}&j=j^{\prime}\end{cases},$$

$${\mathcal{V}}_{j{j}^{\prime }}^{(g)}(\xi )=\sum _{\ell =1}^k{\mu }_{\ell }^{(g)}\cdot {\bar{\omega }}_{j\ell }(\xi )\cdot {\bar{\omega }}_{j^{\prime }\ell }(\xi )$$

$$\frac{1}{k}\Vert f^{(g)}-f^{(r)}{\Vert }_2$$

$$\sum _{\ell =1}^k{\int }_0^{\infty }|{\omega }_{j\ell }(u)|du\le \gamma .$$

$$|L(g,r)-L(g,r^{\prime })|\le {\lambda }_{\text{Lip}}{\rho }_{\mathcal{G}}(r,r^{\prime }),$$

$$\text{diam}(\mathcal{G})=\underset{g,r\in \mathcal{G}}{\sup }{\rho }_{\mathcal{G}}(g,r).$$

$$h_{\text{oracle}}(x)={\mathrm{argmin}}_{g\in \mathcal{G}}\Vert f^{(g)}-x{\Vert }_2^2.$$

$$\mathcal{E}_{\text{approx.}}\big{(}\mathcal{H}\big{)}=\inf_{h\in\mathcal{H}}\sup_{x\in\mathbb{R}^{k}}\rho_{\mathcal{G}}\big{(}h(x),h_{\text{oracle}}(x)\big{)},$$

$$\mathcal{R}(L\circ\mathcal{H})=\mathbb{E}\sup_{h\in\mathcal{H}}\Big{|}\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}L\big{(}G_{i},h(\widehat{f}_{i})\big{)}\Big{|},$$

$$\|\mu^{(g)}-\mu^{(r)}\|_{2}\geq C\frac{d}{\sqrt{k}\sigma_{k}^{2}(\mathbf{W})}T^{-1/5}+3\frac{\sigma_{1}^{2}(\mathbf{W})}{k}(\|\mu^{(g)}\|_{2}\land\|\mu^{(r)}\|_{2})\cdot(\rho_{gr}\lor\rho^{2}_{gr})\big{)},$$

$$\mathbb{E}\Big{[}L\big{(}G,\widehat{h}(\widehat{f})\big{)}\Big{]}\leq\lambda_{\text{Lip}}\mathcal{E}_{\text{approx.}}(\mathcal{H})+2c_{2}d^{2}T^{6/5}\exp\big{(}c_{3}T^{-1/5}\big{)}\cdot\text{diam}(\mathcal{G})+4\cdot\mathcal{R}(L\circ\mathcal{H}).$$

$$\mathbb{E}\Big{[}L\big{(}G,\widehat{h}(\widehat{f})\big{)}\Big{]}\leq\lambda_{\text{Lip}}\mathcal{E}_{\text{approx.}}(\mathcal{H})+2c_{2}d^{2}T^{6/5}\exp\big{(}c_{3}T^{-1/5}\big{)}+C^{\prime}\sqrt{\frac{\text{VC}(\mathcal{H})}{n}}.$$

$$\|\mu^{(g)}-\mu^{(r)}\|_{2}\geq C\frac{d|\mathcal{G}|}{\sigma_{k}^{2}(\mathbf{W})}\sqrt{\frac{n}{k}}\left(\frac{1}{\sqrt{n_{g}}}+\frac{1}{\sqrt{n_{r}}}\right)T_{\min}^{-1/5}+3\frac{\sigma_{1}^{2}(\mathbf{W})}{k}(\|\mu^{(g)}\|_{2}\land\|\mu^{(r)}\|_{2})\cdot(\rho_{gr}\lor\rho^{2}_{gr})\big{)},$$

$$\underset{\sigma \in {\mathcal{S}}_{\mathcal{G}}}{\inf }\left[\sum _{i=1}^{n}1(G_i\ne \sigma ({\hat{G}}_i))\right]=0$$

Full text

Keywords: High-dimensional point process, EHR, Spectral methods, Embedding.

1 Introduction

The growing accessibility of electronic health record (EHR) datasets has beget new paradigms of public health research and precision medicine. Recent studies have demonstrated the utility of EHRs for building risk prediction models (Goldstein et al., 2016; Rasmy et al., 2021; Zhang et al., 2019; Corey et al., 2018), discovering clinically meaningful patient subgroups (Guo et al., 2023; Jiang et al., 2024; Xu et al., 2020; Mullin et al., 2021), and designing optimal individual treatment plans (Wang et al., 2016; Reddy et al., 2019; Zhou et al., 2024; Chen et al., 2025), all of which carry substantial potential to advance clinical decision making (Abul-Husn and Kenny, 2019; Giordano et al., 2021). These successes have spurred the development of new statistical methods designed to address challenges unique to EHR analysis, such as multi-institutional integration (Duan et al., 2020; Zhou et al., 2024; Xiong et al., 2023a; Tong et al., 2022; Li and Duan, 2024; Hu et al., 2024), semi-supervised learning due to limited gold-standard labels (Gronsbell and Cai, 2018; Cai et al., 2025; Zhang et al., 2023; Wang et al., 2025; Nogues et al., 2022; Zhang et al., 2022; Hou et al., 2023), and robustness under patient distribution shifts (Wang et al., 2022; Zhan et al., 2024; Knight et al., 2025; Lee et al., 2023; Wang et al., 2023; Xiong et al., 2023b; Luo et al., 2022; Li et al., 2024b).

However, a fundamental aspect of structured EHR data remains underutilized: its high-dimensional temporal structure. Structured EHRs consist of event times of thousands of billing codes, representing diagnoses, prescriptions, laboratory results, and procedures. Most existing methods collapse event times into vectors of counts over fixed observation windows, and even those that account for time typically do so using only crude temporal summaries. Such strategies do not effectively capture temporal information that may be strongly associated with clinical outcomes (Yuan et al., 2021). This has inspired a growing field of interdisciplinary research aiming to design models and methods for integrating temporal information into the analysis of EHRs. One recent line of work endeavors to model EHRs as a high-dimensional Hawkes process (Hawkes, 1971; Choi et al., 2015; Bao et al., 2017; Alaa et al., 2017; Sun et al., 2024; Zhao and Engelhard, 2025), which enables researchers to learn time-informed interaction networks between large numbers of codes by estimating a high-dimensional conditional intensity function. A parallel literature leverages the temporal information in EHRs to learn medical concept embeddings using contemporary tools from machine learning (ML) (Si et al., 2021), including recurrent neural networks (Ruan et al., 2019; Rasmy et al., 2022; Zhu et al., 2023; Al Olaimat et al., 2024), transformer architectures (Li et al., 2020; Choi et al., 2015; Yang et al., 2023; Yu et al., 2025), and foundation-scale large language models (Yang et al., 2022; Steinberg et al., 2021; Agrawal et al., 2022; Li et al., 2024a). These Hawkes process and ML based methods are able to learn complex relationships between medical concepts over time from the health record, but do so without providing theoretical guarantees, leaving the landscape of statistical theory for time-informed EHR methods underexplored. Recent work has begun to address this by adapting statistical models from text and sequence analysis (Arora et al., 2016), in particular log-linear hidden Markov models (Lu et al., 2023; Xu et al., 2023; Zhou et al., 2023), which capture temporal dependence between codes with simple co-occurrence statistics for building medical concept embeddings and knowledge graphs (Beam et al., 2020; Hong et al., 2021; Gan et al., 2025). The authors of these works are able to prove precise and informative theorems demonstrating the statistical properties of their methods. However, they do so by imposing restrictive parametric assumptions on the code occurrence probabilities and assuming discrete and evenly spaced event times, at odds with the irregularity of real EHR data. Most importantly, all of these methods are designed to capture concept-level relationships, not patient-level representations, which severely detracts from their utility in downstream tasks of high clinical importance, such as risk prediction and unsupervised patient stratification.

To bridge this gap, we propose a new model for temporal EHR data: the latent factor point process model. In this model, the code event times are represented as a high-dimensional point process whose conditional intensity is driven by a low-dimensional, unobserved latent factor Poisson process. This structure is inspired by clinical reality: although thousands of heterogeneous codes are recorded in EHRs, their co-occurrence is often driven by a small number of underlying medical processes (e.g., disease progression, comorbidity conditions, treatment response). In addition to clinical interpretability, the low-rank structure of our proposed model is statistically essential for estimating parameters of interest in high dimensions.

Motivated by patient-level tasks such as classification and clustering, we assume that patients belong to distinct subgroups and that the latent intensity vector varies systematically across these groups. In contrast to the linear Hawkes process and other parametric point process models (Daley and Vere-Jones, 2003; Bacry et al., 2020), we do not impose specific functional forms on transfer functions, allowing the event probabilities to flexibly evolve with the latent history. We show that the cross-covariance matrix of the observed codes inherits a low-rank structure, and that the spectral properties of this matrix differ systematically across patient subgroups. This motivates our novel Fourier-Eigen embedding, a time-informed patient representation constructed from the spectral density matrix of the observed process.

We provide a rigorous theoretical analysis showing that Fourier-Eigen embeddings serve as efficient patient representations for downstream tasks, including empirical risk minimization in classification and subgroup discovery in clustering. We supplement our analysis with simulation studies and an application to an Alzheimer’s disease cohort, a clinically important setting where capturing temporal progression is critical for understanding heterogeneous disease trajectories.

1.1 Related literature

Our work builds on recent developments in the high-dimensional point process literature (Hansen et al., 2015; Bacry and Muzy, 2016; Hawkes, 2018). The latent factor point process model that we propose is inspired by the multivariate temporal point process regression model of Tang and Li (2023), where the conditional intensity of an outcome process is driven entirely by a high-dimensional predictor process. In our setting, the predictor process corresponds to the latent factor process and the outcome process to the observed code occurrence process. A key difference is that our latent factor process is unobserved, which makes the estimation strategy of Tang and Li (2023) inapplicable. Wang and Shojaie (2021) also considers a high-dimensional point process model in which a subset of the component processes are unobserved yet still affect the conditional intensity of the observed processes. Under a linear assumption, they provide a deconfounded estimator of the intensity coefficients, which is consistent when the bias from unmeasured components is negligible. In contrast, our model assumes that the conditional intensity of the observed process is entirely driven by an unobserved process, and we do not require the intensity to follow a particular functional form.

The construction of our Fourier-Eigen embeddings builds on spectral methods for point processes. The spectral density matrix, defined as the Fourier transform of the cross-covariance matrix, has been widely used in classical work on point processes (Bartlett, 1963; Hawkes, 1971; Bacry and Muzy, 2016) and time series in general (Koopmans, 1995). For example, Bacry and Muzy (2016) derive Wiener–Hopf equations for a multivariate Hawkes process, which they solve in the Laplace domain to obtain a nonparametric estimator of the conditional intensity. More recently, Pinkney et al. (2024) study the multi-taper periodogram for estimating high-dimensional spectral density matrices and propose a regularized inverse spectral estimator, while Hellstern et al. (2025) develop an estimator for the difference between two spectral density matrices in differential network analysis. These advances highlight the utility of spectral approaches for extracting structure from complex stochastic processes, and we adapt this perspective to define low-dimensional patient embeddings.

Finally, our work contributes to the growing literature on statistical models for temporal EHR data. The seminal work of Arora et al. (2016) analyzed a log-linear hidden Markov model for text generation, showing that the pointwise mutual information matrix captures latent embeddings of words up to a tolerable error. This idea has since been adapted to EHR analysis: for instance, Lu et al. (2023) allow the code embeddings to admit a block-wise covariance structure for estimating structured knowledge graphs of medical concepts, while Xu et al. (2023) provide edge-wise inference on such graphs. Related extensions include multi-institutional adaptations (Zhou et al., 2023) and code-level clinical knowledge extraction (Beam et al., 2020; Hong et al., 2021; Gan et al., 2025). More recently, Yu et al. (2025) propose a marked point process model with time-aware attention to learn code embeddings. These works are primarily designed for capturing relationships between codes. By contrast, our goal is to develop representations tailored to patient-level tasks such as classification and clustering, while flexibly incorporating the temporal structure of EHR data.

1.2 Notation

For two real numbers $a,b\in\mathbb{R}$, we let $a\land b=\min(a,b)$ and $a\lor b=\max(a,b)$. For a matrix $\mathbf{A}\in\mathbb{C}^{n\times m}$, we let $\{\sigma_{i}(\mathbf{A})\}_{i=1}^{n\land m}$ denote its singular values in decreasing order. We use $\|\mathbf{A}\|_{F}=\sqrt{\sum_{i,j}|A_{ij}|^{2}}$ and $\|\mathbf{A}\|_{\text{op}}=\sigma_{1}(\mathbf{A})$ to denote the Frobenius and operator norms of $\mathbf{A}$ respectively. We will use $\mathbf{A}^{\dagger}$ to denote the conjugate transpose of $\mathbf{A}$, and we say that the matrix $\mathbf{A}$ is Hermitian if $\mathbf{A}=\mathbf{A}^{\dagger}$. For two functions $f,g:\mathbb{R}\to\mathbb{R}$, we write their convolution as $(f*g)(t)=\int f(\tau)g(t-\tau)\textrm{d}\tau$. If $\mathbf{A}(t)\in\mathbb{R}^{n\times m}$ and $\mathbf{B}(t)\in\mathbb{R}^{m\times q}$ are matrix-valued functions with entries $A_{ij}(t)$ and $B_{j\ell}(t)$, we define \big{(}\mathbf{A}\star\mathbf{B}\big{)}(t) as the matrix with entries $(A\star B)_{i\ell}(t)=\sum_{j=1}^{m}(A_{ij}*B_{j\ell})(t)$. For an absolutely integrable function $f:\mathbb{R}\to\mathbb{R}$, we let its Fourier transform evaluated at frequency $\xi$ be $\mathsf{F}\{f\}(\xi)=\int_{\mathbb{R}}f(x)\exp(-i2\pi\xi x)dx$. If $\mathbf{A}(t)\in\mathbb{R}^{n\times m}$ has integrable elements $A_{ij}(t)$, then $\mathsf{F}\{\mathbf{A}(t)\}(\xi)$ denote the matrix with entries $\mathsf{F}\{A_{ij}\}(\xi)$.

A univariate temporal point process $N$ is a stochastic process on $\mathbb{R}^{+}$ with event times $\{t_{\ell}:\ell=1,2,\dots\}$. We write $N(A)=\sum_{\ell}{\bm{1}}(t_{\ell}\in A)$ for $A\subseteq\mathbb{R}^{+}$. and $dN(t)=N([t,t+dt])$ for a small increment $dt$. The mean intensity of $N$ is $\Lambda(t)=\lim_{dt\rightarrow 0}\mathbb{E}[dN(t)]/dt$. We say $N$ is stationary if the distribution of $dN(t)$ does not depend on $t$. The conditional intensity function of $N$ is defined by

[TABLE]

where $\mathcal{H}_{t}$ is the history of $N$ up to time $t$ (i.e. the $\sigma$-algebra generated by ${N(u):0<u<t}$).

We define a multivariate point process as $\bm{N}=(N_{1},...,N_{d})^{\top}\in\mathbb{R}^{d}$ with each element being an univariate point processes. If $\bm{N}$ is stationary, we define its cross-covariance matrix as

[TABLE]

where $\tau\in\mathbb{R}$, $\Lambda_{i}$ is the mean intensity of $N_{i}$, and $\delta_{ij}(\tau)={\bm{1}}(i=j)\delta(\tau)$ with $\delta(\cdot)$ the Dirac delta function.

Finally, for $\mu\in\mathbb{R}^{k}$, we write $\bm{M}\sim\text{PoisProc}(\mu,T)$ if $\bm{M}$ is a homogeneous Poisson process on $\mathbb{R}^{k}$ with intensity vector $\mu$ observed over $[0,T]$. Throughout, $c_{1},c_{2},c_{3},c_{4},C,C^{\prime}$ denote positive constants whose values may vary across uses.

2 Model and methods

2.1 Latent factor point process model

We let $\bm{N}$ denote a $d$-dimensional point process, $\bm{X}\in\mathbb{R}^{p}$ a vector of baseline covariates, $\bm{M}$ an unobserved $k$-dimensional point process whose meaning will be made clear in the sequel, and $G$ is a class indicator that takes values in a finite set $\mathcal{G}$. The full dataset consists of $n$ i.i.d. quartets $\mathcal{D}=\{(\bm{N}_{i},\bm{X}_{i},\bm{M}_{i},G_{i})\}^{n}_{i=1}$, but we only assume that we will have access to either $\mathcal{O}_{\text{sup}}=\{(\bm{N}_{i},\bm{X}_{i},G_{i})\}^{n}_{i=1}$ or $\mathcal{O}_{\text{unsup}}=\{(\bm{N}_{i},\bm{X}_{i})\}^{n}_{i=1}$; in either case, we observe the process $\bm{N}_{i}$ from time [math] to $T_{i}>0$ for each $i\in\{1,\ldots,n\}$. The subscripts “sup” and “unsup” refer to the supervised and unsupervised settings, respectively, and indicate whether or not the labels $\{G_{i}\}_{i=1}^{n}$ are observed. In the context of EHR data, $\bm{N}_{i}$ denotes the occurrence times for $d$ EHR codes observed for patient $i$, and $G_{i}$ may denote a disease status indicator. Our overarching goal is to accurately distinguish the patients belonging to the classes in $\mathcal{G}$ using the observed data, which may or may not contain labels $\{G_{i}\}_{i=1}^{n}$. Our first step in this direction is to articulate a new generative model for high-dimensional point processes that simultaneously describes the correlation structure between codes that we observe in real-world EHR data and captures differences between the subgroups of patients. We will impose the following model on the process $\bm{N}_{i}$:

[TABLE]

Here, $\pi$ is an arbitrary distribution on $\mathcal{G}$ that assigns nonzero mass to each $g\in\mathcal{G}$. As is standard in the point process literature (Hansen et al., 2015), we state our model for $\bm{N}_{i}$ in (2.3) in terms of the conditional intensity function $\lambda^{N_{ij}}$. We let $\bm{M}_{i}$ denote an $k$-dimensional unobserved, latent homogeneous Poisson process whose intensity vector $\mu^{(g)}$ is determined by the latent class $G_{i}$. The choice of a homogeneous Poisson latent driver is a parsimonious assumption that induces a low-rank second-order structure while keeping estimation feasible in high dimensions; extensions to inhomogeneous or self-exciting latent processes can introduce additional analytical challenges, and we leave their exploration to future work.

In (2.3), $\nu_{j}:\mathbb{R}^{p}\rightarrow\mathbb{R}$ is the background intensity of feature $j$ which is a function of baseline covariates, and $\omega_{j\ell}:\mathbb{R}^{+}\mapsto\mathbb{R}^{+}$ is a transfer function that controls the extent to which the history of the unobserved process $M_{i\ell}$ affects the intensity of feature $j$. We can gain intuition for the role that $\bm{M}_{i}$ plays in the conditional intensity of $\bm{N}_{i}$ by unwinding the convolution and the Riemann-Stieltjes integral in (2.3) as follows:

[TABLE]

In (2.6), we can see that each occurrence of $M_{i\ell}$ triggers an increase in the probability of $N_{ij}$ occurring through the transfer function $\omega_{j\ell}$. Analogous to the structure of latent factor models used in the time series analysis literature (Lam et al., 2011; Molenaar, 1985; Bartholomew et al., 2011), the process $\bm{M}_{i}$ acts as the latent factor process that drives the occurrences observed in $\bm{N}_{i}$. For this reason, we say that $\bm{N}_{i}$ follows the latent factor point process model with latent factors $\bm{M}_{i}$. We emphasize that the transfer functions $\{\omega_{j\ell}\}$ are not dependent on the index $i$, and are shared in common across all patients. The history of patient $i$ affects $\lambda^{N_{ij}}$ only through $\bm{M}_{i}$, whose intensity is determined by $G_{i}$.

Model (2.3) is inspired by the structure of electronic health record data. The occurrence of EHR codes is caused not by the prior occurrence of other codes, but by real-world medical events that are represented in the health record only by proxy. For instance, a myocardial infraction (a heart attack) is denoted in the health record by several different ICD-10 codes (Rockenschaub et al., 2020), which will in turn co-occur with codes denoting medical procedures and treatments for patients who have experienced a heart attack. The latent process $\bm{M}_{i}$ can be interpreted as these real-world events whose occurrence is observed only through the ambient process $\bm{N}_{i}$. In light of this interpretation, equation (2.2) states that the rates of the real-world latent events differ systematically between the classes of patients, while the transferring functions that map these events to the health record are shared across all patients. The intuition for the role of $\bm{M}_{i}$ as representing real-world events also motivates our specification of the latent process as a Poisson process rather than, say, a Gaussian process as used in Lu et al. (2023) and Xu et al. (2023).

Furthermore, Model (2.3) grants a low rank structure in the cross-covariance matrix of the process $\bm{N}_{i}$ that interpretably encodes heterogeneity between the groups of patients. The following proposition, which is a special case of Proposition 1 in (Tang and Li, 2023), describes this formally. We provide a proof in the supplement for completeness.

Proposition 2.1**.**

Suppose that $\bm{N}_{i}\in\mathbb{R}^{d}$ and $\bm{M}_{i}\in\mathbb{R}^{k}$ are stationary point processes that satisfy (2.3). Let $\mu_{i}\in\mathbb{R}^{k}$ denote the mean intensity of the process $\bm{M}_{i}$ with entries $\mu_{i\ell}=\lim_{dt\rightarrow 0}\mathbb{E}[dM_{i\ell}(t)]/dt$, and assume that $M_{i\ell}$ and $M_{i\ell^{\prime}}$ are independent for $\ell\neq\ell^{\prime}$. Then the elements of the instantaneous cross-covariance matrix of $\bm{N}_{i}(t)$ at lag $\tau\neq 0$ are given by

[TABLE]

Defining $t\mapsto\mathbf{D}_{i}(t)$ as the constant function that returns $\text{diag}(\mu_{i})$, we can thus express $V^{\bm{N}_{i}}$ as

[TABLE]

This low-rank form is the key property that enables us to compress patient trajectories into a small number of spectral features, which underpins the Fourier-Eigen embedding developed below. The form of $V^{\bm{N}_{i}}$ given in Equation (2.8) is analogous to its spectral decomposition at $\tau$ with ‘spectrum’ $\mu_{i}$. We can make this analogy explicit with elementary tools from harmonic analysis. With suitable regularity conditions on the transfer functions $\{\omega_{j\ell}\}$, we invoke the Convolution Theorem (Dym and McKean, 1972) to convert the structure given in (2.8) to that of a standard matrix product using a Fourier transform. This gives us

[TABLE]

where we define $\bar{\omega}_{j\ell}(\xi)=\mathsf{F}\{\omega_{j\ell}\}(\xi)$ as the Fourier transform of $\omega_{j\ell}$ at evaluated at $\xi$. Letting $\mathbf{W}{(\xi)}\in\mathbb{C}^{d\times k}$ denote the matrix with entries $W{(\xi)}_{j\ell}=\bar{\omega}_{j\ell}(\xi)$, we can thus write

[TABLE]

Although the matrix $\mathbf{D}_{i}$ is diagonal, this factorization of \mathsf{F}\big{\{}V^{\bm{N}_{i}}\big{\}}(\xi) does not constitute an eigendecomposition, as we have no guarantee that $\mathbf{W}{(\xi)}(\mathbf{W}{(\xi)})^{\top}=\bm{I}_{d}$ and $(\mathbf{W}{(\xi)})^{\top}\mathbf{W}{(\xi)}=\bm{I}_{k}$. However, since $\mathbf{W}(\xi)$ spans the column space of \mathsf{F}\big{\{}V^{\bm{N}_{i}}\big{\}}(\xi), the leading eigenvalues retain the subgroup-specific information carried by $\mathbf{D}_{i}$. In particular, the top $k$ eigenvalues of \mathsf{F}\big{\{}V^{\bm{N}_{i}}\big{\}}(\xi) can be expressed as $\mathbf{H}\mathbf{D}_{i}$ for a suitable invertible matrix $\mathbf{H}$ as long as $\mathbf{W}{(\xi)}$ has full column rank. Recall that $\mathbf{D}_{i}=\text{diag}(\mu_{i})$. Under (2.2), this indicates that the eigenvalues of \mathsf{F}\big{\{}V^{\bm{N}_{i}}\big{\}}(\xi) and \mathsf{F}\big{\{}V^{\bm{N}_{i^{\prime}}}\big{\}}(\xi) will differ systematically when $G_{i}\neq G_{i^{\prime}}$. This explains our use of the term Fourier-Eigen in the sequel.

2.2 Classification and clustering with Fourier-Eigen embeddings

This exploration of the structure of $V^{\bm{N}_{i}}$ motivates our proposed procedure for building a classifier or clustering algorithm. We propose to estimate a low-dimensional patient embedding that captures the signal needed to distinguish the classes of patients, and use these embeddings as features in the downstream tasks of classification and clustering. The embeddings are defined as follows. Given a fixed frequency value $\xi_{0}\in\mathbb{R}$, for each $i\in[n]$ we perform the eigendecomposition

[TABLE]

where $\Lambda_{i}$ is the diagonal matrix of eigenvalues $\lambda_{i1},\ldots,\lambda_{id}$. We then define

[TABLE]

as the $k$-dimensional Fourier-Eigen embedding of patient $i$. In Section 3, we will show formally that if patients $i$ and $i^{\prime}$ belong to different classes then their embeddings $f_{i}$ and $f_{i^{\prime}}$ are well-separated, and therefore a standard classifier trained on the pairs $\{(f_{i},G_{i})\}^{n}_{i=1}$ or a suitable clustering algorithm will be able to perfectly distinguish the classes of patients. Of course, $f_{i}$ is a population-level quantity, as its construction depends on knowledge of the true cross-covariance matrix $V^{\bm{N}_{i}}$. We estimate $f_{i}$ by repeating this Fourier-eigendecomposition procedure with an estimator of $V^{\bm{N}_{i}}$. Letting $\widehat{V}^{\bm{N}_{i}}$ denote such an estimator, we compute the eigendecomposition

[TABLE]

and define

[TABLE]

where $\widehat{\lambda}_{i1},\ldots,\widehat{\lambda}_{id}$ denote the eigenvalues of \mathsf{F}\big{\{}\widehat{V}^{\bm{N}_{i}}\big{\}}(\xi_{0}). The vector $\widehat{f}_{i}$ will serve as our estimate of the Fourier-Eigen embedding of patient $i$. We note that the vectors $\{\widehat{f}_{i}\}_{i=1}^{n}$ obtained from Algorithm 1 are guaranteed to be in $\mathbb{R}^{k}$ as long as $\tau\mapsto\widehat{V}^{\bm{N}}(\tau)$ is symmetric in $\tau$, as this ensures that the matrix \mathsf{F}\big{\{}\widehat{V}^{\bm{N}_{i}}\big{\}}(\xi_{0}) is Hermitian. We summarize this procedure in Algorithm 1. From here, we proceed according to whether or not we observe class labels for each of the $n$ patients. If we observe the labeled dataset $\mathcal{O}_{\text{sup}}$ we train a classifier on the pairs $\{(\widehat{f}_{i},G_{i})\}^{n}_{i=1}$ with Empirical Risk Minimization, computing

[TABLE]

where $\mathcal{H}$ is a function class on $\mathbb{R}^{k}$ with outputs in $\mathcal{G}$. The loss function $L:\mathcal{G}^{\otimes 2}\rightarrow\mathbb{R}$ is chosen by the statistician. Otherwise, if we only have access to $\mathcal{O}_{\text{unsup}}$, we run K-means clustering with a spectral initialization (Kumar and Kannan, 2010) on the estimated embeddings $\{\widehat{f}_{i}\}_{i=1}^{n}$. We state the classification procedure in Algorithm 2 and the clustering procedure in Algorithm 3. To estimate the patient-level cross covariance matrix $V^{\bm{N}_{i}}$, we use the kernel smoothing estimator proposed by (Chen et al., 2017) as described in Algorithm 4.

Remark 2.2**.**

Our proposed embedding procedure in Algorithm 1 depends on a choice of frequency value $\xi_{0}$. The choice of $\xi_{0}$ does not matter in general, since under (2.2) the vectors $\{\mu^{(g)}:g\in\mathcal{G}\}$ are assumed to be constant in time and hence are unchanged by taking the Fourier transform. However, we do need the matrix $\mathbf{W}{(\xi_{0})}$ to be full rank to ensure that the Fourier-Eigen embeddings corresponding to distinct classes of patients are linearly separated; see Lemma 3.4. In practice, choosing $\xi_{0}$ at a moderate frequency yields stable embeddings, while extremely low or high frequencies can be noisier.

Remark 2.3**.**

Implicitly, Algorithm 1 assumes knowledge of $k$, the dimension of the latent processes $\{\bm{M}_{i}\}_{i=1}^{n}$. We do not explore estimation of $k$ in the present work, but in practice, users of Algorithm 1 can specify $k$ with standard methods for estimating the number of nonzero eigenvalues of a matrix (Fan et al., 2022).

Remark 2.4**.**

Since our embedding method only relies on the cross-covariance matrices of the observed processes $\{\bm{N}_{i}\}_{i=1}^{n}$, we avoid the issue of estimating the baseline intensity functions $\nu_{j}:\mathbb{R}^{p}\rightarrow\mathbb{R},j=1,\ldots,d$. We treat $\nu_{j}(\bm{X}_{i})$ as a nuisance component that absorbs baseline covariate effects and focus on temporal dependence; explicitly incorporating $\bm{X}_{i}$ into the embedding is an interesting extension left for future work.

3 Theoretical analysis

In this section, we give a theoretical analysis of our Fourier-Eigen embeddings and their use in downstream classification and classification tasks.

Assumption 3.1**.**

The full data $\mathcal{D}=\{(\bm{N}_{i},\bm{X}_{i},\bm{M}_{i},G_{i})\}^{n}_{i=1}$ are drawn i.i.d. from the latent factor point process model (2.1 - 2.3). For each $i\in[n]$, the processes $M_{i\ell}$ and $M_{i\ell^{\prime}}$ are independent for $\ell\neq\ell^{\prime}$. For each $i\in[n]$, the processes $\bm{N}_{i}$ and $\bm{M}_{i}$ are stationary. The transfer functions $\{\omega_{j\ell}\}$ are Lebesgue measurable and absolutely integrable.

Lemma 3.2**.**

Suppose that Assumption 3.1 holds. For $g\in\mathcal{G}$, let $\mathcal{V}^{(g)}(\xi)=\mathsf{F}\{V^{\bm{N}_{i}}\}(\xi)$, where $\bm{N}_{i}$ is drawn from (2.3) with $G_{i}=g$. Then at each $\xi\in\mathbb{R}$, the matrix $\mathcal{V}^{(g)}(\xi)$ is Hermitian with entries given by

[TABLE]

Lemma 3.2 guarantees that the population-level Fourier-Eigen embeddings are real-valued, and formalizes our heuristic discussion after Proposition 2.1. This allows us to guarantee that the population level Fourier-Eigen embeddings for distinct classes of patients are well-separated in $\ell_{2}$ norm under the following assumption.

Assumption 3.3**.**

There exists a frequency value $\xi^{*}\in\mathbb{R}$ at which the matrix $\mathbf{W}{(\xi^{*})}$ has full column rank.

Lemma 3.4**.**

Suppose that Assumptions 3.1 and 3.3 hold. For $g\in\mathcal{G}$, let $f^{(g)}$ denote the population level Fourier-Eigen embedding obtained at frequency $\xi^{*}$ for a patient in class $g$ and $\mathbf{U}^{(g)}\in\mathbb{C}^{d\times k}$ denotes the first $k$ eigenvectors of $\mathcal{V}^{(g)}(\xi^{*})$. Let $g,r\in\mathcal{G}$ denote an arbitrary pair of classes and define $\rho_{gr}=\|(\mathbf{U}^{(g)})^{\dagger}\mathbf{U}^{(r)}-\mathbf{I}\|_{F}$. Then the following inequality holds:

[TABLE]

Remark 3.5**.**

Lemma 3.4 allows us to translate signal between classes in the latent process intensity vectors to differences in the corresponding population-level Fourier-Eigen embeddings. This result formalizes the connection between our latent factor point process model in (2.1 - 2.3) and our construction of the Fourier-Eigen embeddings as a statistical tool. The inequality stated in the lemma provides a sufficient condition such that $\|f^{(g)}-f^{(r)}\|_{2}$ is large. The right-hand side is positive for general $\mu^{(g)}$ and $\mu^{(r)}$ as long as the eigenvectors of $\mathcal{V}^{(g)}(\xi^{*})$ and $\mathcal{V}^{(r)}(\xi^{*})$ are well-aligned as captured by the quantity $\rho_{gr}$. The only identification condition that we need to prove this lemma is Assumption 3.3, which asserts the existence of a frequency $\xi^{*}\in\mathbb{R}$ at which $\mathbf{W}(\xi^{*})$ is full rank. This is satisfied if, for instance, the transfer functions satisfy $\omega_{j\ell}(t)=a_{j\ell}\kappa(t)$, where $\kappa$ is a Lebesgue measurable function of time and the matrix $\mathbf{A}\in\mathbb{R}^{d\times k}$ with entries $A_{j\ell}=a_{j\ell}$ is full rank.

Assumption 3.6**.**

There exists a constant $b_{0}>0$ such that $\omega_{j\ell}(b)=0$ for $|b|\geq b_{0}$ for all $j\in[d]$ and $\ell\in[k]$.

Assumption 3.7**.**

There exists a positive constant $\gamma<1$ such that for all $j\in[d]$,

[TABLE]

Assumption 3.8**.**

The elements of the cross-covariance matrix $V^{\bm{N}}$ are $\theta$-Lipschitz functions.

Remark 3.9**.**

Assumptions 3.6 - 3.8 are required to apply the concentration properties of the kernel-smoothing estimator of the cross-covariance matrix given in Algorithm 4. These are standard in the analysis of high-dimensional point processes; see for instance Assumptions 2 and 3 of Chen et al. (2017). We remark that our Assumption 3.6 is strictly stronger than Assumption 2 of Chen et al. (2017), as we require that the tails of the transfer functions $\omega_{j\ell}$ are identically 0 at long time horizons. This is used to control the fluctuations of the Fourier-transformed cross-covariance matrix; see the appendix for details.

We now introduce the following notation to present our first theorem, which states a bound on the generalization error of a classifier trained on the estimated Fourier-Eigen embeddings with Algorithm 2 over the function class $\mathcal{H}$. Let $\rho_{\mathcal{G}}(\cdot,\cdot)$ denote a real-valued dissimilarity measure on the set of classes $\mathcal{G}$; for instance, $\rho_{\mathcal{G}}(g,r)={\bm{1}}(g\neq r)$. For a constant $\lambda_{\text{Lip}}>0$, we say that the loss function $L(\cdot,\cdot)$ is $\lambda_{\text{Lip}}$-Lipschitz if, for any fixed $g\in\mathcal{G}$, it holds

[TABLE]

for $r,r^{\prime}\in\mathcal{G}$. We define the diameter of the set $\mathcal{G}$ with respect to $\rho_{\mathcal{G}}$ as

[TABLE]

Additionally, we define the oracle classifier $h_{\text{oracle}}:\mathbb{R}^{k}\rightarrow\mathcal{G}$ as follows:

[TABLE]

Finally, we define

[TABLE]

as the approximation error of the function class $\mathcal{H}$. We will measure the complexity of $\mathcal{H}$ in terms of its Rademacher complexity (Bartlett and Mendelson, 2002) under composition with the loss function $L(\cdot,\cdot)$, which is defined as

[TABLE]

where $\varepsilon_{1},\ldots,\varepsilon_{n}$ are i.i.d. Rademacher random variables, and the pairs $\{(G_{i},\widehat{f}_{i})\}_{i=1}^{n}$ constitute the training data passed to Algorithm 2. The outer expectation is taken jointly over the Rademacher random variables and the training data.

Theorem 3.10**.**

Suppose that Assumptions 3.1, 3.3, 3.6, 3.7, and 3.8 hold. Given the data $\mathcal{O}_{\text{sup}}=\{(\bm{N}_{i},\bm{X}_{i},G_{i})\}_{i=1}^{n}$,, let $\{\widehat{f}_{i}\}_{i=1}^{n}$ denote the estimated Fourier-Eigen embeddings returned from Algorithm 1 at frequency $\xi^{*}$, where Algorithm 4 is implemented with bandwidth $h=c_{1}T_{i}^{-1/5}$ for each $i\in\{1,\ldots,n\}$ and threshold $C_{\text{tr}}\leq 2b_{0}$ for a constant $c_{1}>0$. Suppose that there exists a constant $\lambda_{\text{Lip}}>0$ such that $L(\cdot,\cdot)$ is $\lambda_{\text{Lip}}$-Lipschitz. Let $\widehat{h}$ denote the classifier trained on the pairs $\{(\widehat{f}_{i},G_{i})\}_{i=1}^{n}$ with Algorithm 2 using the function class $\mathcal{H}$ and loss function $L(\cdot,\cdot)$ and let $(\widehat{f},G)$ denote a new pair independently drawn from the same distribution as the training data with observation time $T>0$. Then there exist constants $c_{2},c_{3},C>0$ such that if for any pair of classes $g,r\in\mathcal{G}$, $\mu^{(g)}$ and $\mu^{(r)}$ satisfy

[TABLE]

then it holds

[TABLE]

The conclusion of Theorem 3.10 takes the form of standard risk bounds for empirical risk minimizers (Vapnik, 2013; Bartlett and Mendelson, 2002; Hajek and Raginsky, 2021). The first term on the right hand side of (3.10) encodes the ability of functions in $\mathcal{H}$ to approximate an oracle classifier that relies on knowledge of the population-level Fourier-Eigen embeddings, the second term denotes the worst-case error of the oracle classifier, and the final term accounts for the complexity of the function class $\mathcal{H}$. We can refine this result in the case of binary classification, as given in the following corollary. Here, we let $\text{VC}(\mathcal{H})$ denote the Vapnik-Chervonenkis dimension (Devroye et al., 1996) of $\mathcal{H}$.

Corollary 3.11**.**

Suppose that all of the conditions of Theorem 3.10 hold. Furthermore, suppose that $\mathcal{G}=\{0,1\}$ and $\rho_{\mathcal{G}}(g,r)=|g-r|$ for $g,r\in\mathcal{G}$. Then there exist constants $c_{2},c_{3},C,C^{\prime}>0$ such that if (3.9) holds for any pair of classes $g,r\in\mathcal{G}$, then

[TABLE]

Remark 3.12**.**

Corollary 3.11 allows us to evaluate the statistical complexity of binary classification for standard choices of the function class $\mathcal{H}$. For instance, it is well known that the class of linear classifiers $\mathcal{H}_{\text{lin}}$ on $\mathbb{R}^{k}$ satisfies $\text{VC}(\mathcal{H}_{\text{lin}})\asymp k$ (Hajek and Raginsky, 2021), and the class of feedforward neural networks $\mathcal{H}_{\text{FF}}$ with $w$ weights has $\text{VC}(\mathcal{H}_{\text{FF}})\lesssim w\log w$ (Baum and Haussler, 1988). This also emphasizes one of the advantages of using our Fourier-Eigen embeddings for classification rather than the raw EHR trajectory data, as the reduction in dimension granted by our embeddings can significantly reduce the size of $\text{VC}(\mathcal{H})$, especially if $k\ll d$.

We now turn our attention to the performance of Algorithm 3. In the following theorem, we let $n_{g}=\sum_{i=1}^{n}{\bm{1}}(G_{i}=g)$ for each $g\in\mathcal{G}$ and we define ${\mathcal{S}}_{\mathcal{G}}$ as the set of permutations on $\mathcal{G}$, meaning that $\sigma\in{\mathcal{S}}_{\mathcal{G}}$ implies that $\sigma$ is a bijection from $\mathcal{G}$ to $\mathcal{G}$.

Theorem 3.13**.**

Suppose that Assumptions 3.1, 3.3, 3.6, 3.7, and 3.8 hold. Given the data $\mathcal{O}_{\text{unsup}}=\{(\bm{N}_{i},\bm{X}_{i})\}_{i=1}^{n}$, let $\{\widehat{f}_{i}\}_{i=1}^{n}$ denote the estimated Fourier-Eigen embeddings returned from Algorithm 1 at frequency $\xi^{*}$, where Algorithm 4 is implemented with bandwidth $h=c_{1}T_{i}^{-1/5}$ for each $i\in\{1,\ldots,n\}$ and threshold $C_{\text{tr}}\leq 2b_{0}$ for a constant $c_{1}>0$. Finally, let $\{\widehat{G}_{i}\}_{i=1}^{n}$ denote the estimated cluster labels obtained from Algorithm 3 and $T_{\min}=\min_{i=1}^{n}T_{i}$. Then the following claim holds conditional on the labels $\{G_{i}\}_{i=1}^{n}$: there exist constants $c_{2},c_{3},C>0$ such that if for any pair of classes $g,r\in\mathcal{G}$, $\mu^{(g)}$ and $\mu^{(r)}$ satisfy

[TABLE]

then it holds

[TABLE]

with probability at least $1-2c_{2}nd^{2}T_{\min}^{6/5}\exp\left(-c_{3}T_{\min}^{1/5}\right)$.

Theorem 3.13 shows that a standard K-means clustering algorithm can perfectly recover the latent group membership of each patient, as long as the signal separating the groups is sufficiently strong. We point out that the signal strength condition (3.12) is stronger than (3.9) as required in Theorem 3.10. We interpret this gap as the “cost” of unsupervised learning, which is minimal when we observe patients from a small number of balanced classes.

4 Simulation studies

We conduct simulation studies to evaluate the performance of classification and clustering with our proposed Fourier-Eigen embeddings. We simulate data $\{\bm{N}_{i}\}_{i=1}^{n}$ from our latent factor point process model using the thinning algorithm (Daley and Vere-Jones, 2003) with two classes of patients of equal size. We fix the ambient dimension $d=100$ and the latent dimension $k=2$, and vary the total sample size $n\in\{250,500,750\}$ and observation time $T\in\{50,100,150\}$. We specify the transfer functions as $\omega_{j\ell}(t)=a_{j\ell}e^{-t^{2}/2}$, where the coefficients $a_{j\ell}$ are drawn from a $\text{Unif}(0,1/2)$ distribution. We vary the signal strength between the two groups, as encoded by the difference between the latent intensity vectors $\mu^{(0)}$ and $\mu^{(1)}$, as $\|\mu^{(0)}-\mu^{(1)}\|_{2}\in\{0.3,0.4,0.5,0.6,0.7,0.8\}$.

We compute our Fourier-Eigen embeddings with Algorithm 1 at frequency $\xi_{0}=1$, where the cross-covariance matrix is estimated with Algorithm 4 using bandwidth $h=1$, $C_{\text{tr}}=5$, and a Gaussian kernel function. We calculate the Fourier transform of the cross-covariance matrix by numerical integration. We compare the performance of the Fourier-Eigen embeddings on downstream tasks to the follow two embedding methods:

Code counts: For each $i\in\{1,\ldots,n\}$, we compute a feature vector in $\mathbb{R}^{d}$ by summing the number of occurrences of each code. 2. 2.

Point-wise mutual information (PMI) embeddings: We calculate the PMI matrix for each simulated process $N_{i}$ according to the formulation in Arora et al. (2016) with window size $T/20$. We use the first $k=2$ eigenvalues of the PMI matrix as the feature vector.

First we evaluate the performance of each of these three embedding methods for classification. We train three sets of standard logistic regression models using the class membership indicators as outcomes and the Fourier-Eigen embeddings, code counts, and PMI embeddings as features respectively. We evaluate the prediction performance of each embedding via out-of-sample AUC on a testing set of size $n_{\text{test}}=50$. The results are given for each configuration of $n$, $T$, and signal strength in Figure 1. The AUC is averaged over 100 replications.

We can see that the discriminative performance of Fourier-Eigen embeddings becomes especially apparent as the observation time $T$ increases. At $T=50$, the three feature embeddings admit very similar classification performance, with the code counts and Fourier-Eigen embeddings slightly outperforming the PMI embeddings. However, the logistic regression models trained on the PMI embeddings and code counts do not improve with increasing $T$, as these methods are unable to take advantage of increased temporal information due to model misspecification. Our Fourier-Eigen embeddings increasingly out-perform the competing methods for $T\in\{100,150\}$, since the estimated embedding vectors concentrate more sharply around the population-level embeddings which are provably well-separated due to Lemma 3.4.

Next, we run K-means clustering on each set of embeddings with $2$ clusters. We evaluate the performance of each method by the adjusted Rand index, using the true class membership indicators as our reference labels. The results for each configuration of $n$, $T$, and signal strength are presented in Figure 2, averaged over 100 replications. Here, we see that our method is able to recover the true class membership labels with greater accuracy than either the code counts or PMI embeddings, especially as the signal strength increases. The Fourier-Eigen embeddings perform poorly when the signal between the classes is too weak, which is consistent with Theorem 3.13. As the signal strength increases, the Fourier-Eigen embeddings recover the true class labels more efficiently than either of the competing methods, with the count-based clustering results aligning with the Fourier-Eigen results in the very high signal setting. The PMI embeddings return lackluster clustering performance across all settings, which we attribute to the restrictive assumption of a log-linear intensity function that PMI-based methods enforce on the code occurrence probabilities.

Finally, we evaluate the robustness of our Fourier-Eigen embeddings to the functional form of the transfer functions $\{\omega_{j\ell}\}$. We set $\omega_{j\ell}(t)=a_{j\ell}\beta(t)$ for $a_{j\ell}\sim\text{Unif}(0,1/2)$ and we vary the transfer kernel function $\beta(t)$ in the set \beta(t)\in\big{\{}\frac{|\sin(t)|_{{\bm{1}}(t<\pi)}}{t+1},(1-\sqrt{t})_{{\bm{1}}(t<1)},(1-t)_{{\bm{1}}(t<1)},4^{-t}_{{\bm{1}}(t<2)}\big{\}}. We simulate data from two groups with signal strength $\|\mu^{(0)}-\mu^{(1)}\|_{2}=1.6$, ambient dimension $d=100$, latent dimension $k=2$, observation time $T=100$, and total sample size $n=500$. For each choice of $\beta(t)$, we compare the efficacy our Fourier-Eigen embeddings to the PMI embeddings and code counts in classification with logistic regression and clustering, again measuring the performance of each method with AUC and the adjusted Rand index. The results are given (averaged over 100 replications) in Table 1.

The Fourier-Eigen embeddings give the best performance for a majority of the choices of $\beta(t)$ in terms of AUC. For the clustering task, our method gives the highest adjusted Rand index in half of the cases, with the count-based clustering achieving the best performance for $\beta(t)=(1-\sqrt{t})_{{\bm{1}}(t<1)}$ and $\beta(t)=(1-t)_{{\bm{1}}(t<1)}$. However, the Fourier-Eigen embeddings perform nearly as well as the counts in both of these settings, resulting in an average adjusted Rand score that is far greater than the competing methods. These results indicate that our Fourier-Eigen embeddings can provide strong statistical perform in downstream tasks across a variety of forms for the conditional intensity function of the observed data, as suggested by our theoretical results.

5 Real data analysis: Alzheimer’s disease severity

We demonstrate the real-world utility of our proposed Fourier-Eigen embedding method on an analysis of electronic health records data from the University of Pittsburgh Medical Center (UPMC). Our objective is to cluster patients with Alzheimer’s disease (AD) according to severity of AD progression using EHR code occurrence data. Alzheimer’s disease is widely regarded to be the most common cause of dementia, which affects over 50 million people worldwide and leads to severe declines in quality of life and increased rates of mortality (Korczyn and Grinberg, 2024). As the rate of progression of AD has consequences for long-term health outcomes (Tahami Monfared et al., 2022), there is a great interest in leveraging the temporal information in EHR patient trajectories to predict progression patterns and identify at-risk patient groups (Fu et al., 2025; Kumar et al., 2021; Akter et al., 2025; Li et al., 2023; Venkatesh et al., 2025). In this analysis, we aim to evaluate the ability of our Fourier-Eigen embeddings to distinguish patients with different AD progression patterns, in comparison to competing methods. We use nursing home admission time as a proxy for AD progression to compensate for a lack of gold-standard labels from chart review.

Our dataset consists of 92,723 electronic health records of patients who have at least one occurrence of an AD diagnosis code (PheCode:290.11), including 48,454 patients who were admitted to a nursing home and 44,269 who were not. We define the observation period for each patient as the two years prior to the date of the first occurrence of an AD-related code. We remove codes that were present in only $<5\%$ of individuals, which leaves us with $d=1,375$ codes, including 545 disease diagnosis codes (PheCode), 80 procedure codes (CCS), 273 prescription codes, 262 laboratory codes (LOINC), and 215 local codes that do not belong to any of these four categories. From here, we draw a random subset of $8000$ patients for the analysis, split evenly between patients with a recorded nursing home admission time and patients who do not. Using each patient’s trajectory, we compute the following three embeddings: 1) Fourier-Eigen embeddings of dimension $k=5$ using a Gaussian kernel with a half-week bandwidth, and frequency $\xi_{0}=1$; 2) PMI embeddings of dimension $k=5$ as defined in Section 4, using a co-occurrence matrix computed with a week-long window size; 3) BERT-based embeddings (Devlin et al., 2019). Specifically, we compute the BERT embeddings by pulling a dimension $m=769$ representation of each code from a pretrained BERT model, giving us the code-representation matrix $\mathbf{C}\in\mathbb{R}^{d\times m}$. Letting $\mathbf{Z}\in\mathbb{R}^{n\times d}$ denote the code count matrix (defined in Section 4), we define the BERT embeddings as the first $k=5$ left singular vectors of the matrix $\mathbf{Z}\mathbf{C}$.

With each set of embeddings as our matrix of features, we run K-means clustering with two clusters to obtain an assigned subgroup label for each patient. To evaluate the ability of these cluster assignments to discern subgroups with differing AD progression, we compute a Kaplan-Meier curve within each cluster, using time in weeks from the first code occurrence to nursing home admission as the outcome. For patients who were not admitted to a nursing home, we define the right-censoring time as the final code occurrence time in their trajectory. The Kaplan-Meier curves with a 95% confidence band and the p-value from a standard log-rank test are given in Figure 3. It is clear that our Fourier-Eigen embedding approach is the only embedding method able to discern two significant groups of patients, correspond to fast and slow Alzheimer’s progression. The BERT-based embeddings produce two clusters with different Kaplan-Meier curves, but the log-rank test is not significant due to high variance. Meanwhile, the PMI embeddings are not able to discern any signal in the data, likely due to severe model misspecification.

We supplement these results with a Cox regression analysis. For each of the three embedding methods, we fit a standard Cox model using time to nursing home admission as the outcome and the clustering membership indicator as the covariate. In Table 2, we provide the estimated hazard ratios, the 95% confidence intervals, and p-values from a likelihood ratio test for each of the three Cox models. The model fit on the clusters produced by the Fourier-Eigen embeddings is the only one able to detect a significant hazard ratio between the subgroups. These results strongly suggest that the Fourier-Eigen embeddings are able to encapsulate interpretable signal between the EHR trajectories from heterogeneous groups of patients, while the PMI and BERT embeddings are clearly less powerful.

Finally, we perform an exploration of the most prevalent EHR features in each of the two clusters derived from our Fourier-Eigen embeddings. In Figure 4, we provide word clouds of the EHR code descriptions, sized by relative within-cluster prevalence, for the slow and fast progression clusters. We can see that the fast progression cluster admits higher rates of chronic airway obstruction and chronic renal failure, both of which have been found to be associated with severe Alzheimer’s disease and dementia (Stanciu et al., 2020; Zhang et al., 2020; Tondo et al., 2018). Meanwhile, the slow progression cluster seems to have a relatively higher number of occurrences of non-chronic renal failure, suggesting that the Fourier-Eigen embeddings identify a cluster of patients with more severe kidney disease, leading to faster Alzheimer’s progression. The fast progression cluster also presents a higher prevalence of mood disorders which is well known to be associated with Alzheimer’s risk (Ownby et al., 2006) and anemias, which may complicate AD progression by increasing risk of cognitive decline (Faux et al., 2014).

6 Conclusion and discussion

In this work, we have described a new model for high-dimensional point processes, called the latent factor point process model, whose structure is inspired by electronic health records data. We have shown that interpretable differences between subgroups of patients can be described with a spectral analysis of the cross-covariance matrix under this model, which we leverage to define low-dimensional patient embeddings called Fourier-Eigen embeddings. We have demonstrated through theoretical results, simulations, and real data analysis that our Fourier-Eigen embeddings capture signal that can be used to build classification models or accurately cluster patients that belong to different subgroups.

Our contributions are not without limitations, and we foresee several promising directions for future research. First of all, we have imposed that observed and latent processes $(\bm{N}_{i},\bm{M}_{i})$ are stationary for all $i\in\{1,\ldots,n\}$. This excludes cases where the intensity of the latent process may change over time, such as when a patient changes from one class to another (for instance, “healthy” to “sick”). Relaxing this assumption would greatly improve the real-world applicability of our Fourier-Eigen embeddings, and would open up new use-cases such as changepoint detection and localization. Additionally, our model only allows for excitatory effects, as we assume that the transfer functions satisfy $\omega_{j\ell}(t)\geq 0$ for all $j,\ell$ and all $t\geq 0$. Allowing for latent processes that decrease the probability of a particular code occurrence is of great interest, and extensions in this direction could build on recent work studying point processes with inhibitory effects (Chen et al., 2017; Gao et al., 2024; Pfaffelhuber et al., 2022; Bonnet et al., 2021; Cai et al., 2024).