ScalableGWR

This page documents 1 public symbols. Each entry includes its purpose, import path, full API docstring, and the maintained example that exercises it.

Conceptual guide

ScalableGWR

Scalable GWR using a linear multiscale polynomial kernel.

Property	Value
Type	class
Import	from pygwrx.models import ScalableGWR
Signature	ScalableGWR(bandwidth: 'int' = 100, kernel: 'str' = 'gaussian', polynomial: 'int' = 4, criterion: 'str' = 'cv', optimize_bandwidth: 'bool' = True, scale: 'Optional[float]' = None, penalty: 'Optional[float]' = None, fit_intercept: 'bool' = True, sample_size: 'Optional[int]' = None, random_state: 'Optional[int]' = None, optimizer_maxiter: 'int' = 200, numerical_jitter: 'float' = 0.0, verbose: 'bool' = False) -> 'None'
Maintained example	examples/models/12_scalable_gwr.py

ScalableGWR

ScalableGWR(
    bandwidth: int = 100,
    kernel: str = "gaussian",
    polynomial: int = 4,
    criterion: str = "cv",
    optimize_bandwidth: bool = True,
    scale: Optional[float] = None,
    penalty: Optional[float] = None,
    fit_intercept: bool = True,
    sample_size: Optional[int] = None,
    random_state: Optional[int] = None,
    optimizer_maxiter: int = 200,
    numerical_jitter: float = 0.0,
    verbose: bool = False,
)

Scalable GWR using a linear multiscale polynomial kernel.

ScaGWR approximates a continuous Gaussian or exponential kernel by a weighted sum of polynomially transformed base kernels. Only bandwidth nearest neighbours are used for local cross-products, while penalty adds a global OLS cross-product term. The large local cross-products are computed once, before optimization, so calibration remains linear in the number of observations for fixed feature count, polynomial degree, and neighbour count.

Parameters:

Name	Type	Description	Default
bandwidth	int	Number of nearest neighbours, denoted Q in the paper.	100
kernel	str	Base kernel, either "gaussian" or "exponential".	'gaussian'
polynomial	int	Polynomial degree used to approximate the kernel.	4
criterion	str	Parameter-calibration criterion, "cv" or "aicc".	'cv'
optimize_bandwidth	bool	Optimize scale and penalty parameters. The neighbour count itself is fixed in ScaGWR and is not optimized.	True
scale	Optional[float]	Fixed positive scale parameter when optimization is disabled, or optional initial value when optimization is enabled.	None
penalty	Optional[float]	Fixed non-negative global shrinkage parameter when optimization is disabled, or optional initial value when optimization is enabled.	None
fit_intercept	bool	Add a spatially varying intercept.	True
sample_size	Optional[int]	Optional number of target sites used during CV calibration. All observations remain available as neighbours and in the global shrinkage term. Ignored for AICc calibration.	None
random_state	Optional[int]	Random seed used when sample_size is specified.	None
optimizer_maxiter	int	Maximum L-BFGS-B iterations.	200
numerical_jitter	float	Explicit diagonal stabilization added to every local system after the published global penalty term.	0.0
verbose	bool	Print calibration information.	False

Notes

The historical pyGWRx class with this name was a kNN-truncated ordinary GWR that still formed a full distance matrix. It was not ScaGWR. This class implements the published polynomial-kernel estimator instead.

Source code in src/pygwrx/models/scalable_gwr.py

