Local solvers

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

Conceptual guide

weighted_least_squares

Solve a weighted least-squares problem.

Property	Value
Type	function
Import	from pygwrx.core import weighted_least_squares
Signature	weighted_least_squares(X: 'np.ndarray', y: 'np.ndarray', weights: 'np.ndarray', *, ridge: 'float' = 1e-08) -> 'Tuple[np.ndarray, np.ndarray]'
Maintained example	examples/core/04_solver.py

weighted_least_squares

weighted_least_squares(
    X: ndarray,
    y: ndarray,
    weights: ndarray,
    *,
    ridge: float = _DEFAULT_RIDGE
) -> Tuple[np.ndarray, np.ndarray]

Solve a weighted least-squares problem.

Parameters:

Name	Type	Description	Default
X	ndarray	Regression design matrix.	required
y	ndarray	Single-response target vector.	required
weights	ndarray	Non-negative observation weights. Exact zeros are preserved, so zero-weight observations are genuinely excluded from the local objective.	required
ridge	float	Non-negative diagonal regularization used to stabilize the local normal matrix.	_DEFAULT_RIDGE

Returns:

Name	Type	Description
beta	ndarray	Estimated coefficients.
inverse_normal_matrix	ndarray	Inverse (or pseudo-inverse) of X.T @ W @ X + ridge * I. This is the unscaled coefficient covariance factor; a statistical covariance matrix also requires multiplication by an appropriate residual-variance estimate.

Notes

The function solves

min_beta  sum_i weights[i] * (y[i] - X[i] @ beta)^2
            + ridge * ||beta||_2^2.

Unlike the original implementation, zero weights are not replaced by 1e-10. This preserves compact-support kernels and strict leave-one-out calculations.

Source code in src/pygwrx/core/solver.py

def weighted_least_squares(
    X: np.ndarray,
    y: np.ndarray,
    weights: np.ndarray,
    *,
    ridge: float = _DEFAULT_RIDGE,
) -> Tuple[np.ndarray, np.ndarray]:
    """Solve a weighted least-squares problem.

    Args:
        X: Regression design matrix.
        y: Single-response target vector.
        weights: Non-negative observation weights. Exact zeros are preserved, so zero-weight
            observations are genuinely excluded from the local objective.
        ridge: Non-negative diagonal regularization used to stabilize the local normal matrix.

    Returns:
        beta: Estimated coefficients.
        inverse_normal_matrix: Inverse (or pseudo-inverse) of ``X.T @ W @ X + ridge * I``. This is the
            unscaled coefficient covariance factor; a statistical covariance matrix also
            requires multiplication by an appropriate residual-variance estimate.

    Notes:
        The function solves

            min_beta  sum_i weights[i] * (y[i] - X[i] @ beta)^2
                        + ridge * ||beta||_2^2.

        Unlike the original implementation, zero weights are not replaced by ``1e-10``.
        This preserves compact-support kernels and strict leave-one-out calculations.
    """
    X_arr = _validate_design_matrix(X)
    y_arr = _validate_response(y, X_arr.shape[0])
    weights_arr = _validate_weights(weights, X_arr.shape[0])
    ridge_value = _validate_nonnegative_scalar(ridge, "ridge")

    system, XtWy, _ = _normal_equations(
        X_arr,
        y_arr,
        weights_arr,
        ridge=ridge_value,
    )

    beta = _solve_linear_system(system, XtWy)
    inverse_normal_matrix = _solve_linear_system(
        system,
        np.eye(system.shape[0], dtype=float),
    )

    # Remove tiny asymmetric round-off introduced by the numerical solve.
    inverse_normal_matrix = 0.5 * (inverse_normal_matrix + inverse_normal_matrix.T)

    return beta, inverse_normal_matrix

local_regression

Perform local weighted regression at target locations.

