Geostatistical Modeling Spatial Statistics with Python

Choosing GWR Bandwidth with Golden-Section Search

TL;DR — Call bw = Sel_BW(coords, y, X, kernel="bisquare", fixed=False).search(criterion="AICc", search_method="golden_section"). The golden-section search minimises the corrected AIC over candidate bandwidths and returns one scalar — a neighbour count (adaptive) or a distance (fixed). Cross-check with criterion="CV", and confirm the minimum with a manual bandwidth sweep. AICc is the default because it resists small-bandwidth overfitting.

Why this matters

The bandwidth is the one parameter that decides whether a geographically weighted regression is a faithful map of local variation or an overfit mess of noise. Too small and every local surface chases sampling noise; too large and GWR collapses into the global model it was meant to improve on. Getting it right is the whole game, and mgwr automates it through Sel_BW. This page is the deep dive on that selection step for the Geographically Weighted Regression in Python guide, within the larger Python workflows for spatial modeling and regression pipeline.

Environment

bash
pip install "mgwr>=2.2" "libpysal>=4.8" "spreg>=1.4" \
            "geopandas>=0.13" "numpy>=1.22" "scipy>=1.9"
python
import numpy as np                 # >= 1.22
import geopandas as gpd            # >= 0.13
from mgwr.sel_bw import Sel_BW     # mgwr >= 2.2
from mgwr.gwr import GWR

How golden-section search works

Sel_BW.search() does not test every possible bandwidth. It uses a golden-section search, a one-dimensional optimiser that repeatedly narrows a bracketing interval [a,b][a, b] by evaluating the criterion at two interior points spaced by the golden ratio φ≈1.618\varphi \approx 1.618. At each iteration it discards the sub-interval that cannot contain the minimum, shrinking the search by a constant factor of 1/φ1/\varphi per step. It converges in O(log⁡(range))O(\log(\text{range})) criterion evaluations instead of scanning a full grid — critical because each evaluation fits a complete GWR model.

Golden-section search narrowing the bandwidth bracket A bandwidth axis running from 30 to 300 nearest neighbours sits above four stacked rows, one per search step. Each row draws the live bracketing interval in indigo, the slice thrown away by that step in red, and everything discarded on earlier steps in grey. Probe bandwidths 133, 197, 236, 172 and 212 are marked with their AICc values, and the bracket shrinks from 30-300 to 172-212 around the selected bandwidth of 197 neighbours. Golden-section search: every fit throws away part of the bandwidth range adaptive bisquare kernel, n = 300 tracts, bounds bw_min = 30 and bw_max = 300 live bracket dropped this step already dropped bandwidth (neighbours) 30 100 150 200 250 300 step 1 133 197 AICc(133) = 2841 > AICc(197) = 2794 → the minimum lies right of 133, so 30–133 is discarded step 2 197 (kept) 236 AICc(236) = 2812 > 2794 → discard 236–300; only one new fit is needed, since 197 is re-used step 3 172 197 (kept) AICc(172) = 2803 > 2794 → discard 133–172; the bracket is now 172–236, one fifth of its original width step 4 197 212 AICc(212) = 2808 → bracket 172–212 after five GWR fits, not a 271-point grid: bw = 197 neighbours

The method assumes the criterion is unimodal over the bandwidth range: a single minimum with the curve rising on both sides. For AICc as a function of bandwidth this generally holds, which is why golden-section is the default search_method. The alternative "interval" method scans a user-defined coarse grid and is a safer choice when you suspect multiple local minima.

Step-by-step

Step 1 — Prepare arrays and configure Sel_BW

python
gdf = gpd.read_file("data/tracts.gpkg").to_crs("EPSG:32633").reset_index(drop=True)
coords = np.array(list(zip(gdf.geometry.centroid.x, gdf.geometry.centroid.y)))
y = gdf["median_price"].values.reshape(-1, 1)
X = gdf[["income", "rooms", "dist_cbd"]].values
X = (X - X.mean(axis=0)) / X.std(axis=0)

# kernel: bisquare (compact) | gaussian | exponential
# fixed=False -> adaptive (neighbour count) ; fixed=True -> distance in metres
selector = Sel_BW(coords, y, X, kernel="bisquare", fixed=False)

Step 2 — Run the golden-section search on AICc