def __init__(
    self,
    bandwidth: int = 100,
    kernel: str = "gaussian",
    polynomial: int = 4,
    criterion: str = "cv",
    optimize_bandwidth: bool = True,
    scale: Optional[float] = None,
    penalty: Optional[float] = None,
    fit_intercept: bool = True,
    sample_size: Optional[int] = None,
    random_state: Optional[int] = None,
    optimizer_maxiter: int = 200,
    numerical_jitter: float = 0.0,
    verbose: bool = False,
) -> None:
    if not isinstance(bandwidth, (int, np.integer)) or int(bandwidth) < 2:
        raise ValueError("bandwidth must be an integer neighbour count >= 2.")
    kernel_key = str(kernel).strip().lower()
    aliases = {"gau": "gaussian", "exp": "exponential"}
    kernel_key = aliases.get(kernel_key, kernel_key)
    if kernel_key not in self._SUPPORTED_KERNELS:
        raise ValueError(
            "ScalableGWR supports only continuous Gaussian and exponential kernels."
        )
    criterion_key = str(criterion).strip().lower()
    if criterion_key not in self._SUPPORTED_CRITERIA:
        raise ValueError("criterion must be 'cv' or 'aicc'.")
    if not isinstance(polynomial, (int, np.integer)) or int(polynomial) < 1:
        raise ValueError("polynomial must be a positive integer.")
    if scale is not None and (not np.isfinite(scale) or scale <= 0):
        raise ValueError("scale must be finite and positive.")
    if penalty is not None and (not np.isfinite(penalty) or penalty < 0):
        raise ValueError("penalty must be finite and non-negative.")
    if sample_size is not None and (
        not isinstance(sample_size, (int, np.integer)) or int(sample_size) < 2
    ):
        raise ValueError("sample_size must be an integer >= 2.")
    if (
        not isinstance(optimizer_maxiter, (int, np.integer))
        or optimizer_maxiter < 1
    ):
        raise ValueError("optimizer_maxiter must be a positive integer.")
    if not np.isfinite(numerical_jitter) or numerical_jitter < 0:
        raise ValueError("numerical_jitter must be finite and non-negative.")

    self.bandwidth = int(bandwidth)
    self.adaptive = True
    self.kernel = kernel_key
    self.polynomial = int(polynomial)
    self.criterion = criterion_key
    self.optimize_bandwidth = bool(optimize_bandwidth)
    self.scale = scale
    self.penalty = penalty
    self.fit_intercept = bool(fit_intercept)
    self.sample_size = int(sample_size) if sample_size is not None else None
    self.random_state = random_state
    self.optimizer_maxiter = int(optimizer_maxiter)
    self.numerical_jitter = float(numerical_jitter)
    self.verbose = bool(verbose)
    self._clear_fit_state()

fit

fit(
    X: Union[ndarray, DataFrame],
    y: Union[ndarray, Series],
    coords: Union[ndarray, DataFrame],
) -> "ScalableGWR"

Fit the published ScaGWR estimator.

Source code in src/pygwrx/models/scalable_gwr.py