Property	Value
Type	function
Import	from pygwrx.core import local_regression
Signature	local_regression(X: 'np.ndarray', y: 'np.ndarray', coords: 'np.ndarray', target_coords: 'np.ndarray', kernel_func: 'KernelFunction', bandwidth: 'float', distance_metric: 'str' = 'euclidean', adaptive: 'bool' = False, *, ridge: 'float' = 1e-08) -> 'np.ndarray'
Maintained example	examples/core/04_solver.py

local_regression

local_regression(
    X: ndarray,
    y: ndarray,
    coords: ndarray,
    target_coords: ndarray,
    kernel_func: KernelFunction,
    bandwidth: float,
    distance_metric: str = "euclidean",
    adaptive: bool = False,
    *,
    ridge: float = _DEFAULT_RIDGE
) -> np.ndarray

Perform local weighted regression at target locations.

Parameters:

Name	Type	Description	Default
X	ndarray	Regression design matrix.	required
y	ndarray	Target vector.	required
coords	ndarray	Training coordinates.	required
target_coords	ndarray	Locations where local coefficients are required.	required
kernel_func	KernelFunction	Function with signature kernel_func(distances, bandwidth) -> weights.	required
bandwidth	float	Fixed distance when adaptive=False; neighbour-order bandwidth k when adaptive=True.	required
distance_metric	str	Distance metric forwarded to compute_distance_matrix.	'euclidean'
adaptive	bool	Whether bandwidth is interpreted as a neighbour-order bandwidth.	False
ridge	float	Non-negative regularization shared with compute_hat_matrix.	_DEFAULT_RIDGE

Returns:

Name	Type	Description
local_coefs	ndarray	Local coefficient vectors.

Notes

If a location has fewer positive-weight observations than design-matrix columns, the function emits a warning and returns the deterministic ridge-regularized solution. It never copies coefficients from a preceding location and never falls back silently to global OLS, so results do not depend on target ordering.

Source code in src/pygwrx/core/solver.py

def local_regression(
    X: np.ndarray,
    y: np.ndarray,
    coords: np.ndarray,
    target_coords: np.ndarray,
    kernel_func: KernelFunction,
    bandwidth: float,
    distance_metric: str = "euclidean",
    adaptive: bool = False,
    *,
    ridge: float = _DEFAULT_RIDGE,
) -> np.ndarray:
    """Perform local weighted regression at target locations.

    Args:
        X: Regression design matrix.
        y: Target vector.
        coords: Training coordinates.
        target_coords: Locations where local coefficients are required.
        kernel_func: Function with signature ``kernel_func(distances, bandwidth) -> weights``.
        bandwidth: Fixed distance when ``adaptive=False``; neighbour-order bandwidth ``k`` when
            ``adaptive=True``.
        distance_metric: Distance metric forwarded to ``compute_distance_matrix``.
        adaptive: Whether ``bandwidth`` is interpreted as a neighbour-order bandwidth.
        ridge: Non-negative regularization shared with ``compute_hat_matrix``.

    Returns:
        local_coefs: Local coefficient vectors.

    Notes:
        If a location has fewer positive-weight observations than design-matrix columns,
        the function emits a warning and returns the deterministic ridge-regularized
        solution. It never copies coefficients from a preceding location and never falls
        back silently to global OLS, so results do not depend on target ordering.
    """
    from pygwrx.core.utils import compute_distance_matrix

    X_arr = _validate_design_matrix(X)
    y_arr = _validate_response(y, X_arr.shape[0])
    coords_arr = _validate_coordinates(
        coords,
        name="coords",
        expected_rows=X_arr.shape[0],
    )
    target_arr = _validate_coordinates(
        target_coords,
        name="target_coords",
        expected_dimension=coords_arr.shape[1],
    )
    kernel = _validate_kernel(kernel_func)
    ridge_value = _validate_nonnegative_scalar(ridge, "ridge")

    if not isinstance(adaptive, (bool, np.bool_)):
        raise TypeError("adaptive must be a boolean.")
    adaptive_value = bool(adaptive)

    if adaptive_value:
        bandwidth_value: Union[float, int] = _validate_adaptive_k(
            bandwidth,
            X_arr.shape[0],
        )
    else:
        bandwidth_value = _validate_fixed_bandwidth(bandwidth)

    distances = np.asarray(
        compute_distance_matrix(
            target_arr,
            coords_arr,
            metric=distance_metric,
        ),
        dtype=float,
    )

    expected_shape = (target_arr.shape[0], coords_arr.shape[0])
    if distances.shape != expected_shape:
        raise ValueError(
            "compute_distance_matrix returned an unexpected shape; "
            f"expected {expected_shape}, got {distances.shape}."
        )
    if not np.all(np.isfinite(distances)) or np.any(distances < 0):
        raise ValueError("Distance calculation produced invalid values.")

    local_coefs = np.empty(
        (target_arr.shape[0], X_arr.shape[1]),
        dtype=float,
    )

    for i, dists in enumerate(distances):
        try:
            weights = _compute_kernel_weights(
                dists,
                kernel,
                bandwidth_value,
                adaptive=adaptive_value,
            )

            n_positive = int(np.count_nonzero(weights > 0))
            if n_positive < X_arr.shape[1]:
                warnings.warn(
                    f"Location {i}: only {n_positive} positive-weight observations "
                    f"are available for {X_arr.shape[1]} design-matrix columns. "
                    "Returning a ridge-regularized local solution; consider increasing "
                    "the bandwidth.",
                    RuntimeWarning,
                    stacklevel=2,
                )

            beta, _ = weighted_least_squares(
                X_arr,
                y_arr,
                weights,
                ridge=ridge_value,
            )
        except (ValueError, TypeError, np.linalg.LinAlgError) as exc:
            raise RuntimeError(
                f"Local regression failed at target location {i}: {exc}"
            ) from exc

        local_coefs[i] = beta

    return local_coefs

