A stubborn problem in geodesy—the science of measuring Earth’s shape and gravity—may have received a major computational upgrade. Researchers have introduced a Nyström-based version of least-squares collocation, a powerful mathematical technique used to model the geoid, the irregular surface that best represents global mean sea level beneath the continents. Their approach allows regional geoid models to be calculated from hundreds of thousands of measurements without confronting the enormous computational burden that has traditionally limited the method. In a demonstration using the well-known Auvergne geoid computation test dataset in France, the researchers processed 243,889 measurements and produced a geoid model with a standard deviation of 3.0 centimeters when checked against 75 GNSS/leveling benchmarks. The result suggests that a technique long regarded as computationally expensive could become practical for increasingly dense gravity and positioning datasets.
The geoid is not a perfectly smooth sphere. Variations in underground rock density, mountains, ocean trenches, and other geological structures subtly alter Earth’s gravitational field, causing the level surface used to define heights to rise and fall by many meters relative to a reference ellipsoid. Accurate geoid models are essential because satellite positioning systems such as GPS determine heights relative to a mathematical ellipsoid, while engineering, surveying, flood-risk analysis, and infrastructure construction generally require heights related to mean sea level. Converting between these systems depends on knowing the geoid. In regional modeling, scientists combine gravity anomalies, terrain information, satellite-derived geopotential models, and precise height observations. The challenge is to estimate the gravitational signal between measured locations while separating it from noise and measurement errors.
Least-squares collocation, or LSC, is one of the most flexible tools for this task. It treats the unknown geophysical field as a correlated signal and uses covariance functions to describe how measurements at different locations are statistically related. In simplified terms, nearby observations usually provide more information about one another than distant observations, although the exact relationship depends on the assumed gravity-field model. To predict the geoid at a new location, LSC compares that point with every observation through a covariance vector and applies weights obtained from the full covariance matrix. This framework can incorporate different types of observations and account explicitly for their uncertainties, making it attractive for regional geoid determination. But its strength has traditionally been tied to its greatest weakness: the need to manipulate a matrix whose dimensions equal the number of observations.
For a dataset containing n measurements, the covariance matrix has n by n elements. Constructing and directly inverting such a matrix can require computational work that increases approximately as the cube of n, or O(n³). Memory requirements also grow as O(n²), quickly becoming prohibitive. With nearly a quarter of a million observations, a conventional implementation would need to store and process a matrix containing tens of billions of entries. Even when modern computers can technically handle portions of that workload, the time, memory, and data-transfer demands can make the calculation impractical. The new study addresses this bottleneck by approximating the large covariance matrix with a lower-rank representation. Instead of preserving every possible interaction independently, the method seeks to capture the dominant structure of the matrix using a much smaller set of representative observations.
The strategy is based on the Nyström method, an idea originally developed for the numerical solution of integral equations and later adopted in large-scale machine learning. In the geoid application, a subset of p observations is sampled from the full dataset, where p is much smaller than n. Covariances involving these sampled points are used to construct a low-rank approximation of the original matrix. The resulting calculation replaces the conventional O(n³) cost with approximately O(p²n), while reducing the main memory burden associated with the full matrix. The approximation is not simply a crude shortcut: it is designed to preserve the most important covariance patterns governing the predicted field. In the Auvergne experiment, the researchers selected 24,223 sampled points—roughly 10 percent of the 243,889 available measurements—and used them to represent the much larger system.
A second innovation was necessary because even a low-rank method can encounter a practical obstacle when the intermediate matrices are too large to form in memory. The researchers therefore developed a blockwise algorithm. Rather than processing every observation simultaneously, the data are divided into manageable blocks. Each block contributes to the required matrix products or summations, and the partial results are accumulated, potentially in parallel. This approach changes the problem from “fit the entire dataset into memory” to “process a sequence of smaller pieces while preserving the same global calculation.” Blockwise parallelization is particularly important for large geospatial datasets, which are often assembled from irregularly distributed measurements and may contain different observation types, coordinate systems, and uncertainty levels. It also creates a route toward implementation on high-performance computing systems without requiring a single machine to hold the complete covariance structure.
The Auvergne test was designed to reveal whether the computational savings would come at the expense of geoid accuracy. After data curation, all 243,889 measurements were included in the experiment, and the resulting solution was evaluated on a 2-arc-minute by 2-arc-minute grid. The model was independently checked against 75 GNSS/leveling benchmarks, which provide a demanding test because they connect satellite-derived positions with conventional height differences. The reported standard deviation was 3.0 centimeters. That level of agreement places the low-rank solution within the range needed for many regional geodetic applications and indicates that retaining the dominant covariance structure can be sufficient for high-quality prediction. The result is especially notable because the method used only about one-tenth as many sampled points as total observations in its reduced representation.
The researchers also compared their model with four existing geoid models for the Auvergne region. Across the three models considered most consistent, the root-mean-square differences ranged from 2.83 to 3.28 centimeters. RMS values summarize the typical magnitude of discrepancies and are useful for identifying whether competing models describe broadly the same regional gravity signal. No geoid model is expected to match another exactly: differences can arise from the gravity data used, terrain corrections, reference systems, covariance assumptions, filtering choices, and treatment of long-wavelength signals. The close range reported in the comparison suggests that the Nyström-LSC model behaves like a conventional regional solution rather than producing a mathematically convenient but geophysically distorted surface.
Controlled simulations provided another test of the method’s reliability. The researchers varied data density and measurement noise to examine whether the low-rank approximation became unstable under less favorable conditions. According to the study, Nyström-LSC maintained nearly identical prediction accuracy and noise robustness to standard LSC across the tested scenarios. This matters because real-world gravity datasets are rarely uniform. Some areas may have dense terrestrial measurements, while others depend on airborne or satellite observations; noise levels can also differ substantially between instruments and observation campaigns. A method that performs well only under ideal sampling would have limited value. The simulation results instead indicate that, when the sampled points adequately represent the dataset’s spatial and statistical structure, the approximation can preserve the main benefits of full collocation while sharply reducing computational demands.
The advance could have implications beyond the Auvergne test region. Gravity observations are becoming more abundant as airborne surveys, satellite missions, mobile sensors, and national geodetic programs generate increasingly detailed measurements. Higher resolution can improve models, but it also intensifies the matrix problem at the heart of LSC. A scalable low-rank formulation could make it easier to update regional geoids, combine new observations with legacy datasets, and support near-real-time or repeated modeling in areas affected by infrastructure development, groundwater change, or vertical land motion. The method does not eliminate the need for careful covariance design, representative sampling, data quality control, or validation against independent benchmarks. Its success depends on how well the reduced matrix captures the original system. Nevertheless, by pairing Nyström approximation with blockwise parallel computation, the researchers have transformed a classical geodetic method into a more practical tool for the era of big spatial data. The study’s central message is clear: the future of centimeter-level Earth measurement may depend not only on collecting more observations, but also on finding smarter ways to calculate with them.
Subject of Research: Regional geoid modeling and computational geodesy
Article Title: Nyström least-squares collocation for regional geoid modeling with a large dataset
Article References: Zhang, X., Chang, G., Bian, H. et al. “Nyström least-squares collocation for regional geoid modeling with a large dataset.” Earth Science Informatics, volume 19, article 164 (2026).
Image Credits: AI Generated
DOI: 10.1007/s12145-026-02217-4
Keywords: Geoid modeling; least-squares collocation; Nyström method; Auvergne; geodesy
