Derivative-Constrained GPR

Exact RBF Gaussian Process Regression with function-value and derivative observations.

matgpr.derivative_gpr

DerivativeObservationSet dataclass

Observed or physics-implied derivatives for derivative-constrained GPR.

Each row represents one derivative observation:

derivative_values[i] = df(X[i]) / dX[feature_indices[i]].

Derivatives must be expressed with respect to the same feature scale used in the GPR feature matrix. If the feature matrix is standardized, provide derivatives with respect to the standardized feature values.

Source code in matgpr/derivative_gpr.py

@dataclass(frozen=True)
class DerivativeObservationSet:
    """Observed or physics-implied derivatives for derivative-constrained GPR.

    Each row represents one derivative observation:

    ``derivative_values[i] = df(X[i]) / dX[feature_indices[i]]``.

    Derivatives must be expressed with respect to the same feature scale used
    in the GPR feature matrix. If the feature matrix is standardized, provide
    derivatives with respect to the standardized feature values.
    """

    X: np.ndarray
    feature_indices: int | Sequence[int] | np.ndarray
    derivative_values: Sequence[float] | np.ndarray
    noise_std: float | Sequence[float] | np.ndarray = 0.1
    labels: Sequence[str] | np.ndarray | None = None

    def __post_init__(self) -> None:
        x_values = to_2d_finite(self.X, "X")
        derivative_values = to_1d_finite(self.derivative_values, "derivative_values")
        validate_same_row_count(x_values, derivative_values, "X", "derivative_values")

        feature_indices = _feature_index_array(
            self.feature_indices,
            n_observations=derivative_values.shape[0],
            n_features=x_values.shape[1],
        )
        noise = noise_std_array(self.noise_std, derivative_values.shape[0])
        labels = label_array(self.labels, derivative_values.shape[0], default="derivative")

        object.__setattr__(self, "X", x_values)
        object.__setattr__(self, "feature_indices", feature_indices)
        object.__setattr__(self, "derivative_values", derivative_values)
        object.__setattr__(self, "noise_std", noise)
        object.__setattr__(self, "labels", labels)

    @property
    def n_observations(self) -> int:
        """Number of derivative observations."""
        return int(self.X.shape[0])

    @property
    def n_features(self) -> int:
        """Number of feature columns."""
        return int(self.X.shape[1])

    @property
    def alpha(self) -> np.ndarray:
        """Derivative observation variances."""
        return self.noise_std**2

n_observations property

Number of derivative observations.

n_features property

Number of feature columns.

alpha property

Derivative observation variances.

MonotonicDerivativeConstraint dataclass

Create derivative observations for a local monotonic trend.

For an increasing trend, generated derivative observations have positive target value minimum_slope. For a decreasing trend, generated derivative observations have target value -minimum_slope.

These are soft equality observations of the derivative. They encourage the posterior to follow the trend near the anchor points, but they are not a formal guarantee that the function is globally monotonic everywhere.

Source code in matgpr/derivative_gpr.py

@dataclass(frozen=True)
class MonotonicDerivativeConstraint:
    """Create derivative observations for a local monotonic trend.

    For an increasing trend, generated derivative observations have positive
    target value ``minimum_slope``. For a decreasing trend, generated
    derivative observations have target value ``-minimum_slope``.

    These are soft equality observations of the derivative. They encourage the
    posterior to follow the trend near the anchor points, but they are not a
    formal guarantee that the function is globally monotonic everywhere.
    """

    feature: int | str
    direction: str
    minimum_slope: float = 0.0
    noise_std: float | Sequence[float] | np.ndarray = 0.1
    feature_min: float | None = None
    feature_max: float | None = None
    label: str | None = None

    def generate(self, X_reference) -> DerivativeObservationSet:
        """Generate derivative observations from reference feature rows."""
        x_values, columns = _as_numeric_matrix(X_reference, "X_reference")
        feature_index = _resolve_feature_index(self.feature, columns, x_values.shape[1])
        direction_sign = response_sign(self.direction)
        minimum_slope = nonnegative_scalar(self.minimum_slope, "minimum_slope")
        derivative_value = direction_sign * minimum_slope

        keep_mask = feature_bounds_mask(
            x_values[:, feature_index],
            feature_min=self.feature_min,
            feature_max=self.feature_max,
        )
        if not np.any(keep_mask):
            raise ValueError("No derivative observations remain after applying feature bounds")

        x_virtual = x_values[keep_mask]
        derivative_values = np.full(x_virtual.shape[0], derivative_value, dtype=float)
        feature_indices = np.full(x_virtual.shape[0], feature_index, dtype=int)
        noise = noise_std_array(self.noise_std, x_values.shape[0])[keep_mask]
        label = self.label or f"derivative_{normalized_direction(self.direction)}"
        labels = np.full(x_virtual.shape[0], label, dtype=object)

        return DerivativeObservationSet(
            X=x_virtual,
            feature_indices=feature_indices,
            derivative_values=derivative_values,
            noise_std=noise,
            labels=labels,
        )

