Kernel Bw

bw

FUNCTION	DESCRIPTION
bw_compute	Compute bandwidth (squared) for Gaussian kernel from pairwise distance

Functions

bw_compute(D: torch.Tensor, bw_method: str = 'median', bw_kwargs: dict | None = None, weights: torch.Tensor | None = None)

Compute bandwidth (squared) for Gaussian kernel from pairwise distance matrix.

When weights are provided the bandwidth is estimated under the reweighted distribution. Each pair \((i, j)\) is weighted by \(w_i w_j\). Implemented only for bandwidths estimated using quantiles.

PARAMETER	DESCRIPTION
D	Pairwise distance matrix of shape (n, m). TYPE: Tensor
bw_method	Method to compute bandwidth. "median": Median of all distances. "median_1d": Median per row and column, returned as outer product. "quantile": Quantile of distances specified by "q" in bw_kwargs. "quantile_1d": Quantile per row and col specified by "q" in bw_kwargs. "quantile_search": Bandwidth from quantile search method. "knn": Bandwidth from k-nearest neighbors specified by "k_neighbor". "bgh": Bandwidth using BGH method [1]. "median_scaled": Median with dimension scaling d^alpha. "median_1d_scaled": Per-dimension median with dimension scaling. "silverman": Silverman's rule. "scott": Scott's rule. TYPE: str, optional, by default "median" DEFAULT: 'median'
bw_kwargs	Additional parameters for specific bandwidth methods. Options are: For "quantile" and "quantile_1d": "q" (float): Quantile value in [0.0, 1.0], default 0.5. For "quantile_search" "n_quantiles" (int): Number of quantiles to evaluate, default 100. "eps" (float): Epsilon to exclude from quantiles, default 1e-3. For "knn": "k_neighbor" (int): Number of nearest neighbors, default 5. For "bgh": "bw_values" (torch.Tensor | None): Bandwidth values to test, default None. For "median_scaled" and "median_1d_scaled": "d" (int): Feature dimensionality (required). "alpha" (float): Scaling exponent, default 0.5. Typical range: 0.2-1.0. For "silverman" and "scott": "d" (int): Feature dimensionality (required). "n" (int | None): Sample size, default None (inferred from D). TYPE: dict DEFAULT: None
weights	Sample weights, shape (n,). TYPE: torch.Tensor | None, optional, by default None DEFAULT: None

RETURNS	DESCRIPTION
Tensor	Squared bandwidth value or matrix depending on the selected method.

Examples:

Standard bandwidth methods:

>>> import torch
>>> from spectre.kernel import bw_compute
>>> D = torch.tensor([[0.0, 1.0, 2.0], [1.0, 0.0, 1.5], [2.0, 1.5, 0.0]])
>>> bw = bw_compute(D, bw_method="median")
>>> print(bw)
tensor(1.0000)

>>> bw_matrix = bw_compute(D, bw_method="median_1d")
>>> print(bw_matrix)
tensor([[1.0000, 1.0000, 1.5000],
        [1.0000, 1.0000, 1.5000],
        [1.5000, 1.5000, 2.2500]])

>>> bw = bw_compute(D, bw_method="quantile", bw_kwargs={"q": 0.9})
>>> bw = bw_compute(D, bw_method="knn", bw_kwargs={"k_neighbor": 2})

Dimension-aware methods for feature selection:

>>> # Data has 50 features
>>> bw = bw_compute(
...     D,
...     bw_method="median_scaled",
...     bw_kwargs={"d": 50, "alpha": 0.5},
... )
>>> bw = bw_compute(D, bw_method="silverman", bw_kwargs={"d": 50})
>>> bw = bw_compute(D, bw_method="scott", bw_kwargs={"d": 50, "n": 100})

Source code in spectre/kernel/bw.py