compute_hat_matrix

Compute the GWR hat matrix S such that y_hat = S @ y.

Property	Value
Type	function
Import	from pygwrx.core import compute_hat_matrix
Signature	compute_hat_matrix(X: 'np.ndarray', coords: 'np.ndarray', kernel_func: 'KernelFunction', bandwidth: 'float', distance_metric: 'str' = 'euclidean', adaptive: 'bool' = False, *, ridge: 'float' = 1e-08) -> 'np.ndarray'
Maintained example	examples/core/04_solver.py

compute_hat_matrix

compute_hat_matrix(
    X: ndarray,
    coords: ndarray,
    kernel_func: KernelFunction,
    bandwidth: float,
    distance_metric: str = "euclidean",
    adaptive: bool = False,
    *,
    ridge: float = _DEFAULT_RIDGE
) -> np.ndarray

Compute the GWR hat matrix S such that y_hat = S @ y.

Parameters:

Name	Type	Description	Default
X	ndarray	Regression design matrix.	required
coords	ndarray	Training coordinates.	required
kernel_func	KernelFunction	Spatial kernel function.	required
bandwidth	float	Fixed distance or adaptive neighbour-order bandwidth.	required
distance_metric	str	Distance metric forwarded to compute_distance_matrix.	'euclidean'
adaptive	bool	Whether bandwidth is an adaptive neighbour-order bandwidth.	False
ridge	float	Non-negative regularization. The same value and normal-system construction are used by weighted_least_squares and local_regression.	_DEFAULT_RIDGE

Returns:

Name	Type	Description
hat_matrix	ndarray	Full smoother matrix.

Notes

The full matrix requires 8 * n_samples**2 bytes for float64 storage, excluding the distance matrix and temporary arrays. Large-data models should eventually use trace-only or chunked diagnostics instead.

Source code in src/pygwrx/core/solver.py