generate(X_reference)

Generate derivative observations from reference feature rows.

Source code in matgpr/derivative_gpr.py

def generate(self, X_reference) -> DerivativeObservationSet:
    """Generate derivative observations from reference feature rows."""
    x_values, columns = _as_numeric_matrix(X_reference, "X_reference")
    feature_index = _resolve_feature_index(self.feature, columns, x_values.shape[1])
    direction_sign = response_sign(self.direction)
    minimum_slope = nonnegative_scalar(self.minimum_slope, "minimum_slope")
    derivative_value = direction_sign * minimum_slope

    keep_mask = feature_bounds_mask(
        x_values[:, feature_index],
        feature_min=self.feature_min,
        feature_max=self.feature_max,
    )
    if not np.any(keep_mask):
        raise ValueError("No derivative observations remain after applying feature bounds")

    x_virtual = x_values[keep_mask]
    derivative_values = np.full(x_virtual.shape[0], derivative_value, dtype=float)
    feature_indices = np.full(x_virtual.shape[0], feature_index, dtype=int)
    noise = noise_std_array(self.noise_std, x_values.shape[0])[keep_mask]
    label = self.label or f"derivative_{normalized_direction(self.direction)}"
    labels = np.full(x_virtual.shape[0], label, dtype=object)

    return DerivativeObservationSet(
        X=x_virtual,
        feature_indices=feature_indices,
        derivative_values=derivative_values,
        noise_std=noise,
        labels=labels,
    )

DerivativeConstrainedGPRPrediction dataclass

Predictions from a derivative-constrained exact GPR model.

Source code in matgpr/derivative_gpr.py

@dataclass(frozen=True)
class DerivativeConstrainedGPRPrediction:
    """Predictions from a derivative-constrained exact GPR model."""

    mean: np.ndarray
    std: np.ndarray | None = None
    lower: np.ndarray | None = None
    upper: np.ndarray | None = None

DerivativeConstrainedGPRResult dataclass

Fitted exact GPR model with derivative observations.

Source code in matgpr/derivative_gpr.py

@dataclass(frozen=True)
class DerivativeConstrainedGPRResult:
    """Fitted exact GPR model with derivative observations."""

    X_train: np.ndarray
    y_train: np.ndarray
    derivative_observations: DerivativeObservationSet | None
    length_scale: np.ndarray
    signal_variance: float
    value_noise_std: np.ndarray
    standardize_y: bool
    target_mean: float
    target_scale: float
    lower_cholesky: np.ndarray
    alpha: np.ndarray
    log_marginal_likelihood: float
    optimizer_success: bool | None = None
    optimizer_message: str | None = None

    def predict(
        self,
        X,
        *,
        return_std: bool = True,
        confidence_level: float | None = None,
        include_observation_noise: bool = False,
    ) -> DerivativeConstrainedGPRPrediction:
        """Predict function values at new feature rows."""
        validate_confidence_level(confidence_level)
        x_test, _ = _as_numeric_matrix(X, "X")
        if x_test.shape[1] != self.X_train.shape[1]:
            raise ValueError(
                f"X has {x_test.shape[1]} features, but the model was fitted with "
                f"{self.X_train.shape[1]} features"
            )

        cross_covariance = _cross_covariance_to_training(
            x_test,
            self.X_train,
            self.derivative_observations,
            length_scale=self.length_scale,
            signal_variance=self.signal_variance,
        )
        mean_model = cross_covariance @ self.alpha
        mean = mean_model * self.target_scale + self.target_mean

        std = None
        if return_std or confidence_level is not None:
            solved = _solve_lower_triangular(self.lower_cholesky, cross_covariance.T)
            variance_model = self.signal_variance - np.sum(solved * solved, axis=0)
            if include_observation_noise:
                variance_model = variance_model + float(
                    np.mean((self.value_noise_std / self.target_scale) ** 2)
                )
            variance_model = np.maximum(variance_model, 0.0)
            std = np.sqrt(variance_model) * self.target_scale

        lower = None
        upper = None
        if confidence_level is not None:
            if std is None:
                raise ValueError("confidence_level requires return_std=True")
            z_value = NormalDist().inv_cdf(0.5 + confidence_level / 2.0)
            lower = mean - z_value * std
            upper = mean + z_value * std

        return DerivativeConstrainedGPRPrediction(mean=mean, std=std, lower=lower, upper=upper)