def bw_compute(
    D: torch.Tensor,
    bw_method: str = "median",
    bw_kwargs: dict | None = None,
    weights: torch.Tensor | None = None,
):
    """
    Compute bandwidth (squared) for Gaussian kernel from pairwise distance
    matrix.

    When `weights` are provided the bandwidth is estimated under the reweighted
    distribution. Each pair $(i, j)$ is weighted by $w_i w_j$. Implemented only
    for bandwidths estimated using quantiles.

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix of shape (n, m).

    bw_method : str, optional, by default "median"
        Method to compute bandwidth.

        - "median": Median of all distances.
        - "median_1d": Median per row and column, returned as outer product.
        - "quantile": Quantile of distances specified by "q" in `bw_kwargs`.
        - "quantile_1d": Quantile per row and col specified by "q" in `bw_kwargs`.
        - "quantile_search": Bandwidth from quantile search method.
        - "knn": Bandwidth from k-nearest neighbors specified by "k_neighbor".
        - "bgh": Bandwidth using BGH method [@berry2015nonparametric].
        - "median_scaled": Median with dimension scaling d^alpha.
        - "median_1d_scaled": Per-dimension median with dimension scaling.
        - "silverman": Silverman's rule.
        - "scott": Scott's rule.

    bw_kwargs : dict, optional
        Additional parameters for specific bandwidth methods. Options are:

        - For "quantile" and "quantile_1d":
            - "q" (float): Quantile value in [0.0, 1.0], default 0.5.

        - For "quantile_search"
            - "n_quantiles" (int): Number of quantiles to evaluate, default
              100.
            - "eps" (float): Epsilon to exclude from quantiles, default 1e-3.

        - For "knn":
            - "k_neighbor" (int): Number of nearest neighbors, default 5.

        - For "bgh":
            - "bw_values" (torch.Tensor | None): Bandwidth values to test,
              default None.

        - For "median_scaled" and "median_1d_scaled":
            - "d" (int): Feature dimensionality (required).
            - "alpha" (float): Scaling exponent, default 0.5. Typical range:
              0.2-1.0.

        - For "silverman" and "scott":
            - "d" (int): Feature dimensionality (required).
            - "n" (int | None): Sample size, default None (inferred from D).

    weights : torch.Tensor | None, optional, by default None
        Sample weights, shape (n,).

    Returns
    -------
    torch.Tensor
        Squared bandwidth value or matrix depending on the selected method.

    Examples
    --------
    Standard bandwidth methods:

    >>> import torch
    >>> from spectre.kernel import bw_compute
    >>> D = torch.tensor([[0.0, 1.0, 2.0], [1.0, 0.0, 1.5], [2.0, 1.5, 0.0]])
    >>> bw = bw_compute(D, bw_method="median")
    >>> print(bw)
    tensor(1.0000)

    >>> bw_matrix = bw_compute(D, bw_method="median_1d")
    >>> print(bw_matrix)
    tensor([[1.0000, 1.0000, 1.5000],
            [1.0000, 1.0000, 1.5000],
            [1.5000, 1.5000, 2.2500]])

    >>> bw = bw_compute(D, bw_method="quantile", bw_kwargs={"q": 0.9})
    >>> bw = bw_compute(D, bw_method="knn", bw_kwargs={"k_neighbor": 2})

    Dimension-aware methods for feature selection:

    >>> # Data has 50 features
    >>> bw = bw_compute(
    ...     D,
    ...     bw_method="median_scaled",
    ...     bw_kwargs={"d": 50, "alpha": 0.5},
    ... )
    >>> bw = bw_compute(D, bw_method="silverman", bw_kwargs={"d": 50})
    >>> bw = bw_compute(D, bw_method="scott", bw_kwargs={"d": 50, "n": 100})
    """
    check_2d(D)

    if bw_kwargs is None:
        bw_kwargs = {}

    if bw_method == "median":
        return bw_median(D, weights=weights)

    elif bw_method == "median_1d":
        return bw_median_1d(D, weights=weights)

    elif bw_method in ["quantile", "quantile_1d"]:
        q = bw_kwargs.get("q", 0.5)
        if not isinstance(q, float):
            raise TypeError(f"q must be a float, got {type(q)}.")
        check_in_interval(q, "[0.0, 1.0]")

        return (
            bw_quantile(D, q=q, weights=weights)
            if bw_method == "quantile"
            else bw_quantile_1d(D, q=q, weights=weights)
        )

    elif bw_method == "quantile_search":
        n_quantiles = bw_kwargs.get("n_quantiles", 100)
        if not isinstance(n_quantiles, int):
            raise TypeError(f"n_quantiles must be an integer, got {type(n_quantiles)}.")
        check_in_interval(n_quantiles, "[10, inf)")

        eps = bw_kwargs.get("eps", 1e-3)
        if not isinstance(eps, float):
            raise TypeError(f"eps must be a float, got {type(eps)}.")
        check_in_interval(eps, "(0.0, 0.5)")

        return bw_quantile_search(D, n_quantiles=n_quantiles, eps=eps, weights=weights)

    elif bw_method == "knn":
        k_neighbor = bw_kwargs.get("k_neighbor", 5)
        if not isinstance(k_neighbor, int):
            raise TypeError(f"k must be an integer, got {type(k_neighbor)}.")
        n, m = D.shape
        check_in_interval(k_neighbor, f"[1, {min(m, n)})")

        return bw_k_nearest_neighbors(D, k_neighbor=k_neighbor)

    elif bw_method == "bgh":
        bw_values = bw_kwargs.get("bw_values", None)
        if bw_values is not None and not isinstance(bw_values, torch.Tensor):
            raise TypeError(f"bw_values must be a torch.Tensor, got {type(bw_values)}.")

        return bw_bgh(D, bw_values=bw_values)

    elif bw_method in ["median_scaled", "median_1d_scaled"]:
        d = bw_kwargs.get("d", None)
        if d is None:
            raise ValueError(
                f"'{bw_method}' requires 'd' (feature dim.) in `bw_kwargs`."
            )
        if not isinstance(d, int):
            raise TypeError(f"d must be an integer, got {type(d)}.")
        check_in_interval(d, "[1, inf)")

        alpha = bw_kwargs.get("alpha", 0.5)
        if not isinstance(alpha, (int, float)):
            raise TypeError(f"alpha must be a float, got {type(alpha)}.")
        check_in_interval(float(alpha), "(0.0, inf)")

        return (
            bw_median_scaled(D, d=d, alpha=alpha, weights=weights)
            if bw_method == "median_scaled"
            else bw_median_1d_scaled(D, d=d, alpha=alpha, weights=weights)
        )

    elif bw_method in ["silverman", "scott"]:
        d = bw_kwargs.get("d", None)
        if d is None:
            raise ValueError(
                f"Statistical rule '{bw_method}' requires 'd' (feature dim.) "
                f"in `bw_kwargs`."
            )
        if not isinstance(d, int):
            raise TypeError(f"d must be an integer, got {type(d)}.")
        check_in_interval(d, "[1, inf)")

        n = bw_kwargs.get("n", None)
        if n is not None:
            if not isinstance(n, int):
                raise TypeError(f"n must be an integer, got {type(n)}.")
            check_in_interval(n, "[1, inf)")
        return (
            bw_silverman(D, d=d, n=n, weights=weights)
            if bw_method == "silverman"
            else bw_scott(D, d=d, n=n, weights=weights)
        )

    else:
        raise ValueError(f"Invalid `bw_method` '{bw_method}'.")