def fit(
    self,
    X: Union[np.ndarray, pd.DataFrame],
    y: Union[np.ndarray, pd.Series],
    coords: Union[np.ndarray, pd.DataFrame],
) -> "ScalableGWR":
    """Fit the published ScaGWR estimator."""
    self._clear_fit_state()
    try:
        X_array, columns = self._as_2d_numeric(X, name="X")
        y_array = self._as_1d_numeric(y)
        coords_array = validate_coords(coords)
        if (
            X_array.shape[0] != y_array.size
            or X_array.shape[0] != coords_array.shape[0]
        ):
            raise ValueError(
                "X, y, and coords must contain the same number of rows."
            )
        n = X_array.shape[0]
        if self.bandwidth >= n:
            raise ValueError("bandwidth must be smaller than n_samples.")

        design, design_names = self._prepare_design(X_array, columns)
        if self.bandwidth < design.shape[1] + 1:
            raise ValueError(
                "bandwidth must exceed the number of design columns for "
                "stable local fitting."
            )
        self.n_features_in_ = X_array.shape[1]
        self.feature_names_in_ = columns or tuple(
            f"x{i}" for i in range(X_array.shape[1])
        )
        self.design_feature_names_ = design_names
        self.X_train_ = design
        self.y_train_ = y_array
        self.coords_train_ = coords_array
        self.bandwidth_ = self.bandwidth
        self.global_cross_product_ = design.T @ design
        self.global_response_product_ = design.T @ y_array
        self._tree_ = cKDTree(coords_array)

        all_indices = np.arange(n, dtype=int)
        nonself_indices, nonself_distances = self._training_neighbors(
            coords_array, all_indices, exclude_self=True
        )
        rank = min(50, self.bandwidth) - 1
        reference_distance = float(np.median(nonself_distances[:, rank]))
        if reference_distance <= 0 or not np.isfinite(reference_distance):
            positive = nonself_distances[nonself_distances > 0]
            if positive.size == 0:
                raise ValueError(
                    "Coordinates do not provide positive neighbour distances."
                )
            reference_distance = float(np.median(positive))
        self.base_bandwidth_ = (
            reference_distance / np.sqrt(3.0)
            if self.kernel == "gaussian"
            else reference_distance / 3.0
        )

        if self.criterion == "aicc" and self.sample_size is not None:
            warnings.warn(
                "sample_size is ignored for AICc calibration.",
                RuntimeWarning,
                stacklevel=2,
            )
            calibration_indices = all_indices
        elif self.sample_size is not None and self.sample_size < n:
            rng = np.random.default_rng(self.random_state)
            calibration_indices = np.sort(
                rng.choice(n, size=self.sample_size, replace=False)
            )
        else:
            calibration_indices = all_indices

        if self.criterion == "cv":
            cv_neighbors, cv_distances = self._training_neighbors(
                coords_array, calibration_indices, exclude_self=True
            )
            calibration_moments = self._compress(
                design[calibration_indices],
                coords_array[calibration_indices],
                cv_neighbors,
                cv_distances,
                target_y=y_array[calibration_indices],
                need_squared=False,
            )
        else:
            fit_neighbors, fit_distances = self._training_neighbors(
                coords_array, calibration_indices, exclude_self=False
            )
            calibration_moments = self._compress(
                design[calibration_indices],
                coords_array[calibration_indices],
                fit_neighbors,
                fit_distances,
                target_y=y_array[calibration_indices],
                need_squared=True,
            )

        self.scale_, self.penalty_ = self._calibrate(calibration_moments)
        log_scale = float(np.log(self.scale_))
        log_penalty = float(np.log(max(self.penalty_, 1.0e-300)))

        final_neighbors, final_distances = self._training_neighbors(
            coords_array, all_indices, exclude_self=False
        )
        final_moments = self._compress(
            design,
            coords_array,
            final_neighbors,
            final_distances,
            target_y=y_array,
            need_squared=True,
        )
        systems, rhs, second = self._assemble(
            final_moments, log_scale, log_penalty, inference=True
        )
        beta = self._solve_batch(systems, rhs)
        fitted = np.einsum("ni,ni->n", design, beta)
        residuals = y_array - fitted
        rss = float(residuals @ residuals)

        inverse_x = np.linalg.solve(systems, design[..., None])[..., 0]
        trace_s = float(
            (1.0 + self.penalty_) * np.einsum("ni,ni->", design, inverse_x)
        )
        trace_sts = float(
            np.einsum("ni,nij,nj->", inverse_x, second, inverse_x, optimize=True)
        )
        enp = float(2.0 * trace_s - trace_sts)
        edf = float(n - enp)
        if edf <= 0:
            raise ValueError(
                "ScaGWR effective residual degrees of freedom are non-positive."
            )
        sigma = float(np.sqrt(rss / edf))

        inverse_systems = np.linalg.inv(systems)
        covariance = np.einsum(
            "nij,njk,nlk->nil",
            inverse_systems,
            second,
            inverse_systems,
            optimize=True,
        )
        standard_errors = sigma * np.sqrt(
            np.maximum(np.diagonal(covariance, axis1=1, axis2=2), 0.0)
        )
        with np.errstate(divide="ignore", invalid="ignore"):
            t_values = beta / standard_errors
        p_values = 2.0 * student_t.sf(np.abs(t_values), df=edf)

        full_cv_neighbors, full_cv_distances = self._training_neighbors(
            coords_array, all_indices, exclude_self=True
        )
        full_cv_moments = self._compress(
            design,
            coords_array,
            full_cv_neighbors,
            full_cv_distances,
            target_y=y_array,
            need_squared=False,
        )
        cv_rss = self._evaluate_params(
            full_cv_moments,
            np.array([log_scale, log_penalty]),
            criterion="cv",
        )
        self.cv_score_ = float(np.sqrt(cv_rss / n))

        sigma_ml = np.sqrt(max(rss / n, np.finfo(float).tiny))
        log_likelihood = float(
            -n * np.log(sigma_ml) - n / 2.0 * np.log(2.0 * np.pi)
        )
        aic = float(-2.0 * log_likelihood + n + trace_s)
        aicc = (
            float(-2.0 * log_likelihood + n * (n + trace_s) / (n - 2.0 - trace_s))
            if n - 2.0 - trace_s > 0
            else np.inf
        )
        tss = float(np.sum(np.square(y_array - np.mean(y_array))))
        rss_r2 = 1.0 - rss / tss if tss > 0 else np.nan
        corr = np.corrcoef(y_array, fitted)[0, 1]
        r2 = float(corr * corr) if np.isfinite(corr) else float(rss_r2)
        adjusted_r2 = float(1.0 - (1.0 - r2) * (n - 1.0) / (n - enp - 1.0))

        self.coefficients_ = beta
        if self.fit_intercept:
            self.intercept_ = beta[:, 0]
            self.coef_ = beta[:, 1:]
            self.intercept_standard_errors_ = standard_errors[:, 0]
            self.coef_standard_errors_ = standard_errors[:, 1:]
        else:
            self.intercept_ = np.zeros(n)
            self.coef_ = beta
            self.intercept_standard_errors_ = np.zeros(n)
            self.coef_standard_errors_ = standard_errors
        self.standard_errors_ = standard_errors
        self.t_values_ = t_values
        self.p_values_ = p_values
        self.fitted_values_ = fitted
        self.residuals_ = residuals
        self.trace_S_ = trace_s
        self.trace_StS_ = trace_sts
        self.effective_n_params_ = enp
        self.effective_df_ = edf
        self.sigma_ = sigma
        self.aic_ = aic
        self.aicc_ = aicc
        self.r2_ = r2
        self.adjusted_r2_ = adjusted_r2
        self.diagnostics_ = {
            "rss": rss,
            "sigma": sigma,
            "trace_S": trace_s,
            "trace_StS": trace_sts,
            "enp": enp,
            "edf": edf,
            "r2": r2,
            "rss_r2": float(rss_r2),
            "adjusted_r2": adjusted_r2,
            "aic": aic,
            "aicc": aicc,
            "cv_rmse": self.cv_score_,
            "scale": self.scale_,
            "penalty": self.penalty_,
            "base_bandwidth": self.base_bandwidth_,
            "n_neighbors": self.bandwidth_,
        }
        self._is_fitted = True

        if self.verbose:
            print(
                "ScalableGWR fitted: "
                f"n={n}, Q={self.bandwidth_}, scale={self.scale_:.6g}, "
                f"penalty={self.penalty_:.6g}, CV_RMSE={self.cv_score_:.6g}"
            )
        return self
    except Exception:
        self._clear_fit_state()
        raise