predict(X, *, return_std=True, confidence_level=None, include_observation_noise=False)

Predict function values at new feature rows.

Source code in matgpr/derivative_gpr.py

def predict(
    self,
    X,
    *,
    return_std: bool = True,
    confidence_level: float | None = None,
    include_observation_noise: bool = False,
) -> DerivativeConstrainedGPRPrediction:
    """Predict function values at new feature rows."""
    validate_confidence_level(confidence_level)
    x_test, _ = _as_numeric_matrix(X, "X")
    if x_test.shape[1] != self.X_train.shape[1]:
        raise ValueError(
            f"X has {x_test.shape[1]} features, but the model was fitted with "
            f"{self.X_train.shape[1]} features"
        )

    cross_covariance = _cross_covariance_to_training(
        x_test,
        self.X_train,
        self.derivative_observations,
        length_scale=self.length_scale,
        signal_variance=self.signal_variance,
    )
    mean_model = cross_covariance @ self.alpha
    mean = mean_model * self.target_scale + self.target_mean

    std = None
    if return_std or confidence_level is not None:
        solved = _solve_lower_triangular(self.lower_cholesky, cross_covariance.T)
        variance_model = self.signal_variance - np.sum(solved * solved, axis=0)
        if include_observation_noise:
            variance_model = variance_model + float(
                np.mean((self.value_noise_std / self.target_scale) ** 2)
            )
        variance_model = np.maximum(variance_model, 0.0)
        std = np.sqrt(variance_model) * self.target_scale

    lower = None
    upper = None
    if confidence_level is not None:
        if std is None:
            raise ValueError("confidence_level requires return_std=True")
        z_value = NormalDist().inv_cdf(0.5 + confidence_level / 2.0)
        lower = mean - z_value * std
        upper = mean + z_value * std

    return DerivativeConstrainedGPRPrediction(mean=mean, std=std, lower=lower, upper=upper)

combine_derivative_observations(*derivative_observations)

Combine compatible derivative-observation sets.

Source code in matgpr/derivative_gpr.py

def combine_derivative_observations(
    *derivative_observations: DerivativeObservationSet,
) -> DerivativeObservationSet:
    """Combine compatible derivative-observation sets."""
    sets = _flatten_derivative_observation_sets(derivative_observations)
    if not sets:
        raise ValueError("At least one DerivativeObservationSet is required")

    n_features = sets[0].n_features
    for observation_set in sets:
        if observation_set.n_features != n_features:
            raise ValueError(
                "All derivative observation sets must have the same number of features"
            )

    return DerivativeObservationSet(
        X=np.vstack([observation_set.X for observation_set in sets]),
        feature_indices=np.concatenate(
            [observation_set.feature_indices for observation_set in sets]
        ),
        derivative_values=np.concatenate(
            [observation_set.derivative_values for observation_set in sets]
        ),
        noise_std=np.concatenate([observation_set.noise_std for observation_set in sets]),
        labels=np.concatenate([observation_set.labels for observation_set in sets]),
    )

fit_derivative_constrained_gpr(X_train, y_train, derivative_observations=None, *, length_scale=None, signal_variance=1.0, value_noise_std=1e-06, standardize_y=True, optimize_hyperparameters=False, length_scale_bounds=(0.0001, 10000.0), signal_variance_bounds=(1e-08, 100000000.0), maxiter=100, jitter=1e-08)