bw_median(D: torch.Tensor, weights: torch.Tensor | None = None) -> torch.Tensor

Calculate bandwidth using the median of all pairwise distances.

PARAMETER	DESCRIPTION
D	Pairwise distance matrix. TYPE: Tensor
weights	Sample weights, shape (n,). TYPE: torch.Tensor | None, optional, by default None DEFAULT: None

RETURNS	DESCRIPTION
Tensor	Squared bandwidth value.

Source code in spectre/kernel/bw.py

def bw_median(D: torch.Tensor, weights: torch.Tensor | None = None) -> torch.Tensor:
    """
    Calculate bandwidth using the median of all pairwise distances.

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix.

    weights : torch.Tensor | None, optional, by default None
        Sample weights, shape (n,).

    Returns
    -------
    torch.Tensor
        Squared bandwidth value.
    """
    if weights is not None:
        return quantile(D.flatten(), q=0.5, weights=outer(weights, flatten=True)) ** 2
    else:
        return torch.median(D) ** 2

bw_median_1d(D: torch.Tensor, weights: torch.Tensor | None = None) -> torch.Tensor

Calculate bandwidth using the median per row and column.

PARAMETER	DESCRIPTION
D	Pairwise distance matrix. TYPE: Tensor
weights	Sample weights, shape (n,). TYPE: torch.Tensor | None, optional, by default None DEFAULT: None