python
bw_aicc = selector.search(criterion="AICc", search_method="golden_section")
print(f"AICc-optimal adaptive bandwidth: {bw_aicc:.0f} neighbours")

criterion accepts "AICc", "AIC", "BIC", or "CV". The returned scalar is the neighbour count for an adaptive kernel or a metric distance for a fixed one.

Step 3 — Cross-check against the CV criterion

python
selector_cv = Sel_BW(coords, y, X, kernel="bisquare", fixed=False)
bw_cv = selector_cv.search(criterion="CV", search_method="golden_section")
print(f"CV-optimal bandwidth: {bw_cv:.0f}   AICc-optimal: {bw_aicc:.0f}")

If AICc and CV land close together, the choice is robust. A large gap — CV often selects a smaller bandwidth — is a signal to inspect the criterion curve directly (Step 5).

Step 4 — Compare fixed vs adaptive under the same kernel

python
# Adaptive: neighbour count expands/contracts with local density
bw_adaptive = Sel_BW(coords, y, X, kernel="bisquare", fixed=False).search(criterion="AICc")

# Fixed: constant radius in metres — only for evenly spaced data
bw_fixed = Sel_BW(coords, y, X, kernel="bisquare", fixed=True).search(criterion="AICc")

res_a = GWR(coords, y, X, bw_adaptive, fixed=False).fit()
res_f = GWR(coords, y, X, bw_fixed,    fixed=True ).fit()
print(f"adaptive AICc={res_a.aicc:.1f}  |  fixed AICc={res_f.aicc:.1f}  (lower wins)")
Fixed versus adaptive bandwidth over the same uneven point pattern Two side-by-side maps show an identical scatter of 38 tract centroids, dense on the left and sparse on the right of each map. In the left map a constant 4 km search radius is drawn around an urban regression point A and a rural regression point B: A encloses 12 tracts, B encloses only 2 and is drawn in red as a failed local fit. In the right map the adaptive kernel shrinks the radius to 3.3 km at A and grows it to 7.2 km at B, so both local regressions run on exactly 10 tracts. One radius cannot serve two densities — which is why fixed=False is the default 38-tract detail of the study area, same bisquare kernel in both maps; only the neighbourhood rule changes in bandwidth outside, weight 0 fixed=True — every point searches the same 4 km radius A B 4 km fixed=False — every point searches its 10 nearest tracts A B 4 km A urban core: 12 tracts inside the 4 km radius — well posed B rural edge: 2 tracts for 3 covariates — the local fit is rank deficient A the radius contracts to 3.3 km to reach 10 tracts B it expands to 7.2 km — still 10 tracts, still estimable

Step 5 — Sensitivity sweep to confirm the minimum

python
from mgwr.gwr import GWR

# Manual grid around the selected bandwidth to see the AICc curve shape
grid = range(int(bw_aicc) - 40, int(bw_aicc) + 60, 10)
for k in grid:
    if k < len(y):
        aicc_k = GWR(coords, y, X, k, kernel="bisquare", fixed=False).fit().aicc
        marker = "  <-- selected" if abs(k - bw_aicc) < 5 else ""
        print(f"bw={k:4d}  AICc={aicc_k:9.2f}{marker}")

A clear, single-basined dip around the selected bandwidth confirms the golden-section search found a well-identified optimum. A flat valley means the result is insensitive over a range — report that range rather than a single number.

Interpreting the selected bandwidth

The bandwidth is not just a tuning knob; it is a result. Read it relative to the sample size nn:

  • Near the maximum (approaching nn, or a very large distance): relationships are effectively stationary. GWR is converging on global OLS, and a global spatial regression model is the more parsimonious choice.
  • Small (a low neighbour count or short radius): strong local non-stationarity. The covariate effects genuinely vary over short distances.
  • Intermediate: the typical, informative case. The value estimates the spatial scale over which the modelled relationships hold roughly constant.