Fit exact RBF GPR with function and derivative observations.

The model uses the joint Gaussian-process covariance over function values and partial derivatives. This is useful when physics provides slope information, for example positive diffusivity change with temperature, decreasing viscosity with temperature, or a known zero-slope saturation region.

Source code in matgpr/derivative_gpr.py

def fit_derivative_constrained_gpr(
    X_train,
    y_train,
    derivative_observations: DerivativeObservationSet
    | Sequence[DerivativeObservationSet]
    | None = None,
    *,
    length_scale: float | Sequence[float] | np.ndarray | None = None,
    signal_variance: float | None = 1.0,
    value_noise_std: float | Sequence[float] | np.ndarray = 1e-6,
    standardize_y: bool = True,
    optimize_hyperparameters: bool = False,
    length_scale_bounds: tuple[float, float] = (1e-4, 1e4),
    signal_variance_bounds: tuple[float, float] = (1e-8, 1e8),
    maxiter: int = 100,
    jitter: float = 1e-8,
) -> DerivativeConstrainedGPRResult:
    """Fit exact RBF GPR with function and derivative observations.

    The model uses the joint Gaussian-process covariance over function values
    and partial derivatives. This is useful when physics provides slope
    information, for example positive diffusivity change with temperature,
    decreasing viscosity with temperature, or a known zero-slope saturation
    region.
    """
    x_train, _ = _as_numeric_matrix(X_train, "X_train")
    y_values = to_1d_finite(y_train, "y_train")
    validate_same_row_count(x_train, y_values, "X_train", "y_train")

    derivative_set = _prepare_derivative_observations(derivative_observations, x_train.shape[1])
    value_noise = noise_std_array(value_noise_std, y_values.shape[0])

    if standardize_y:
        target_mean = float(np.mean(y_values))
        target_scale = float(np.std(y_values, ddof=0))
        if target_scale <= 0:
            raise ValueError("Cannot standardize a constant target")
    else:
        target_mean = 0.0
        target_scale = 1.0

    y_model = (y_values - target_mean) / target_scale
    derivative_values_model = None
    derivative_noise_model = None
    if derivative_set is not None:
        derivative_values_model = derivative_set.derivative_values / target_scale
        derivative_noise_model = derivative_set.noise_std / target_scale
    value_noise_model = value_noise / target_scale

    length_scale_array = _initial_length_scale(length_scale, x_train)
    signal_variance_value = _initial_signal_variance(signal_variance, y_model)

    optimizer_success = None
    optimizer_message = None
    if optimize_hyperparameters:
        length_scale_array, signal_variance_value, optimizer_success, optimizer_message = (
            _optimize_hyperparameters(
                x_train,
                y_model,
                derivative_set,
                derivative_values_model,
                derivative_noise_model,
                value_noise_model,
                initial_length_scale=length_scale_array,
                initial_signal_variance=signal_variance_value,
                length_scale_bounds=length_scale_bounds,
                signal_variance_bounds=signal_variance_bounds,
                maxiter=maxiter,
                jitter=jitter,
            )
        )

    training_target = _combined_training_target(y_model, derivative_values_model)
    training_covariance = _training_covariance(
        x_train,
        derivative_set,
        length_scale=length_scale_array,
        signal_variance=signal_variance_value,
        value_noise_std=value_noise_model,
        derivative_noise_std=derivative_noise_model,
    )
    lower_cholesky = _stable_cholesky(training_covariance, jitter=jitter)
    alpha = _cho_solve(lower_cholesky, training_target)
    log_marginal_likelihood = _log_marginal_likelihood(
        training_target,
        lower_cholesky,
        alpha,
    )

    return DerivativeConstrainedGPRResult(
        X_train=x_train,
        y_train=y_values,
        derivative_observations=derivative_set,
        length_scale=length_scale_array,
        signal_variance=signal_variance_value,
        value_noise_std=value_noise,
        standardize_y=standardize_y,
        target_mean=target_mean,
        target_scale=target_scale,
        lower_cholesky=lower_cholesky,
        alpha=alpha,
        log_marginal_likelihood=log_marginal_likelihood,
        optimizer_success=optimizer_success,
        optimizer_message=optimizer_message,
    )