RETURNS	DESCRIPTION
Tensor	Bandwidth matrix as outer product of per-row and per-column medians.

Source code in spectre/kernel/bw.py

def bw_median_1d(D: torch.Tensor, weights: torch.Tensor | None = None) -> torch.Tensor:
    """
    Calculate bandwidth using the median per row and column.

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix.

    weights : torch.Tensor | None, optional, by default None
        Sample weights, shape (n,).

    Returns
    -------
    torch.Tensor
        Bandwidth matrix as outer product of per-row and per-column medians.
    """
    if weights is not None:
        bw_rows = quantile(D, q=0.5, weights=weights.unsqueeze(0).expand_as(D), dim=1)
        bw_cols = quantile(D, q=0.5, weights=weights.unsqueeze(1).expand_as(D), dim=0)

        return outer(bw_rows, bw_cols)
    else:
        bw_rows = D.median(dim=1).values
        bw_cols = D.median(dim=0).values

        return outer(bw_rows, bw_cols)

bw_quantile(D: torch.Tensor, q: float, weights: torch.Tensor | None = None) -> torch.Tensor

Calculate bandwidth using the specified quantile of all pairwise distances.

PARAMETER	DESCRIPTION
D	Pairwise distance matrix. TYPE: Tensor
q	Quantile value in [0.0, 1.0]. TYPE: float
weights	Sample weights, shape (n,). TYPE: torch.Tensor | None, optional, by default None DEFAULT: None

RETURNS	DESCRIPTION
Tensor	Squared bandwidth value.

Source code in spectre/kernel/bw.py

def bw_quantile(
    D: torch.Tensor, q: float, weights: torch.Tensor | None = None
) -> torch.Tensor:
    """
    Calculate bandwidth using the specified quantile of all pairwise distances.

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix.

    q : float
        Quantile value in [0.0, 1.0].

    weights : torch.Tensor | None, optional, by default None
        Sample weights, shape (n,).

    Returns
    -------
    torch.Tensor
        Squared bandwidth value.
    """
    if weights is not None:
        w = outer(weights, flatten=True)
    else:
        w = None
    return quantile(D.flatten(), q=q, weights=w) ** 2

bw_quantile_1d(D: torch.Tensor, q: float, weights: torch.Tensor | None = None) -> torch.Tensor

Calculate bandwidth using the specified quantile per row and column.

PARAMETER	DESCRIPTION
D	Pairwise distance matrix. TYPE: Tensor
q	Quantile value in [0.0, 1.0]. TYPE: float
weights	Sample weights, shape (n,). TYPE: torch.Tensor | None, optional, by default None DEFAULT: None

RETURNS	DESCRIPTION
Tensor	Bandwidth matrix as outer product of per-row and per-column quantiles.

Source code in spectre/kernel/bw.py