def compute_hat_matrix(
    X: np.ndarray,
    coords: np.ndarray,
    kernel_func: KernelFunction,
    bandwidth: float,
    distance_metric: str = "euclidean",
    adaptive: bool = False,
    *,
    ridge: float = _DEFAULT_RIDGE,
) -> np.ndarray:
    """Compute the GWR hat matrix ``S`` such that ``y_hat = S @ y``.

    Args:
        X: Regression design matrix.
        coords: Training coordinates.
        kernel_func: Spatial kernel function.
        bandwidth: Fixed distance or adaptive neighbour-order bandwidth.
        distance_metric: Distance metric forwarded to ``compute_distance_matrix``.
        adaptive: Whether ``bandwidth`` is an adaptive neighbour-order bandwidth.
        ridge: Non-negative regularization. The same value and normal-system construction are
            used by ``weighted_least_squares`` and ``local_regression``.

    Returns:
        hat_matrix: Full smoother matrix.

    Notes:
        The full matrix requires ``8 * n_samples**2`` bytes for float64 storage, excluding
        the distance matrix and temporary arrays. Large-data models should eventually use
        trace-only or chunked diagnostics instead.
    """
    from pygwrx.core.utils import compute_distance_matrix

    X_arr = _validate_design_matrix(X)
    coords_arr = _validate_coordinates(
        coords,
        name="coords",
        expected_rows=X_arr.shape[0],
    )
    kernel = _validate_kernel(kernel_func)
    ridge_value = _validate_nonnegative_scalar(ridge, "ridge")

    if not isinstance(adaptive, (bool, np.bool_)):
        raise TypeError("adaptive must be a boolean.")
    adaptive_value = bool(adaptive)

    if adaptive_value:
        bandwidth_value: Union[float, int] = _validate_adaptive_k(
            bandwidth,
            X_arr.shape[0],
        )
    else:
        bandwidth_value = _validate_fixed_bandwidth(bandwidth)

    distances = np.asarray(
        compute_distance_matrix(
            coords_arr,
            coords_arr,
            metric=distance_metric,
        ),
        dtype=float,
    )

    expected_shape = (X_arr.shape[0], X_arr.shape[0])
    if distances.shape != expected_shape:
        raise ValueError(
            "compute_distance_matrix returned an unexpected shape; "
            f"expected {expected_shape}, got {distances.shape}."
        )
    if not np.all(np.isfinite(distances)) or np.any(distances < 0):
        raise ValueError("Distance calculation produced invalid values.")

    hat_matrix = np.empty(expected_shape, dtype=float)

    # A dummy y is used only to reuse the identical weighted normal-system builder.
    dummy_y = np.zeros(X_arr.shape[0], dtype=float)

    for i, dists in enumerate(distances):
        try:
            weights = _compute_kernel_weights(
                dists,
                kernel,
                bandwidth_value,
                adaptive=adaptive_value,
            )

            n_positive = int(np.count_nonzero(weights > 0))
            if n_positive < X_arr.shape[1]:
                warnings.warn(
                    f"Location {i}: only {n_positive} positive-weight observations "
                    f"are available for {X_arr.shape[1]} design-matrix columns. "
                    "The corresponding hat-matrix row is ridge regularized.",
                    RuntimeWarning,
                    stacklevel=2,
                )

            system, _, XtW = _normal_equations(
                X_arr,
                dummy_y,
                weights,
                ridge=ridge_value,
            )

            # X_i @ inv(system) @ XtW, computed without explicitly inverting system.
            left = _solve_linear_system(system.T, X_arr[i])
            hat_row = left @ XtW

            if not np.all(np.isfinite(hat_row)):
                raise np.linalg.LinAlgError(
                    "The hat-matrix row contains non-finite values."
                )
        except (ValueError, TypeError, np.linalg.LinAlgError) as exc:
            raise RuntimeError(
                f"Hat-matrix computation failed at location {i}: {exc}"
            ) from exc

        hat_matrix[i] = hat_row

    return hat_matrix

adaptive_bandwidth_weights

Convert an adaptive neighbour-order bandwidth into a distance scale.

Property	Value
Type	function
Import	from pygwrx.core import adaptive_bandwidth_weights
Signature	adaptive_bandwidth_weights(distances: 'np.ndarray', k_nearest: 'int') -> 'float'
Maintained example	examples/core/04_solver.py

adaptive_bandwidth_weights

adaptive_bandwidth_weights(
    distances: ndarray, k_nearest: int
) -> float

