Changelog
Source:NEWS.md
fastconley (development version)
Performance (bit-identical)
A second optimisation round speeds up the pair loop itself, the balanced-panel path, and the fixed cost of every call. Every result is bit-identical to 0.11.1.
- Branch-free pair loop. The fused grid loop screens each row’s candidates without branching (the accept flag advances a write cursor into a bounded batch), evaluates Bartlett weights only for the survivors, and folds them into the row’s accumulator with the running sums held in registers, up to 16 columns per pass. The accept/reject branch at the cutoff boundary mispredicted often enough to dominate the old loop. Pair work is about 1.7-2x faster from k = 8 upward (both kernels) and 1.35-1.65x at k = 1 to 3.
-
Balanced panels. The per-period neighbour-list stream uses the same register-held accumulation. Building the neighbour list no longer evaluates Bartlett weights (
asin/sqrt) just to count neighbours, except within a relative 1e-12 of the cutoff, where a weight can round to zero; the fill pass checks the count. - Cell order. Rows are put in cell order by a stable, parallel radix sort instead of a serial comparison sort: the same order, 0.15 to 0.04 s at one million rows on one thread (0.02 s on four) and 0.93 to 0.16 s at four million.
-
No-op aggregation skipped. With
pixel = 0and no two rows sharing (time, lat, lon), as in scattered cross-sections and panels of distinct unit locations,vcovSpHAC()no longer builds, groups, sorts, and copies a table that merges nothing. It orders the rows as the aggregation would and passes that order to the engine, which composes it into the score gather it performs anyway (the internalFastSpatialMeat()gainedrows =). Shared locations aggregate exactly as before. This is most of the per-call saving for small samples and cross-section loops. - Many regressors. The per-row k x k update of the meat walks the column-major accumulator column by column instead of striding across it. Every element still receives one addition per row in row order, so the result is unchanged bit for bit; the spatial meat is about 1.1x faster at k = 30 and 1.2x at k = 100 (200,000 points, one thread) and unchanged at small k. Engine version 0.11.3.
End to end, vcovSpHAC() on fixest fits with ten regressors and an intercept (five for the small samples), post-estimation only, on a 4-core VM with four threads unless noted; minimum of two interleaved runs, with identical matrices before and after:
| data | cutoff | 0.11.1 | this version | speedup |
|---|---|---|---|---|
| global cross-section, 1,000,000 points, uniform | 100 km | 1.26 s | 0.61 s | 2.05x |
| CONUS cross-section, 100,000 points, uniform, 1 thread | 500 km | 3.67 s | 2.06 s | 1.78x |
| CONUS cross-section, 100,000 points, uniform | 500 km | 1.04 s | 0.57 s | 1.83x |
| CONUS cross-section, 100,000 points, Bartlett, 1 thread | 500 km | 6.94 s | 3.90 s | 1.78x |
| CONUS cross-section, 100,000 points, Bartlett | 500 km | 1.84 s | 1.08 s | 1.71x |
| balanced panel, 50,000 units x 10, uniform, lag 1 | 500 km | 1.78 s | 0.96 s | 1.87x |
| balanced panel, 50,000 units x 10, Bartlett, lag 1 | 500 km | 2.61 s | 1.66 s | 1.57x |
| cross-section, 2,000 points, per call | 100 km | 7.0 ms | 4.0 ms | 1.75x |
| cross-section, 20,000 points, per call | 100 km | 23.5 ms | 12.0 ms | 1.96x |
-
CONLEY_CORE_VERSIONis 0.11.2 because the engine’sspatial_meat()entry point gained the optional score row map (numbers are unchanged). The Stata ados still expect 0.11.1 and load the committed 0.11.1 plugins; the expectation moves with the CI-built 0.11.2 plugins, which bring these engine speedups to Stata (seestata/CHANGELOG.md).
Tests and tooling
- testthat grew by 41 expectations: the unique-location path against the table aggregation (cross-sections and panels with permuted rows), shared and negative-zero keys, the row map against a gathered score matrix on every engine path, and identical covariances from the two aggregation paths end to end.
-
tests/manual/engine-ab-check.shcompares the engine header bit for bit with any git revision on randomized configurations: every kernel and distance, k = 1 to 20, zero, tiny, capped, and antipodal cutoffs, duplicate locations, lattice points exactly at the cutoff, balanced (double and float weights) and band paths, 1 versus 4 threads, and row maps. - This round was checked with the bitwise battery (60 of 60 configurations
identical()to a baseline install), the plugin golden check in strict mode under GCC and Clang, the standalone header check, the edge probes, and the A/B check (2,257 cases under GCC, 1,038 under Clang).
fastconley 0.11.1
CRAN release: 2026-09-26
Two regressions of 0.11.0 in vcovSpHAC.felm(), found by a replication pipeline within days of the release, are fixed. Results for fixest fits and for felm fits with explicit unit/time are unchanged.
-
vcovSpHAC.felm()no longer infersunitandtimefrom the absorbed fixed effects. 0.11.0 used the first two absorbed effects as unit and time whenever both arguments were omitted, which silently restricted spatial pairs to within the second effect’s groups (a region, an ethnic group) and moved Conley standard errors by percentages for anyone relying on the default; every earlier release, and the fixest method, treated such a fit as one cross-sectional block. That behaviour is restored: with both arguments omitted every pair within the cutoff enters and each row is its own unit. Absorbed-effect names remain valid explicitunit/timevalues. Regression test added. Reported from a replication pipeline whose published values were computed with the pre-0.11.0 behaviour. -
vcovSpHAC.felm()no longer takes coordinate columns directly from a data frame it found by name in the caller’s or the formula’s environment: whendata =is not supplied, rows are always recovered through the model-frame alignment (expand.model.felm()), as in versions before 0.11.0. In 0.11.0 a frame with the same name as the one in the felm call and the same row count, for example a per-outcome frame reassigned in a loop or a pooled frame in the global environment, could be used unaligned, moving Conley standard errors silently, or produce NA coordinates that the new validation then rejected. A user-supplieddata =with as many rows as the fit is still taken as aligned. Passdata =explicitly inside helper functions and loops. Regression test added. - Documentation: the
vcovSpHAC()generic’s help page gained a Value section describing the returned matrix (CRAN request).
fastconley 0.11.0
Engine extraction, a Stata port sharing the same C++ engine, and a review round (six independent code reviews of the engine, the R layer, the Stata command, the plugin build, the reghdfe proposal, and the validation coverage) whose findings are fixed below.
Corrections that change numbers
-
Clustered
felmfits got the wrong defaultssc.lfestores the cluster-adjusted degrees of freedom (number of clusters minus one) indf.residualwhen a fit carries a cluster specification, and the small-sample scalen / df.residualinherited that value: on a 20-row fixture the covariance was five times too large. The scale now uses the regression’s own residual degrees of freedom (N - p), so clustered and unclustered fits of the same model give identical Conley covariances. -
vcovSpHAC.felmignored theunitandtimearguments and always used the first two absorbed fixed effects. The names are now matched againstnames(reg$fe)(with the model data as a fallback), a missing name is an error, and the first-two-FE default applies only when nothing was passed. A model with| region + unit + timeused the wrong keys before. -
Non-numeric time collapsed gaps in the serial HAC. Character or factor time values were recoded to consecutive integers, so “2000”, “2002” became lag 1. Values that parse as numbers now keep their numeric spacing (both methods, including the absorbed-FE factor in
felmfits); a time key that does not parse is accepted for spatial blocking only and is an error whenlag_cutoff > 0. -
timewithoutunitis honoured. The fixest method silently dropped the time blocking when no unit was given; it now blocks by time with each row as its own unit, matching the Stata command.lag_cutoff > 0without a unit is an error. -
fixest observation selection uses
fixest::obs(), so models extracted fromsplit =estimations and fits with several stacked selections align their coordinates correctly.fixest_multiobjects,lean = TRUEfeolsfits, and fixed-effects-only models now fail with clear messages instead of misleading ones. - Grid engine Bartlett weights use the same chord form as the pairwise engine (haversine and spherical). Four raster/grid Bartlett battery configurations move by at most 3.3e-14 relative; everything else is bit-identical to 0.10.0 apart from the Bartlett haversine/spherical pairwise change described under “Engine”.
- Edge cases that were silently wrong: exact duplicates at a zero cutoff keep their cross terms, antipodal points with a cutoff above the maximum distance are accepted, Bartlett boundary weights can no longer come out as -2e-16, and a raster lattice at a tiny or zero cutoff keeps its self terms.
New validation and errors
- Non-finite or missing
lat,lon,unit,timevalues, latitudes outside [-90, 90], and longitudes outside [-180, 360] are rejected in R before anything reaches C++; the engine repeats the checks (including scores and cutoff) for its other front-ends. ANAtime used to split co-located observations into separate blocks. - Both
latandlatitude(or any two aliases) present in the data is now an auto-detection error, as documented. - The fixest method accepts
data =either as the full original data (rows addressed throughfixest::obs()) or as a frame already aligned with the fit, such as the fitted model frame aftersubset =or NA removal. A frame with as many rows as the original data is always read as original data (so subsets that only permute rows stay correct); any other row count is an error. - Unknown arguments in
...are errors (typos such aspsd_fxi =were swallowed).maxobsmemwarns once that it is ignored;neighbor = "band"is documented as deprecated but still works. -
lfefits made through the k-class path (anykclass =argument, evenkclass = 1) are rejected: lfe stores the raw endogenous regressors for them instead of the projected 2SLS design, so the sandwich would be silently wrong (akclass = 1fit gave 0.0028 where 2SLS gives 0.0097). Only OLS and ordinary 2SLS are supported. -
ncoresmust be a single finite positive integer. The default remains all logical cores, butgetOption("fastconley.ncores")overrides it and a true-ish_R_CHECK_LIMIT_CORES_(CRAN’s check convention) caps it at 2. - Long computations can be interrupted: the engine polls a front-end hook on the calling thread once per block (
Rcpp::checkUserInterrupt()in R,SF_pollin the Stata plugin) and joins its workers before raising. - Thread-pool construction is exception-safe (a failed
std::threadconstructor used to terminate the process),ncoresis capped at 1024, and the sub-metre-cutoff cell cap no longer collapses all points into one cell. - Raster geometry is validated (
dlat,dlonfinite and positive, finitelat0in range): a zero longitude step used to spin forever.
Engine
- Haversine weights are computed from the cached 3D unit vectors (
sin^2(theta/2) = |u_i - u_j|^2 / 4, distance2R asin(sqrt a)), and the spherical Bartlett weight uses the same asin form instead ofacos(dot), which is ill-conditioned at small angles. No per-pairsin/atan2calls remain, which speeds up the pairwise engine everywhere (about 2x on Linux, 6x on mingw-built Windows binaries, whose libm implements those functions in x87 microcode). Results move at the 1e-14 relative level versus 0.10.0 on Bartlett haversine/spherical configurations. - The C++ engine lives in a front-end-agnostic header,
src/conley_core.h(namespaceconley,CONLEY_CORE_VERSION0.11.1, bumped on every change that alters results or entry points), shared with the Stata plugin. It compiles without R against plain Armadillo under C++14 (tests/manual/core_standalone_check.cpp+.Rverify bit-identity). - RcppParallel dependency dropped: threading is a
std::threadpool inside the engine with the same deterministic chunked reduction, soncoresstill never affects results.Depends: R (>= 4.0)declared for the posix-threads Rtools toolchain on Windows. - Dead code from the pre-chord screens removed; the balanced path builds its coordinate cache from the first period only.
Performance (bit-identical)
Three optimisation reviews (2026-09-05) led to preparation and engine changes that leave every result bit-identical (the 60-configuration bitwise battery, the standalone header check, and the plugin golden check all pass unchanged):
- Balanced-panel coordinate validation compares each period’s coordinates with period 1 instead of grouping by unit twice; the 10,000-unit panel call drops from 0.29 s to 0.07 s and a 500k-row panel from 3.9 s to 2.1 s.
- The fitted design stays a matrix (no data.table round trip), felm’s bread is computed from the fit’s design, and the serial HAC reuses the raw scores through a row map instead of re-sorting the wide table; allocations fall by a third at one million rows.
- Raster detection rejects scattered data on latitude before sorting longitude.
- Engine: the raster Bartlett engine batches its weight FFTs so KissFFT’s worker and twiddle table are shared across ring pairs; the raster window prepass and prefix build parallelise per ring; neighbour-range construction uses monotone cursor sweeps instead of ten binary searches per cell.
Tests and tooling
- testthat grew from 48 to 106 expectations: independent dense references for
fepois/feglm, felm and fixest IV, the serial HAC, and pixel aggregation; 1-vs-2-coreidentical()checks for every engine path; and regression tests for each item above. -
tests/manual/bitwise-battery.Rgainedwrite-ref/checkmodes with a committed reference (bitwise-battery-reference.txt), andtests/manual/core_edge_probes.cppexercises the edge cases in C++. - A GitHub Actions
R CMD check --as-cranworkflow runs on Linux, macOS, and Windows. - Documentation reconciled with the code: Conley (1999) throughout,
method = "auto"as the default, chord versus great-circle wording, the balanced-panel error (not a fallback), and the pixel-aggregation distance bound (sqrt(2) * pixelfor a pair, notpixel / 2).
Stata port
-
stata/holdsfastconley, a reghdfe-family command with a pure-Mata engine and a compiled plugin built fromsrc/conley_core.hfor Linux, Windows, and macOS (seestata/README.mdandstata/CHANGELOG.md), andstata/upstream/two alternative reghdfe proposals: a genericvce(external PROVIDER, ...)hook with fastconley as its first provider (primary), and a self-contained nativevce(conley ...)patch.
CRAN preparation.
-
DESCRIPTION: package names quoted, Conley (1999) DOI added, copyright holder role declared;LICENSEupdated to match. -
.Rbuildignorekeeps development-only files (the agent instructions,notes/,tests/manual/,stata/,.github/, andcran-comments.md) out of the source tarball. - Tests: the tiny
felmfixture in the balanced-panel validation test no longer emits “NaNs produced” warnings.
fastconley 0.10.0
GLM support and documented IV support.
fixest::feglm() / fixest::fepois() fits supported
vcovSpHAC.fixest now accepts GLM fits (any feglm family, including fepois). The variance is the M-estimation sandwich H^{-1} B H^{-1}, built from the maximum-likelihood score matrix and inverse Hessian that fixest stores on every (non-lean) fit — the same construction fixest’s own vcov_conley() uses for GLMs, verified against it at the distance-formulation tolerance and against an exact same-distance yardstick at ~1e-15. No estimation flag is needed (demeaned = TRUE is only required for feols); weights, offsets, and the fixed-effect profiling are already folded into the stored scores. Because the scores ride through the existing engines unchanged, everything composes: pairwise and grid/FFT engines, pixel aggregation, ssc, psd_fix, and — beyond what fixest offers — the panel spatial + serial HAC via lag_cutoff, now available for Poisson/GLM panels. lean = TRUE fits (no stored scores) and femlm()/feNmlm() fits are rejected with clear errors.
IV/2SLS support documented and tested
IV fits have in fact always produced the correct 2SLS Conley sandwich through both methods — lfe and fixest store the projected (second-stage) design in cX / X_demeaned and the structural residuals in residuals, which is exactly what the score construction needs. This is now documented and covered by the validation suite: felm IV and feols IV agree with each other and with the exact yardstick at ~1e-15, including weighted IV, multiple endogenous regressors, and IV panels with serial HAC. See tests/manual/test-glm-iv-parity.R.
fastconley 0.9.0
fixest feature parity: weighted fits, small-sample correction, PSD repair, lat/lon auto-detection. Breaking: ssc and psd_fix default to TRUE to match fixest’s defaults out of the box; pass ssc = FALSE, psd_fix = FALSE to reproduce earlier fastconley versions and rbluhm/conley bit-for-bit.
Weighted (WLS) fits supported
vcovSpHAC now accepts weighted felm() and feols() fits. The meat scores become s_i = w_i * e_i * x_i and the bread (X'WX)^{-1} — the formula fixest’s own weighted Conley vcov uses (verified exactly against it with a self-pairs-only cutoff, rel. err ~1e-15). Weights enter only the scores and the bread, so every engine — pairwise, grid/FFT, serial HAC, pixel aggregation, balanced CSR reuse — works unchanged. Note lfe stores sqrt(w) on the fit; vcovSpHAC squares it back.
ssc: small-sample correction (default TRUE)
Scales the variance matrix by n / (n - K), with K counting all estimated parameters including absorbed fixed-effect levels (taken from the fit’s residual degrees of freedom). This is exactly fixest’s default Conley correction — its cluster adjustment (G.adj / cluster.adj) is a no-op for Conley vcovs, so this one factor reproduces fixest defaults. ssc = FALSE applies no correction, matching rbluhm/conley and previous fastconley versions.
psd_fix: positive semi-definite repair (default TRUE)
Conley spatial kernels do not guarantee a PSD variance matrix. With psd_fix = TRUE (default, as in fixest) negative eigenvalues are clamped to 1e-16 — the same semantics as fixest’s vcov_fix — with a warning when the fix noticeably changed the matrix (> 1e-8). With psd_fix = FALSE the matrix is returned as computed and a warning is emitted when it is noticeably non-PSD.
lat/lon auto-detection
lat and lon now default to NULL and are auto-detected from the data’s column names (lat/latitude and lon/long/longitude/lng, case-insensitive exact matches). A message reports the pick; ambiguous or missing matches error with instructions.
Validation
- 24-config battery (balanced/unbalanced/cross-section x kernels x distances x lags x engines) bitwise-identical to v0.8.0 with
ssc = FALSE, psd_fix = FALSE(the C++ engines and the default-off code paths are untouched). - Weighted fits: exact match to
fixestat self-pairs-only cutoff; ~1e-10 against a dense brute-force reference at 200 km for both kernels; felm and fixest paths agree to ~1e-10. - New testthat coverage:
tests/testthat/test-weights-ssc.R(14 tests).
fastconley 0.8.0
Grid engine, part two: bartlett support (ring-FFT) and dateline wrap.
method = "grid" / "auto" now covers kernel = "bartlett"
On a lattice the bartlett weight varies with the longitude offset, so the per-ring-pair inner sum is a true 1D convolution rather than a boxcar. FastGridMeat computes it via FFT (arma::fft): per ring pair, the even-symmetric weight vector’s (real) spectrum multiplies cached per-ring score spectra, with one inverse FFT per target ring. Score spectra are cached for a sliding latitude band plus a cutoff halo, and the reduction is deterministically chunked — results remain bit-identical across ncores.
Weights use the same per-distance arithmetic as the pairwise engine (atan2 haversine, acos spherical, sqrt chord), and the same-cell distance is hard-set to 0, so agreement with the pairwise engine is ~1e-15 (haversine) to ~1e-12 (spherical/chord — inherent conditioning of acos/sqrt near zero distance, not algorithm error). The "auto" rule uses an FFT-aware cost model for bartlett.
- C1-study config (2.25M cells at ~1.1 km, 250 km cutoff, ~1.8e11 pairs, bartlett/spherical): 2.9 s vs 710 s for the 16-core pairwise engine (242x); single-threaded (20.5 s) it still beats 16-core pairwise 35x. The kernel’s evenness halves the work (T(r2,r1) = T(r1,r2)’, so only upper-triangle ring pairs are computed, with self-ring weights halved and the half-meat symmetrized).
Dateline wrap (bug fix for global rasters)
v0.7.0’s grid engine clamped longitude windows at the lattice edges, so on a raster spanning the full 360° circle it silently missed pairs that are close “the short way” across the dateline (observed ~5e-3 relative error on a global test raster). The engine now detects when the accept window reaches across the dateline gap and switches both kernels to circular windows (modular prefix-sum arcs for uniform; circular convolution with period n_col_full for bartlett) — exact, validated against the pairwise engine. When wrap would be needed but the lon step does not tile 360° evenly (no consistent circular lattice exists), method = "grid" stops with an informative error and method = "auto" falls back to the pairwise engine. Non-wrapping rasters are unaffected: results are bitwise identical to v0.7.0 (verified on the 30-config battery).
fastconley 0.7.0
Workstream C2: exact grid-native meat for raster data.
New: method = c("auto", "pairwise", "grid")
For the uniform kernel on a regular lat/lon lattice (raster cell centers, gridded covariates), the within-cutoff accept set between two latitude rings is a longitude-index interval, so the spatial meat reduces to sliding-window sums over per-ring prefix sums — FastGridMeat. Cost is O(n_ring * window * n_col * k), independent of the pair count, and the accept threshold is the same dot-product constant the pairwise engine uses, so the result is exact (agrees to FP summation order; no approximation anywhere for natively gridded data).
- C1-study config (2.25M cells at ~1.1 km, 250 km cutoff, ~1.8e11 pairs): 0.81 s vs 244 s for the 16-core pairwise engine — 302x, error 5e-15.
- 9M-cell continental raster (3000 x 3000, 250 km, ~7e11 pairs): 3.5 s.
- Supports all three
dist_fns (monotone in the chord), multiple time blocks (panels apply it per period), sparse occupancy, duplicate cells, and is deterministic acrossncores.
method = "auto" (the new default) switches to the grid engine only when a lattice is detected, the kernel is uniform, and a flop-balance estimate says it wins; otherwise the pairwise engine runs as before. Scattered (non-lattice) data is unaffected. method = "grid" errors informatively when its requirements are not met.
The bartlett ring-FFT variant (exact per the C1 study) is planned as a follow-up; bartlett rasters currently stay on the pairwise engine.
fastconley 0.6.1
Minor-backlog items M1-M4 from notes/OPTIMIZATION_PLAN.md.
- Serial HAC is now O(T * k) per unit instead of O(T^2 * k): the Bartlett-in-time weight decomposes over a sliding two-pointer window (with a per-block time shift for conditioning). A 100k-row panel with T = 2000 and lag = 50 computes in ~6 ms. Rows must be time-sorted within unit blocks (the R layer always sorts; direct callers get a clear error).
-
vcovSpHAC.felmgainsdata =: when the passed frame has the same row count as the fit (no NAs dropped, no subset), coordinates are taken by direct column access with no model-frame re-evaluation; otherwise the frame overrides the call-recovered data in the aligned model-frame path. -
Prep path parallelized (grid sort via TBB parallel_sort, coordinate cache trig, permutation gathers) — all order-preserving, so results stay bit-identical across
ncores. 1M-row global cross-section at 16 threads: 0.46 -> 0.21 s (scaling 2.2x -> 4.8x). - Spatial results are bitwise identical to v0.6.0; the serial meat differs by ~1e-15 relative from the new summation algebra.
- README benchmarks refreshed against fixest 0.14.1 at scale: 10.7-45.8x on cross-sections up to n = 1M, 1.4-2.9x on panel SHAC up to 1.5M obs (where fastconley computes ~T x more spatial work by construction). A brute-force arbitration shows fastconley matches the exact great-circle uniform meat to ~1e-15 while fixest 0.14.1’s “spherical” distance deviates ~2e-2.
fastconley 0.6.0
Phase 2 of notes/OPTIMIZATION_PLAN.md: memory diet, deterministic reduction, and screen work. Results remain exact (same pairs, same weights); summation order changed, so values differ from v0.5.0 by <= ~5e-15 relative.
Results are now invariant to ncores
All meat accumulations use a deterministic chunked reduction (fixed-size row/block chunks, partials summed in chunk order). ncores = 1 and ncores = 16 produce bit-identical matrices — the old multicore tolerance caveat is gone.
Memory
- The C++ entry points take a pre-computed score matrix and alias R memory directly (no Rcpp input copies, no intermediate unsorted score buffer; one gather pass builds the sorted layout).
vcovSpHACpasses the aggregated scores straight through;X/earguments remain on the internal wrappers for convenience. Measured peak RSS: 4M-row cross-section 2321 -> 1558 MB (-33%); 40k x 40 balanced panel 1028 -> 748 MB (-27%). Same speed, identical checksums. - New
csr_weight = c("double", "float"): opt-in float storage for the balanced-path bartlett weights (8 -> 4 bytes per pair, <= ~6e-8 relative error per weight). Default stays exact double.
Speed
- Haversine gains a conservative 3D dot-product pre-screen (exact identity
a = (1 - dot)/2); the exact a-test remains the arbiter, so results are bit-identical. CONUS 100k / 500 km bartlett/haversine: 15.1 -> 10.5 s single-threaded (1.4x). - Balanced 16-thread uniform meat: 2.30 -> 1.56 s on a 40k x 40 panel (better task granularity from the chunked reduction).
Tried and reverted (documented for the record)
A unit-major T*k stacked score layout that streams the balanced CSR once instead of once per period was implemented and benchmarked. It lost to the period-major layout: spatial sorting makes neighbor gathers a sliding window of k-wide rows that stays L2-resident per period; any wider stacking pushes the window past L2 and thrashes the shared L3 at high thread counts (2.6 s vs 2.0 s at 16 cores). The period-major traversal stays, now deterministic and float-capable.
fastconley 0.5.0
Phase 1 of notes/OPTIMIZATION_PLAN.md: 3D cell-grid neighbor search.
New: neighbor = c("grid", "band")
The spatial meat’s candidate enumeration now defaults to a 3D cell grid. Points are bucketed by their unit vectors into a cubic grid whose edge is the unit-sphere chord equivalent of the cutoff; every supported distance is monotone in the chord, so accepted pairs are never more than one cell apart per axis — no pole or dateline special cases. Each row scans its own cell plus five contiguous row ranges covering the 13 forward neighbor cells: ~3–4 candidates per accepted pair, independent of geographic extent, versus the latitude band scan’s 2*L_lon/(pi*r).
Both strategies call the identical per-pair accept test, so pair sets and weights are exactly the same; results differ only by floating-point summation order (observed <= 6e-15 relative across the validation matrix; neighbor = "band" remains bitwise identical to v0.4.1). The band path is kept for one release and will be removed in v0.6.0.
Measured single-threaded speedups (FastSpatialMeat, k = 10):
- global extent, n = 1M, 100 km: 8.6x (uniform/spherical) to 25.3x (bartlett/haversine)
- CONUS, n = 100k, 100 km: 4.5x; 500 km: 2.1-2.3x (that config is dominated by true-pair accumulation, not candidate search)
- balanced CSR build, 200k units x 4 periods, 100 km: 3.6x
- 16-thread scaling on CONUS 100k/500 km improves from ~5.3x to ~7.9x (less memory-bandwidth pressure from candidate streaming)
fastconley 0.4.1
Phase 0 quick wins from notes/OPTIMIZATION_PLAN.md (Q1–Q6). Numerical results are unchanged (verified bitwise against v0.4.0 on the balanced / general / unbalanced × kernel × distance matrix at ncores = 1).
Memory
- The balanced-path CSR no longer allocates its weight array for
kernel = "uniform"(it was filled but never read), and column indices are now 32-bit. Per stored pair: 16 bytes -> 4 (uniform) / 12 (bartlett). A guard errors if a single period exceeds 2^32 - 1 units.
fastconley 0.3.0.9000
Development build for testing a faster spatial HAC path.
Spatial meat changes
- Adds a
FastSpatialMeat()internal Rcpp routine. - Replaces dense
n x nspatial distance matrices invcovSpHAC()with a screened spatial edge list. - Computes the meat via cumulative scores:
S = e * X,C_i = 0.5 S_i + sum_j w_ij S_j, andS'C + C'S. - Supports both
kernel = "bartlett"andkernel = "uniform". - Supports
dist_fnchoices: haversine, spherical, chord. The upstreamflatearthoption was dropped — its formula was not symmetric in(i, j)and the equirectangular approximation is no longer worth a separate code path now that spherical+uniform runs as a 3D dot-product threshold with no trig. - In balanced panels, builds one spatial edge list from the first period and reuses it for all periods.
Performance work (most recent)
-
Templated
(dist_fn, kernel)dispatch. The screen + distance + kernel-weight stack is specialised at compile time for each of the six(distance, kernel)combinations, removing per-pair runtime branches. -
3D unit-vector threshold for spherical and chord.
(spherical, uniform)pair inclusion now reduces to a 3D dot product compared againstcos(cutoff / R)— no trig in the inner loop.(chord, uniform)uses squared-Euclidean against(cutoff / R)². Thebartlettvariants only payacos/sqrton accepted pairs. -
Row-major score buffer. Scores
s_i = e_i · X_iare stored row-major in a single flatstd::vector<double>instead of being recomputed against column-major Armadillo memory on every pair. - Sorted-order layout. Both the coordinate cache and the scores buffer are permuted into latitude-sorted order before the meat workers run, so every read in the hot loop is sequential — independent of how the caller laid out the input.
-
pixelargument for score pre-aggregation (R/vcovSpHAC.R). Atpixel = 0rows that share(lat, lon)within a time period are collapsed exactly; atpixel > 0rows are first snapped to a uniformpixel-km grid (speed/accuracy trade-off). - Cross-section vs
fixest::vcov_conley(kernel = "uniform", spherical, 500 km cutoff):fastconleyis now 3-4× faster atpixel = 0and 27-62× faster atpixel = 25acrossn ∈ {5k … 100k}, and scales better than fixest past 4 threads.
Notes
- This is meant as a test branch/source drop, not a CRAN-ready release.
- The implementation uses RcppParallel for CSR construction and meat accumulation.
- Rcpp exports were updated manually. Running
Rcpp::compileAttributes()is recommended after further edits.