predict_result

predict_result(
    X: Optional[Union[ndarray, DataFrame]],
    coords: Union[ndarray, DataFrame],
    *,
    return_standard_errors: bool = False
) -> ScalableGWRPredictionResult

Estimate coefficients and optionally predictions at new locations.

Source code in src/pygwrx/models/scalable_gwr.py

def predict_result(
    self,
    X: Optional[Union[np.ndarray, pd.DataFrame]],
    coords: Union[np.ndarray, pd.DataFrame],
    *,
    return_standard_errors: bool = False,
) -> ScalableGWRPredictionResult:
    """Estimate coefficients and optionally predictions at new locations."""
    if not self._is_fitted:
        raise ValueError("Model not fitted. Call fit() first.")
    coords_array = validate_coords(coords)
    if X is None:
        target_design = np.zeros((coords_array.shape[0], self.X_train_.shape[1]))
        predictions = None
    else:
        target_design = self._check_prediction_features(X)
        if target_design.shape[0] != coords_array.shape[0]:
            raise ValueError("X and coords must contain the same number of rows.")
        predictions = np.empty(coords_array.shape[0], dtype=float)

    distances, indices = self._tree_.query(coords_array, k=self.bandwidth_)
    if self.bandwidth_ == 1:
        distances = distances[:, None]
        indices = indices[:, None]
    moments = self._compress(
        target_design,
        coords_array,
        np.asarray(indices, dtype=int),
        np.asarray(distances, dtype=float),
        target_y=None,
        need_squared=return_standard_errors,
    )
    systems, rhs, second = self._assemble(
        moments,
        float(np.log(self.scale_)),
        float(np.log(max(self.penalty_, 1.0e-300))),
        inference=return_standard_errors,
    )
    beta = self._solve_batch(systems, rhs)
    if predictions is not None:
        predictions[:] = np.einsum("ni,ni->n", target_design, beta)

    standard_errors = None
    if return_standard_errors:
        inverse_systems = np.linalg.inv(systems)
        covariance = np.einsum(
            "nij,njk,nlk->nil",
            inverse_systems,
            second,
            inverse_systems,
            optimize=True,
        )
        standard_errors = self.sigma_ * np.sqrt(
            np.maximum(np.diagonal(covariance, axis1=1, axis2=2), 0.0)
        )
    return ScalableGWRPredictionResult(
        predictions=predictions,
        coefficients=beta,
        standard_errors=standard_errors,
        coords=coords_array,
        feature_names=self.design_feature_names_,
    )

