A collapsed Gibbs sampler for Latent Dirichlet Allocation that weights topic-word counts by inverse temporal density, so sparse historical periods contribute equally to learned topics despite having far fewer documents.
atweight.pdf contains a full description of the method and simple validation
tests of standard LDA's failure modes in unevenly sampled regimes. It also shows
how Temporal LDA's weighting method addresses those failures; see Appendices A
and B.
Temporal LDA solves a key problem in the use of topic models on historical archives: uneven sampling. If the number of documents in an archive is increasing as a function of time (for example), a standard topic model will, in optimizing the overall log-Likelihood, tend to make finer distinctions in the later period compared to the earlier: we overmodel (or overfit) the later, more densely-sampled era and undermodel (or underfit) earlier periods.
There are ad hoc solutions, such as simply repeating earlier documents verbatim, but these exact repeats create their own problems -- the model can, for example, try to model these exact repeats using specialized topics.
Temporal LDA solves this problem by weighting the word-by-word counts in localized temporal regions. It also has adaptive prior learning for the document-topic distribution that allows for the possibility that early documents (say) tend to have broader, or narrower, distributions.
In some cases, the time-span of the archive itself causes biases. If you model documents from England in the span 1055 to 1500, you will only have a decade or so of pre-conquest data. Temporal LDA can't correct for this effect! It is possible to create your own weighting files but, for simplicity, our current implementation of the code does not do this. If you have suggestions for updates to Temporal LDA, please let me (sdedeo[at]andrew.cmu.edu) know.
Temporal LDA was developed for projects stemming from the Proofs and Reasons project, https://proofsandreasons.io, and was supported by Grant Number 63750 from the John Templeton Foundation, and by the Survival and Flourishing Fund. This code is released under CC-by-4.0; if you modify or reuse it, it's fine just to mention that fact in your own version in a README.
If you use Temporal LDA in your academic work, feel free to cite it; a good BibTeX entry is:
@misc{temporal_lda,
author = {DeDeo, Simon},
title = {{temporal-lda}},
year = {2026},
howpublished = {\url{https://github.com/simon-dedeo/temporal-lda}},
note = {Accessed: 2026-03-30}
}
The underlying code is written in optimized C for speed; it does no preprocessing of the data (for, e.g., stopwords or anything else) so you will need to write your own processing code to create the necessary input files.
gcc -O3 -std=c11 temporal_lda.c -lm -o temporal_lda
Or use the Makefile:
make
For large corpora, an OpenMP-parallel build runs the Gibbs sweep across all cores of a shared-memory machine:
make omp # builds temporal_lda_omp
OMP_NUM_THREADS=24 ./temporal_lda_omp --docs ... (same CLI as the serial binary)
The parallel sampling path is guarded by #ifdef _OPENMP: compiling without
-fopenmp (that is, plain make) uses the serial sampling path and does not
allocate the thread-local delta arrays.
Design. The sweep is parallelized over documents using a deterministic, synchronous AD-LDA-style update:
n_dk(document–topic counts) are owner-private — each document is touched by exactly one thread — and stay exact.- Global
n_wkandn_kare immutable during a sweep. Each thread samples against the sweep-start snapshot plus its own evolvingn_wk/n_kdeltas. - Thread-local deltas are merged in fixed thread order after the sweep. This avoids concurrent read/write races and makes floating-point merge order repeatable.
- Documents use deterministic static scheduling, so a fixed thread count assigns each document to the same RNG stream on every run.
- Each thread has its own xorshift64* RNG stream, splitmix-decorrelated from
--seed. - The per-year alpha(y) Minka update is parallelized over years (its digamma(alpha_k) terms are hoisted out of the document loop, and the Gaussian kernel is cut at 1e-8 instead of 1e-12 — numerically negligible, and without it the alpha update dominates run time at large D). The log-likelihood and beta updates are parallelized with reductions.
Reproducibility. Runs are byte-for-byte deterministic for a fixed binary,
--seed, thread count, input, and hardware/compiler floating-point behavior.
make test-repro runs the same four-thread fit three times and compares every
model output. Results are not expected to be identical across different thread
counts or between serial and parallel inference, because each thread evolves a
different local view during a sweep.
The deterministic implementation trades memory for correctness: it allocates
threads * vocab_size * K * sizeof(double) bytes for topic-word deltas. For
K=200 and an 86k vocabulary this is about 131 MiB per thread (about 6.7 GiB at
52 threads), in addition to the shared model and corpus arrays. The program
prints the delta allocation at startup; choose a lower thread count when memory
is constrained.
Validation. The historical benchmark numbers below describe the former
atomic/Hogwild implementation and are retained only as provenance. The current
deterministic local-delta implementation must pass make test-repro, the paper
protocol suite, and thread-count equivalence tests before a production binary
is accepted. The calling project records those validation artifacts alongside
its production provenance.
Acceptance results (2026-07-30): three repeated fits were byte-identical at 4, 24, and 52 threads. On a 1,227-document, 1.97M-token Hansard subset (K=50, 200 iterations, three seeds), optimally matched serial-vs-OpenMP mean topic JSD was 0.403--0.487 bits at 1, 24, and 52 threads, below the predeclared 0.573-bit limit (1.25 times the median serial seed-to-seed baseline). The paper's weighted K=4 and K=6 recovery cases and its unweighted K=4 pathology passed. The non-production unweighted K=6 case recovered three semantic groups versus four in the serial reference, so that case is not claimed as an exact recovery pass.
On the same small real-data subset, 24 OpenMP threads took 56.38 seconds and 0.85 GiB peak RSS, versus 78.68 seconds and 0.085 GiB for the serial binary (1.40x wall-time speedup). The deterministic delta arrays are memory-bandwidth bound; high thread counts are primarily for the much larger production fits, not a promise of linear scaling.
A full Commons acceptance benchmark on Orchard (786,047 documents, 558M tokens, K=200, 52 threads) took 101 seconds for input loading, initialization, and ten sweeps, and peaked at 13.6 GiB RSS. For that corpus, production jobs request 16 GiB; a lower memory request needs another measured full-corpus check.
The paper-protocol checks use the shipped test_data/ corpus, eight seeds, and
best-of-eight selection by log-likelihood. Weighted + local-alpha recovered all
6 true topics at K=6 (LL -169,149 versus serial -169,147), and at K=4 both
implementations recovered the intended 2:2 split (LL -196,057 versus serial
-196,056). Unweighted K=4 exactly reproduced the paper's 3:1 pathology at LL
-173,059. Unweighted K=6 did not exactly match the serial rerun: the parallel
fit recovered three unique groups rather than four (LL -170,716 versus serial
-170,514). That non-production comparison is retained as a disclosed limitation,
not counted as an exact recovery pass.
Performance tuning. Two further optimizations are built in: the
per-token probability loop carries an omp simd reduction, and topic
assignments are stored as uint16 rather than int (halving the hottest
per-token memory stream; hence the K ≤ 65535 limit, checked at startup).
For best throughput compile for your machine:
make omp CFLAGS="-O3 -Wall -std=c11 -march=native"
Historical measurement on a 10M-token, K=200 corpus (2 threads, dual Opteron
6274, former atomic implementation):
baseline OpenMP 1.00x -> -march=native 1.15x -> +simd 1.28x -> +uint16
assignments 1.76x. This is a historical optimization ladder, not a scaling
guarantee for the current deterministic sampler. On multi-socket systems,
benchmark numactl --interleave=all at the intended thread count; it can help
large jobs but hurt when all threads and memory fit within one NUMA node.
./temporal_lda \
--docs test_data/documents.txt \
--vocab test_data/vocab.txt \
--metadata test_data/metadata.txt \
--output results/ \
--K 6 --iterations 2000 --seed 99 \
--sigma 10 --optimize-interval 50 --converge 1e-6 \
--local-alpha results/alpha.txt
The --optimize-interval 50 flag is recommended: it learns alpha and beta from the data (the --alpha and --beta values become initial guesses). if you also pass --local-alpha with a filename, it will do a temporal reconstruction of the alpha prior (using the same sigma). Pass --no-weighting to run standard (unweighted) LDA for comparison.
Three plain-text files are required.
One line per document. Each line is space-separated:
year length word_id_1 word_id_2 ... word_id_N
- year — integer publication year of the document.
- length — number of tokens that follow on this line.
- word_id — zero-based integer index into the vocabulary.
Example (3 documents):
1780 5 0 3 7 3 12
1780 4 1 1 5 9
1950 6 2 8 14 14 6 10
One word per line, ordered by ID. Line 0 is word 0, line 1 is word 1, etc.
castle
steam
factory
empire
railway
...
The number of lines must equal the vocab_size declared in metadata.txt.
Key-value pairs describing the corpus dimensions:
num_documents=550
vocab_size=18
year_min=1850
year_max=1950
| Key | Meaning |
|---|---|
num_documents |
Total number of lines in documents.txt |
vocab_size |
Total number of lines in vocab.txt |
year_min |
Earliest year in the corpus |
year_max |
Latest year in the corpus |
All model parameters (K, alpha, beta, sigma, etc.) are set via command-line flags, not in this file.
| Flag | Default | Description |
|---|---|---|
--docs FILE |
(required) | Path to documents.txt |
--vocab FILE |
(required) | Path to vocab.txt |
--metadata FILE |
(required) | Path to metadata.txt |
--output DIR |
(required) | Directory for output files (created if needed) |
--K N |
4 | Number of topics |
--iterations N |
1000 | Maximum number of Gibbs sampling iterations |
--seed N |
current time | Random seed for reproducibility |
--alpha A |
1.0 | Initial document-topic prior per topic (if --optimize-interval is set, this is the starting value for asymmetric Minka) |
--beta B |
1.0 | Symmetric topic-word prior (initial value if --optimize-interval is set) |
--local-alpha FILE |
(off) | Enable per-year asymmetric α(y) and write the K-vector per year to FILE |
--sigma S |
50.0 | Gaussian kernel bandwidth (years) for density estimation |
--converge T |
(off) | Stop early when log-likelihood stabilises (e.g. 1e-6) |
--no-weighting |
(off) | Disable density weighting (run standard LDA) |
--optimize-interval N |
(off; recommended: 50) | Update alpha/beta every N iterations via Minka's fixed-point |
--help, -h |
Print usage and exit |
All outputs are written to the directory specified by --output.
One line per (topic, word) pair:
topic_id word_id probability
For each topic k and word v:
where counts are density-weighted.
One line per (document, topic) pair:
doc_id topic_id probability
For each document d and topic k:
where n_dk counts are unweighted (reflecting actual document composition)
and α_k(y_d) is the per-topic, per-year Dirichlet prior for the document's year.
One line per document:
doc_id year length weight
Documents in sparse periods receive weight > 1; documents in dense periods receive weight < 1. Weights are normalised so the mean across all documents is 1.0.
By default the sampler runs for exactly --iterations iterations. Pass
--converge THRESHOLD to enable adaptive stopping:
./temporal_lda ... --iterations 5000 --converge 1e-6
The per-token predictive log-likelihood is computed every 10 iterations and collected into consecutive windows of 10 values (each spanning 100 iterations). Once two consecutive windows have been filled, convergence is detected when the relative change in window means falls below the threshold:
|mean_curr - mean_prev| / |mean_curr| < threshold
Averaging over windows smooths out per-sample MCMC noise that would defeat
a range-based check. --iterations acts as an upper bound.
Typical thresholds:
1e-3— loose, stops early1e-4— moderate (good default)1e-6— tight, may hit the iteration cap on small corpora
Pass --optimize-interval N to learn alpha and beta from the data using
Minka's fixed-point iteration for asymmetric Dirichlet priors:
./temporal_lda ... --alpha 1.0 --beta 1.0 --optimize-interval 50
Every N iterations (after a 50-iteration burn-in), the sampler updates the K-vector α and scalar β. The asymmetric Minka update for each topic k is:
α_k(y)^new = α_k(y) × Σ_d K(t_d-y) [ψ(n_{d,k}+α_k) - ψ(α_k)]
/ Σ_d K(t_d-y) [ψ(N_d+α_Σ) - ψ(α_Σ)]
where α_Σ = Σ_k α_k is the concentration sum. The shared denominator couples topics through document length; per-topic numerators let inactive topics shrink toward zero.
Without --local-alpha, a single global K-vector is estimated and broadcast
to all years. With --local-alpha FILE, each year gets its own K-vector,
kernel-weighted by temporal proximity.
On the shipped synthetic benchmark, the optimizer typically settles within 100–200 Gibbs iterations from ordinary positive alpha and beta starting values. Other corpora can behave differently, so inspect the optimization trace.
When --local-alpha FILE is passed, the file contains one line per year:
year α_1 α_2 ... α_K
with K columns of per-topic Dirichlet concentrations.
- A Gaussian kernel with bandwidth σ estimates the local word-mass density around each document's year.
- Each document gets weight = 1 / density, normalised to mean 1.
- During Gibbs sampling, document-topic counts (
n_dk) are updated ±1 per token (preserving realistic document composition), while topic-word counts (n_kv,n_k) are updated ±weight (equalising the contribution of sparse and dense periods to learned vocabularies).
A Ruby script generates a synthetic corpus (6 true topics, 1:20 token imbalance between two eras, heterogeneous Dir(0.1)/Dir(0.3) sparsity) and compares weighted+local-alpha vs. unweighted LDA:
ruby gen_test.rb
This writes the corpus to test_data/ (documents.txt, vocab.txt, metadata.txt).
Both conditions use asymmetric Minka optimization (--optimize-interval 50)
and convergence detection (--converge 1e-6). The weighted condition passes
--local-alpha to learn per-year, per-topic α vectors.
See atweight.tex for full methodology and results.