Convert an adaptive neighbour-order bandwidth into a distance scale.

Parameters:

Name	Type	Description	Default
distances	ndarray	Non-negative distances from one regression location to all observations.	required
k_nearest	int	One-based neighbour order used to determine the local distance scale. The current PyGWRx convention includes a zero-distance self observation when the regression location is one of the training locations.	required

Returns:

Name	Type	Description
bandwidth	float	Strictly positive distance scale corresponding to the k-th ordered distance.

Notes

np.partition is used for expected O(n) selection. If duplicate coordinates put the requested ordered distance at zero, the smallest positive distance is used. The result is advanced by one representable float with np.nextafter so compact kernels assign a positive (possibly tiny) weight to the boundary neighbour instead of excluding it exactly at d == bandwidth.

Source code in src/pygwrx/core/solver.py

def adaptive_bandwidth_weights(
    distances: np.ndarray,
    k_nearest: int,
) -> float:
    """Convert an adaptive neighbour-order bandwidth into a distance scale.

    Args:
        distances: Non-negative distances from one regression location to all observations.
        k_nearest: One-based neighbour order used to determine the local distance scale. The
            current PyGWRx convention includes a zero-distance self observation when the
            regression location is one of the training locations.

    Returns:
        bandwidth: Strictly positive distance scale corresponding to the k-th ordered distance.

    Notes:
        ``np.partition`` is used for expected O(n) selection. If duplicate coordinates put
        the requested ordered distance at zero, the smallest positive distance is used. The
        result is advanced by one representable float with ``np.nextafter`` so compact
        kernels assign a positive (possibly tiny) weight to the boundary neighbour instead
        of excluding it exactly at ``d == bandwidth``.
    """
    try:
        distances_arr = np.asarray(distances, dtype=float)
    except (TypeError, ValueError) as exc:
        raise TypeError("distances must contain numeric values.") from exc

    if distances_arr.ndim != 1:
        raise ValueError("distances must be a one-dimensional array.")
    if distances_arr.size == 0:
        raise ValueError("distances must contain at least one value.")
    if not np.all(np.isfinite(distances_arr)):
        raise ValueError("distances contain NaN or infinite values.")
    if np.any(distances_arr < 0):
        raise ValueError("distances must be non-negative.")

    k_value = _validate_adaptive_k(k_nearest, distances_arr.size)
    distance_bandwidth = float(np.partition(distances_arr, k_value - 1)[k_value - 1])

    if distance_bandwidth <= 0:
        positive_distances = distances_arr[distances_arr > 0]
        if positive_distances.size == 0:
            raise ValueError(
                "Adaptive bandwidth is undefined because all distances are zero."
            )
        distance_bandwidth = float(np.min(positive_distances))

    return float(np.nextafter(distance_bandwidth, np.inf))

Runnable examples used on this page

examples/core/04_solver.py

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

"""Run all public local-regression solver utilities."""

# 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))

import numpy as np

from pygwrx.core import (
    adaptive_bandwidth_weights,
    compute_hat_matrix,
    gaussian_kernel,
    local_regression,
    weighted_least_squares,
)

rng = np.random.default_rng(0)
coords = rng.uniform(0.0, 5.0, size=(20, 2))
x = rng.normal(size=20)
X = np.column_stack((np.ones(20), x))
y = 1.0 + 2.0 * x + rng.normal(0.0, 0.05, 20)
distances = np.linalg.norm(coords - coords[0], axis=1)
weights = gaussian_kernel(distances, bandwidth=2.0)
beta, covariance = weighted_least_squares(X, y, weights)
print("beta=", beta)
print("covariance_shape=", covariance.shape)
print("adaptive_scale=", adaptive_bandwidth_weights(distances, 8))
print(
    "local_parameters=",
    local_regression(X, y, coords, coords[:3], gaussian_kernel, 2.0),
)
hat = compute_hat_matrix(X, coords, gaussian_kernel, 2.0)
print("hat_shape_trace=", hat.shape, np.trace(hat))