predict

predict(
    X: Union[ndarray, DataFrame],
    coords: Union[ndarray, DataFrame],
) -> np.ndarray

Predict responses by estimating ScaGWR coefficients at new locations.

Source code in src/pygwrx/models/scalable_gwr.py

def predict(
    self,
    X: Union[np.ndarray, pd.DataFrame],
    coords: Union[np.ndarray, pd.DataFrame],
) -> np.ndarray:
    """Predict responses by estimating ScaGWR coefficients at new locations."""
    result = self.predict_result(X, coords)
    return np.asarray(result.predictions)

to_frame

to_frame() -> pd.DataFrame

Return training-location coefficients, inference, and fit diagnostics.

Source code in src/pygwrx/models/scalable_gwr.py

def to_frame(self) -> pd.DataFrame:
    """Return training-location coefficients, inference, and fit diagnostics."""
    if not self._is_fitted:
        raise ValueError("Model not fitted. Call fit() first.")
    data: Dict[str, np.ndarray] = {
        "coord_0": self.coords_train_[:, 0],
        "coord_1": self.coords_train_[:, 1],
        "observed": self.y_train_,
        "fitted": self.fitted_values_,
        "residual": self.residuals_,
    }
    for index, name in enumerate(self.design_feature_names_):
        data[f"coef_{name}"] = self.coefficients_[:, index]
        data[f"se_{name}"] = self.standard_errors_[:, index]
        data[f"t_{name}"] = self.t_values_[:, index]
        data[f"p_{name}"] = self.p_values_[:, index]
    return pd.DataFrame(data)

summary

summary() -> str

Return fitted diagnostics as a plain-text table.

Source code in src/pygwrx/models/scalable_gwr.py

def summary(self) -> str:
    """Return fitted diagnostics as a plain-text table."""
    if not self._is_fitted:
        raise ValueError("Model not fitted. Call fit() first.")
    return format_summary("Scalable GWR Summary", self.diagnostics_)

Runnable examples used on this page

examples/models/12_scalable_gwr.py

# SPDX-FileCopyrightText: 2026 Jinghao Hu
# SPDX-License-Identifier: MIT

"""Fit scalable GWR with a fixed multiscale-kernel approximation."""

# Allow this script to run directly from any working directory.
import sys
from pathlib import Path

_PROJECT_ROOT = Path(__file__).resolve().parents[2]
_EXAMPLES_ROOT = _PROJECT_ROOT / "examples"
_SRC_ROOT = _PROJECT_ROOT / "src"
for _path in (_SRC_ROOT, _EXAMPLES_ROOT):
    if str(_path) not in sys.path:
        sys.path.insert(0, str(_path))

from _common import print_model_result, spatial_regression

from pygwrx import ScalableGWR

X, y, coords = spatial_regression(n=54, p=2)
model = ScalableGWR(
    bandwidth=24, optimize_bandwidth=False, polynomial=4, random_state=0
).fit(X, y, coords)
print_model_result(model)
print("predictions=", model.predict(X.iloc[:3], coords.iloc[:3]))