def bw_quantile_1d(
    D: torch.Tensor, q: float, weights: torch.Tensor | None = None
) -> torch.Tensor:
    """
    Calculate bandwidth using the specified quantile per row and column.

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix.

    q : float
        Quantile value in [0.0, 1.0].

    weights : torch.Tensor | None, optional, by default None
        Sample weights, shape (n,).

    Returns
    -------
    torch.Tensor
        Bandwidth matrix as outer product of per-row and per-column quantiles.
    """
    if weights is not None:
        w_rows = weights.unsqueeze(0).expand_as(D)
        w_cols = weights.unsqueeze(1).expand_as(D)
    else:
        w_rows = None
        w_cols = None

    bw_rows = quantile(D, q=q, weights=w_rows, dim=1)
    bw_cols = quantile(D, q=q, weights=w_cols, dim=0)

    return outer(bw_rows, bw_cols)

bw_quantile_search(D: torch.Tensor, n_quantiles: int = 100, eps: float = 0.001, weights: torch.Tensor | None = None) -> torch.Tensor

Calculate bandwidth using quantile search method.

Finds the optimal bandwidth by searching across quantiles of the distance distribution and selecting the one with maximum derivative.

PARAMETER	DESCRIPTION
D	Pairwise distance matrix. TYPE: Tensor
n_quantiles	Number of quantiles to evaluate. TYPE: int, optional, by default 100 DEFAULT: 100
eps	Epsilon value to exclude symmetrically from the quantiles. TYPE: float, optional, by default 1e-3 DEFAULT: 0.001
weights	Sample weights, shape (n,). TYPE: torch.Tensor | None, optional, by default None DEFAULT: None

RETURNS	DESCRIPTION
Tensor	Optimal bandwidth value.

Source code in spectre/kernel/bw.py

def bw_quantile_search(
    D: torch.Tensor,
    n_quantiles: int = 100,
    eps: float = 1e-3,
    weights: torch.Tensor | None = None,
) -> torch.Tensor:
    """
    Calculate bandwidth using quantile search method.

    Finds the optimal bandwidth by searching across quantiles of the distance
    distribution and selecting the one with maximum derivative.

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix.

    n_quantiles : int, optional, by default 100
        Number of quantiles to evaluate.

    eps : float, optional, by default 1e-3
        Epsilon value to exclude symmetrically from the quantiles.

    weights : torch.Tensor | None, optional, by default None
        Sample weights, shape (n,).

    Returns
    -------
    torch.Tensor
        Optimal bandwidth value.
    """
    if weights is not None:
        w = outer(weights, flatten=True)
    else:
        w = None

    quantiles_q = torch.linspace(eps, 1 - eps, n_quantiles)
    epsilon = quantile(D.flatten(), q=quantiles_q, weights=w)
    deriv = torch.diff(epsilon) / torch.diff(quantiles_q)
    idx = torch.argmax(deriv)

    return epsilon[idx] ** 2

bw_k_nearest_neighbors(D: torch.Tensor, k_neighbor: int, weights: torch.Tensor | None = None) -> torch.Tensor

Calculate bandwidth matrix for Gaussian kernel from k nearest neighbors.

Computes k-nearest neighbors along both dimensions and returns their outer product, producing a bandwidth matrix of the same shape as the input distance matrix.

PARAMETER	DESCRIPTION
D	Pairwise distance matrix of shape (n_samples, m_samples). TYPE: Tensor
k_neighbor	Number of nearest neighbors. Must be < min(n_samples, m_samples). TYPE: int

RETURNS	DESCRIPTION
Tensor	Bandwidth matrix of shape (n_samples, m_samples) computed as outer product of per-row and per-column k-nearest neighbor distances.

Source code in spectre/kernel/bw.py

