v0.29.12

vectorized kernel

Author: bartzbeielstein

.vscode/settings.json

@@ -43,5 +43,10 @@
     ],
     "cSpell.enableFiletypes": [
         "quarto"
-    ]
+    ],
+    "python.testing.pytestArgs": [
+        "test"
+    ],
+    "python.testing.unittestEnabled": false,
+    "python.testing.pytestEnabled": true
 }

pyproject.toml

@@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "spotpython"
-version = "0.29.11"
+version = "0.29.12"
 authors = [
   { name="T. Bartz-Beielstein", email="tbb@bartzundbartz.de" }
 ]

src/spotpython/surrogate/kriging.py

@@ -316,6 +316,34 @@ def predict(self, X: np.ndarray, return_std=False, return_val: str = "y") -> np.
             predictions = [self._pred(x_i)[0] for x_i in X]
             return np.array(predictions)
 
+    def _kernel(self, X: np.ndarray, theta: np.ndarray, p: float) -> np.ndarray:
+        """
+        Computes the correlation matrix Psi using vectorized operations.
+
+        Args:
+            X (np.ndarray): Input data of shape (n_samples, n_features).
+            theta (np.ndarray): Theta parameters of shape (n_features,).
+            p (float): Power exponent.
+
+        Returns:
+            np.ndarray: The upper triangle of the correlation matrix Psi.
+        """
+        n_samples, n_features = X.shape
+        Psi = np.zeros((n_samples, n_samples), dtype=float)
+        # Calculate all pairwise differences:
+        # X_expanded_rows will have shape (n_samples, 1, n_features)
+        # X_expanded_cols will have shape (1, n_samples, n_features)
+        # diff will have shape (n_samples, n_samples, n_features)
+        diff = np.abs(X[:, np.newaxis, :] - X[np.newaxis, :, :]) ** p
+        # Apply theta and sum over features
+        # dist_matrix will have shape (n_samples, n_samples)
+        dist_matrix = np.sum(theta * diff, axis=2)
+        # Compute Psi using the exponential kernel
+        Psi = np.exp(-dist_matrix)
+        # Return only the upper triangle, as the matrix is symmetric
+        # and the diagonal will be handled later.
+        return np.triu(Psi, k=1)
+
     def likelihood(self, x: np.ndarray) -> Tuple[float, np.ndarray, np.ndarray]:
         """
         Computes the negative of the concentrated log-likelihood for a given set
@@ -359,13 +387,8 @@ def likelihood(self, x: np.ndarray) -> Tuple[float, np.ndarray, np.ndarray]:
         one = np.ones(n)
 
         # Build correlation matrix
-        Psi = np.zeros((n, n), dtype=float)
-        for i in range(n):
-            for j in range(i + 1, n):
-                dist_vec = np.abs(X[i, :] - X[j, :]) ** p
-                Psi[i, j] = np.exp(-np.sum(theta * dist_vec))
-
-        Psi = Psi + Psi.T + np.eye(n) + np.eye(n) * lambda_
+        Psi_upper_triangle = self._kernel(X, theta, p)
+        Psi = Psi_upper_triangle + Psi_upper_triangle.T + np.eye(n) + np.eye(n) * lambda_
 
         try:
             U = np.linalg.cholesky(Psi)

test/test_kernel.py

@@ -0,0 +1,124 @@
+import numpy as np
+import pytest
+from spotpython.surrogate.kriging import Kriging # Assuming Kriging class is in this path
+
+# Create a dummy Kriging instance for testing, as _kernel is a method of the class
+# We can use default parameters for the Kriging constructor as they don't affect _kernel
+@pytest.fixture
+def kriging_model():
+    """Fixture to create a Kriging model instance."""
+    return Kriging()
+
+def test_kernel_simple_1d(kriging_model):
+    """Test _kernel with 2 samples, 1 feature."""
+    X = np.array([[0.0], [1.0]])
+    theta = np.array([1.0])
+    p = 2.0
+    
+    expected_psi_upper_triangle = np.array([
+        [0.0, np.exp(-1.0)],
+        [0.0, 0.0]
+    ])
+    
+    result = kriging_model._kernel(X, theta, p)
+    np.testing.assert_array_almost_equal(result, expected_psi_upper_triangle)
+
+def test_kernel_simple_2d(kriging_model):
+    """Test _kernel with 2 samples, 2 features."""
+    X = np.array([[0.0, 0.0], [1.0, 1.0]])
+    theta = np.array([1.0, 1.0])
+    p = 2.0
+    
+    # diff for [0,1]: [abs(0-1)^2, abs(0-1)^2] = [1,1]
+    # dist_matrix[0,1] = sum(theta * diff) = 1*1 + 1*1 = 2
+    # Psi[0,1] = exp(-2)
+    expected_psi_upper_triangle = np.array([
+        [0.0, np.exp(-2.0)],
+        [0.0, 0.0]
+    ])
+    
+    result = kriging_model._kernel(X, theta, p)
+    np.testing.assert_array_almost_equal(result, expected_psi_upper_triangle)
+
+def test_kernel_three_samples_1d(kriging_model):
+    """Test _kernel with 3 samples, 1 feature."""
+    X = np.array([[0.0], [1.0], [2.0]])
+    theta = np.array([0.5])
+    p = 2.0
+    
+    # For X[0] and X[1]: diff = abs(0-1)^2 = 1. dist = 0.5 * 1 = 0.5. val = exp(-0.5)
+    # For X[0] and X[2]: diff = abs(0-2)^2 = 4. dist = 0.5 * 4 = 2.0. val = exp(-2.0)
+    # For X[1] and X[2]: diff = abs(1-2)^2 = 1. dist = 0.5 * 1 = 0.5. val = exp(-0.5)
+    expected_psi_upper_triangle = np.array([
+        [0.0, np.exp(-0.5), np.exp(-2.0)],
+        [0.0, 0.0,         np.exp(-0.5)],
+        [0.0, 0.0,         0.0]
+    ])
+    
+    result = kriging_model._kernel(X, theta, p)
+    np.testing.assert_array_almost_equal(result, expected_psi_upper_triangle)
+
+def test_kernel_different_p_value(kriging_model):
+    """Test _kernel with a different p value."""
+    X = np.array([[0.0], [1.0]])
+    theta = np.array([1.0])
+    p = 1.0 # Changed p value
+    
+    # diff for [0,1]: abs(0-1)^1 = 1
+    # dist_matrix[0,1] = sum(theta * diff) = 1*1 = 1
+    # Psi[0,1] = exp(-1)
+    expected_psi_upper_triangle = np.array([
+        [0.0, np.exp(-1.0)],
+        [0.0, 0.0]
+    ])
+    
+    result = kriging_model._kernel(X, theta, p)
+    np.testing.assert_array_almost_equal(result, expected_psi_upper_triangle)
+
+def test_kernel_theta_zeros(kriging_model):
+    """Test _kernel when a theta value is zero."""
+    X = np.array([[0.0, 0.0], [1.0, 1.0]])
+    theta = np.array([1.0, 0.0]) # Second feature's theta is 0
+    p = 2.0
+    
+    # diff for [0,1]: [abs(0-1)^2, abs(0-1)^2] = [1,1]
+    # dist_matrix[0,1] = sum(theta * diff) = 1*1 + 0*1 = 1
+    # Psi[0,1] = exp(-1)
+    expected_psi_upper_triangle = np.array([
+        [0.0, np.exp(-1.0)],
+        [0.0, 0.0]
+    ])
+    
+    result = kriging_model._kernel(X, theta, p)
+    np.testing.assert_array_almost_equal(result, expected_psi_upper_triangle)
+
+def test_kernel_no_samples(kriging_model):
+    """Test _kernel with no samples. Should produce an empty array or handle gracefully."""
+    X = np.empty((0, 2)) # No samples, 2 features
+    theta = np.array([1.0, 1.0])
+    p = 2.0
+    
+    # Expected: an empty (0,0) array for the upper triangle
+    expected_psi_upper_triangle = np.empty((0,0)) 
+    
+    result = kriging_model._kernel(X, theta, p)
+    # The current implementation of _kernel will likely raise an error or behave unexpectedly
+    # with X.shape[0] == 0 before np.triu.
+    # np.zeros((0,0)) is valid.
+    # np.abs(X[:, np.newaxis, :] - X[np.newaxis, :, :]) will result in (0,0,2) shape
+    # np.sum on axis 2 will result in (0,0) shape
+    # np.exp will result in (0,0) shape
+    # np.triu will result in (0,0) shape
+    np.testing.assert_array_almost_equal(result, expected_psi_upper_triangle)
+
+def test_kernel_one_sample(kriging_model):
+    """Test _kernel with only one sample."""
+    X = np.array([[1.0, 2.0]]) # One sample, 2 features
+    theta = np.array([1.0, 1.0])
+    p = 2.0
+    
+    # Expected: a (1,1) array with 0 in the upper triangle (as k=1)
+    expected_psi_upper_triangle = np.array([[0.0]])
+    
+    result = kriging_model._kernel(X, theta, p)
+    np.testing.assert_array_almost_equal(result, expected_psi_upper_triangle)

test/test_kriging.py test/test_kriging_1.pytest/test_kriging.py renamed to test/test_kriging_1.py

@@ -0,0 +1,70 @@
+import numpy as np
+import pytest
+from spotpython.surrogate.kriging import Kriging
+
+
+@pytest.mark.parametrize("method, x_shape", [
+    ("interpolation", (1,)),
+    ("regression", (2,)),
+    ("reinterpolation", (2,))
+])
+def test_likelihood_basic(method, x_shape):
+    # Small toy data: 3 points, 1D
+    X = np.array([[0.0], [0.5], [1.0]])
+    y = np.array([1.0, 2.0, 3.0])
+    model = Kriging(method=method, penalty=12345.0)
+    model.X_ = X
+    model.y_ = y
+    model.eps = 1e-6
+    model.penalty = 12345.0
+
+    # log10(theta) = 0, log10(lambda) = -3 for regression/reinterpolation
+    if method == "interpolation":
+        x = np.zeros(x_shape)
+    else:
+        x = np.zeros(x_shape)
+        x[-1] = -3  # log10(lambda)
+
+    negLnLike, Psi, U = model.likelihood(x)
+
+    # Check types and shapes
+    assert isinstance(negLnLike, float)
+    assert Psi.shape == (3, 3)
+    assert np.allclose(Psi, Psi.T)
+    if U is not None:
+        assert U.shape == (3, 3)
+        # U should be upper-triangular (Cholesky returns lower, but code uses U as upper)
+        # Actually, np.linalg.cholesky returns lower-triangular, so U @ U.T == Psi
+        assert np.allclose(U @ U.T, Psi, atol=1e-8)
+    else:
+        # If U is None, penalty should be returned
+        assert negLnLike == model.penalty
+
+    # negLnLike should be finite unless penalty
+    if U is not None:
+        assert np.isfinite(negLnLike)
+
+def test_likelihood_invalid_method():
+    X = np.array([[0.0], [1.0]])
+    y = np.array([1.0, 2.0])
+    model = Kriging(method="regression")
+    model.X_ = X
+    model.y_ = y
+    # Set an invalid method
+    model.method = "invalid"
+    with pytest.raises(ValueError):
+        model.likelihood(np.zeros(1))
+
+def test_likelihood_ill_conditioned_returns_penalty():
+    # Use two identical points to force singular Psi
+    X = np.array([[0.0], [0.0]])
+    y = np.array([1.0, 1.0])
+    model = Kriging(method="interpolation", penalty=99999.0)
+    model.X_ = X
+    model.y_ = y
+    model.eps = 0.0  # No regularization
+    x = np.zeros(1)
+    negLnLike, Psi, U = model.likelihood(x)
+    assert negLnLike == 99999.0
+    assert U is None
+    assert Psi.shape == (2, 2)