Reading the selected bandwidth as an estimate of spatial scale A horizontal ruler of adaptive bandwidths from 30 to 300 neighbours is divided into an over-local band, an informative band and a near-stationary band. Three markers on the ruler, at 40, 197 and 299, connect down to three maps of the local income coefficient over the same 300 tracts. At bandwidth 40 the map is a high-contrast patchwork with a mean jump of 0.71 between adjacent tracts and AICc 2907; at the selected bandwidth 197 it is a smooth west-to-east gradient with jump 0.14 and AICc 2794; at bandwidth 299 the surface is almost uniform with jump 0.02 and AICc 2831. The bandwidth is a result: it names the scale the relationship holds over local coefficient on income, same 300 tracts; jump = mean change in β between adjacent tracts local β −1 +1 over-local: chasing noise informative: a real spatial scale near-stationary 30 100 200 300 bw = 40 neighbours AICc 2907, jump 0.71 — sign flips tract to tract bw = 197 — selected by AICc AICc 2794, jump 0.14 — a west-to-east gradient bw = 299 = n − 1 AICc 2831, jump 0.02 — prefer global OLS

For an adaptive kernel the bandwidth is a count of neighbours; for a fixed kernel it is a distance you can quote in metres or kilometres as the characteristic scale of the process.

Critical best practices

Prefer AICc as the default criterion

AICc penalises the effective degrees of freedom tr(S)\text{tr}(\mathbf{S}), so it actively resists the over-localised solutions that leave-one-out CV tends to favour. Treat CV as a robustness check, or use it only when out-of-sample prediction accuracy is the explicit objective.

Pass bw_min and bw_max so the search cannot wander into degenerate territory:

python
bw = selector.search(criterion="AICc", search_method="golden_section",
                     bw_min=30, bw_max=len(y) - 1)

An unbounded search can select a neighbour count exceeding nn or a radius that encloses the whole study area, both of which defeat the purpose of local modelling.

Watch for multimodal criterion curves

Golden-section assumes unimodality. If the sensitivity sweep in Step 5 reveals more than one dip, switch to search_method="interval" for a coarse global scan, then refine near the true minimum. Silent convergence to the wrong basin is the main failure mode of the default search.

Re-select on the actual model data, never reuse a stale bandwidth

The optimal bandwidth is specific to the coordinates, the sample size, and the covariate set. Adding a predictor, changing the CRS, or subsampling all shift it. Never carry a bandwidth over from a different run; re-run Sel_BW whenever the data change.

Troubleshooting

Symptom Likely cause Fix
Search returns bw_max Stationary relationships; GWR unwarranted Accept a global model; report stationarity
Search returns the smallest allowed value Overfitting or noisy y; CV criterion too aggressive Switch to AICc; raise bw_min; smooth or clean y
AICc and CV disagree widely Multimodal criterion or weak signal Run the manual sweep; prefer AICc; verify global fit
Search is extremely slow Large n, each evaluation a full GWR fit Bound the search; use bisquare; prototype on a subsample
Bandwidth meaningless (huge distance) Coordinates in degrees, not metres Reproject to a projected CRS before building coords
Selected bandwidth exceeds n (adaptive) No bw_max bound Set bw_max=len(y) - 1

Next steps

With a defensible bandwidth in hand, run the full fit and read the surfaces in how to run geographically weighted regression in mgwr. Return to the Geographically Weighted Regression in Python guide for MGWR, where each covariate gets its own bandwidth.

Frequently Asked Questions

Should I select a GWR bandwidth with AICc or cross-validation?

Use AICc as the default. It penalises the effective degrees of freedom directly, so it resists the small-bandwidth overfitting that leave-one-out cross-validation is prone to. Cross-validation (criterion='CV') is a useful robustness check and is preferable when prediction accuracy, rather than model parsimony, is the explicit goal. When AICc and CV disagree strongly, inspect the AICc curve across a manual bandwidth grid before trusting either.

What does the selected bandwidth value actually mean?

For an adaptive kernel it is the number of nearest neighbours in each local regression; for a fixed kernel it is a radius in the coordinate units, normally metres. A bandwidth close to the sample size means relationships are nearly stationary and GWR approaches global OLS; a small bandwidth means strong local variation. The value is itself an interpretable estimate of the spatial scale of the process.

Why does golden-section search sometimes miss the true optimum?

Golden-section search assumes the criterion is unimodal over the bandwidth range. If the AICc curve has multiple local minima, the search can settle in the wrong basin. Supplying explicit bw_min and bw_max bounds and confirming the result with a manual grid sweep guards against this, as does the interval search method for a coarse global scan.

Related

← Back to Geographically Weighted Regression in Python

Related

How to Run Geographically Weighted Regression in mgwr Multiscale GWR vs Classic GWR