def bw_k_nearest_neighbors(
    D: torch.Tensor,
    k_neighbor: int,
    weights: torch.Tensor | None = None,
) -> torch.Tensor:
    """
    Calculate bandwidth matrix for Gaussian kernel from k nearest neighbors.

    Computes k-nearest neighbors along both dimensions and returns their outer
    product, producing a bandwidth matrix of the same shape as the input
    distance matrix.

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix of shape (n_samples, m_samples).

    k_neighbor : int
        Number of nearest neighbors. Must be < min(n_samples, m_samples).

    Returns
    -------
    torch.Tensor
        Bandwidth matrix of shape (n_samples, m_samples) computed as outer
        product of per-row and per-column k-nearest neighbor distances.
    """
    min_dim = min(D.shape[0], D.shape[1])
    check_in_interval(k_neighbor, "[1, inf)")
    check_in_interval(k_neighbor, f"[1, {min_dim})")

    # k-nearest neighbors along rows (for each row, find k nearest in columns)
    bw_rows = D.topk(k_neighbor + 1, dim=1, largest=False)[0][:, -1]  # shape: (n, )

    # k-nearest neighbors along columns (for each column, find k nearest in rows)
    bw_cols = D.topk(k_neighbor + 1, dim=0, largest=False)[0][-1, :]  # shape: (m, )

    return outer(bw_rows, bw_cols)

bw_bgh(D: torch.Tensor, bw_values: torch.Tensor | None = None) -> torch.Tensor

Calculate the bandwidth for a Gaussian kernel matrix using the BGH method [1].

PARAMETER	DESCRIPTION
D	Pairwise distance matrix. TYPE: Tensor
bw_values	Bandwidth values to test. TYPE: torch.Tensor | None, optional, by default None DEFAULT: None

RETURNS	DESCRIPTION
Tensor	Optimal bandwidth value.

Source code in spectre/kernel/bw.py

def bw_bgh(D: torch.Tensor, bw_values: torch.Tensor | None = None) -> torch.Tensor:
    """
    Calculate the bandwidth for a Gaussian kernel matrix using the BGH method
    [@berry2015nonparametric].

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix.

    bw_values : torch.Tensor | None, optional, by default None
        Bandwidth values to test.

    Returns
    -------
    torch.Tensor
        Optimal bandwidth value.
    """
    if bw_values is None:
        bw_values = 2 ** torch.arange(-10.0, 10.0, 1.0)
    bw_values = torch.sort(bw_values).values

    log_T = torch.stack(
        [torch.logsumexp(-(D.flatten() ** 2) / bw**2, dim=-1) for bw in bw_values]
    )
    log_eps = torch.log(bw_values)
    log_deriv = torch.diff(log_T) / torch.diff(log_eps)
    max_loc = torch.argmax(log_deriv)
    epsilon = torch.exp(log_eps[max_loc])

    return epsilon

bw_median_scaled(D: torch.Tensor, d: int, alpha: float = 0.5, weights: torch.Tensor | None = None) -> torch.Tensor

Calculate bandwidth using median of distances with dimension scaling.

Applies dimension-aware scaling to account for curse of dimensionality. Useful for comparing kernel densities across feature subsets with different dimensionalities.

PARAMETER	DESCRIPTION
D	Pairwise distance matrix. TYPE: Tensor
d	Feature dimensionality (number of features in original data). TYPE: int
alpha	Scaling exponent. Bandwidth is multiplied by d^alpha. Typical values: 0.2-1.0. Higher values apply stronger penalty for high dimensions. TYPE: float, optional, by default 0.5 DEFAULT: 0.5
weights	Sample weights, shape (n,). TYPE: torch.Tensor | None, optional, by default None DEFAULT: None

RETURNS	DESCRIPTION
Tensor	Squared bandwidth value scaled by dimensionality.

Source code in spectre/kernel/bw.py

def bw_median_scaled(
    D: torch.Tensor,
    d: int,
    alpha: float = 0.5,
    weights: torch.Tensor | None = None,
) -> torch.Tensor:
    """
    Calculate bandwidth using median of distances with dimension scaling.

    Applies dimension-aware scaling to account for curse of dimensionality.
    Useful for comparing kernel densities across feature subsets with different
    dimensionalities.

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix.

    d : int
        Feature dimensionality (number of features in original data).

    alpha : float, optional, by default 0.5
        Scaling exponent. Bandwidth is multiplied by d^alpha. Typical values:
        0.2-1.0. Higher values apply stronger penalty for high dimensions.

    weights : torch.Tensor | None, optional, by default None
        Sample weights, shape (n,).

    Returns
    -------
    torch.Tensor
        Squared bandwidth value scaled by dimensionality.
    """
    return bw_median(D, weights=weights) * d**alpha

bw_median_1d_scaled(D: torch.Tensor, d: int, alpha: float = 0.5, weights: torch.Tensor | None = None) -> torch.Tensor

Calculate bandwidth matrix using per-dimension median with dimension scaling.

PARAMETER	DESCRIPTION
D	Pairwise distance matrix. TYPE: Tensor
d	Feature dimensionality (number of features in original data). TYPE: int
alpha	Scaling exponent. Bandwidth is multiplied by d^alpha. TYPE: float, optional, by default 0.5 DEFAULT: 0.5
weights	Sample weights, shape (n,). TYPE: torch.Tensor | None, optional, by default None DEFAULT: None

RETURNS	DESCRIPTION
Tensor	Bandwidth matrix as outer product scaled by dimensionality.

Source code in spectre/kernel/bw.py

def bw_median_1d_scaled(
    D: torch.Tensor,
    d: int,
    alpha: float = 0.5,
    weights: torch.Tensor | None = None,
) -> torch.Tensor:
    """
    Calculate bandwidth matrix using per-dimension median with dimension
    scaling.

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix.

    d : int
        Feature dimensionality (number of features in original data).

    alpha : float, optional, by default 0.5
        Scaling exponent. Bandwidth is multiplied by d^alpha.

    weights : torch.Tensor | None, optional, by default None
        Sample weights, shape (n,).

    Returns
    -------
    torch.Tensor
        Bandwidth matrix as outer product scaled by dimensionality.
    """
    return bw_median_1d(D, weights=weights) * d**alpha

bw_silverman(D: torch.Tensor, d: int, n: int | float | None = None, weights: torch.Tensor | None = None) -> torch.Tensor

Calculate bandwidth using Silverman's rule of thumb.

Silverman's rule provides automatic bandwidth selection that accounts for both sample size and dimensionality. Based on optimal bandwidth for Gaussian densities.

When weights are provided, the spread estimate uses the weighted median and the sample size is replaced by Kish's effective sample size: n_eff = (sum w_i)^2 / sum(w_i^2).

PARAMETER	DESCRIPTION
D	Pairwise distance matrix. TYPE: Tensor
d	Feature dimensionality (number of features in original data). TYPE: int
n	Sample size. If None, inferred from D.shape[0] (or Kish's effective sample size when weights are provided). TYPE: int | float | None, optional, by default None DEFAULT: None
weights	Sample weights, shape (n_samples,). TYPE: torch.Tensor | None, optional, by default None DEFAULT: None

RETURNS	DESCRIPTION
Tensor	Squared bandwidth value from Silverman's rule.

Source code in spectre/kernel/bw.py

def bw_silverman(
    D: torch.Tensor,
    d: int,
    n: int | float | None = None,
    weights: torch.Tensor | None = None,
) -> torch.Tensor:
    """
    Calculate bandwidth using Silverman's rule of thumb.

    Silverman's rule provides automatic bandwidth selection that accounts for
    both sample size and dimensionality. Based on optimal bandwidth for
    Gaussian densities.

    When `weights` are provided, the spread estimate uses the weighted median
    and the sample size is replaced by Kish's effective sample size: `n_eff =
    (sum w_i)^2 / sum(w_i^2)`.

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix.

    d : int
        Feature dimensionality (number of features in original data).

    n : int | float | None, optional, by default None
        Sample size. If None, inferred from D.shape[0] (or Kish's effective
        sample size when `weights` are provided).

    weights : torch.Tensor | None, optional, by default None
        Sample weights, shape (n_samples,).

    Returns
    -------
    torch.Tensor
        Squared bandwidth value from Silverman's rule.
    """
    median = bw_median(D, weights=weights)

    if n is None:
        if weights is not None:
            # Kish's effective sample size: n_eff = (sum w_i)^2 / sum(w_i^2)
            n = float(weights.sum() ** 2 / (weights**2).sum())
        else:
            n = D.shape[0]

    # Estimate standard deviation from median pairwise distance
    # For normal data, median(|X-Y|) ≈ sqrt(2) * sigma * 0.6745
    # Therefore: sigma ≈ median(D) / (sqrt(2) * 0.6745)
    sigma = median / (1.4142 * 0.6745)

    # Silverman's rule: h = (n*(d+2)/4)^(-1/(d+4)) * sigma
    h = ((n * (d + 2) / 4) ** (-1.0 / (d + 4))) * sigma

    return h**2

bw_scott(D: torch.Tensor, d: int, n: int | float | None = None, weights: torch.Tensor | None = None) -> torch.Tensor

Calculate bandwidth using Scott's rule.

Scott's rule provides automatic bandwidth selection similar to Silverman's but with different constants. Often produces slightly larger bandwidths.

When weights are provided, the spread estimate uses the weighted median and the sample size is replaced by Kish's effective sample size: n_eff = (sum w_i)^2 / sum(w_i^2).

PARAMETER	DESCRIPTION
D	Pairwise distance matrix. TYPE: Tensor
d	Feature dimensionality (number of features in original data). TYPE: int
n	Sample size. If None, inferred from D.shape[0] (or Kish's effective sample size when weights are provided). TYPE: int | float | None, optional, by default None DEFAULT: None
weights	Sample weights, shape (n_samples,). TYPE: torch.Tensor | None, optional, by default None DEFAULT: None

RETURNS	DESCRIPTION
Tensor	Squared bandwidth value from Scott's rule.

Source code in spectre/kernel/bw.py

def bw_scott(
    D: torch.Tensor,
    d: int,
    n: int | float | None = None,
    weights: torch.Tensor | None = None,
) -> torch.Tensor:
    """
    Calculate bandwidth using Scott's rule.

    Scott's rule provides automatic bandwidth selection similar to Silverman's
    but with different constants. Often produces slightly larger bandwidths.

    When `weights` are provided, the spread estimate uses the weighted median
    and the sample size is replaced by Kish's effective sample size: `n_eff =
    (sum w_i)^2 / sum(w_i^2)`.

    Parameters
    ----------
    D : torch.Tensor
        Pairwise distance matrix.

    d : int
        Feature dimensionality (number of features in original data).

    n : int | float | None, optional, by default None
        Sample size. If None, inferred from D.shape[0] (or Kish's effective
        sample size when `weights` are provided).

    weights : torch.Tensor | None, optional, by default None
        Sample weights, shape (n_samples,).

    Returns
    -------
    torch.Tensor
        Squared bandwidth value from Scott's rule.
    """
    median = bw_median(D, weights=weights)

    if n is None:
        if weights is not None:
            # Kish's effective sample size: n_eff = (sum w_i)^2 / sum(w_i^2)
            n = float(weights.sum() ** 2 / (weights**2).sum())
        else:
            n = D.shape[0]

    # Estimate standard deviation from median pairwise distance
    # For normal data, median(|X-Y|) ≈ sqrt(2) * sigma * 0.6745
    # Therefore: sigma ≈ median(D) / (sqrt(2) * 0.6745)
    sigma = median / (1.4142 * 0.6745)

    # Scott's rule: h = n^(-1/(d+4)) * sigma
    h = (n ** (-1.0 / (d + 4))) * sigma

    return h**2