Fix for multidimensional gaussian process

CHANGELOGS.md

@@ -2,6 +2,8 @@
 
 ## 1.17.0 (development)
 
+* Fix for multidimensional gaussian process
+  [#1097](https://github.com/onnx/sklearn-onnx/pull/1097)
 * Minor fixes to support scikit-learn==1.5.0
   [#1095](https://github.com/onnx/sklearn-onnx/pull/1095)
 * Fix the conversion of pipeline including pipelines,

skl2onnx/operator_converters/gaussian_process.py

@@ -126,12 +126,12 @@ def convert_gaussian_process_regressor(
         if len(mean_y.shape) == 1:
             mean_y = mean_y.reshape(mean_y.shape + (1,))
 
-        if not hasattr(op, "_y_train_std") or op._y_train_std == 1:
+        if not hasattr(op, "_y_train_std") or np.all(op._y_train_std == 1):
             if isinstance(y_mean_b, (np.float32, np.float64)):
                 y_mean_b = np.array([y_mean_b])
             if isinstance(mean_y, (np.float32, np.float64)):
                 mean_y = np.array([mean_y])
-            y_mean = OnnxAdd(y_mean_b, mean_y, op_version=opv)
+            y_mean = OnnxAdd(y_mean_b, mean_y.T, op_version=opv)
         else:
             # A bug was fixed in 0.23 and it changed
             # the predictions when return_std is True.
@@ -145,13 +145,13 @@ def convert_gaussian_process_regressor(
             if isinstance(mean_y, (np.float32, np.float64)):
                 mean_y = np.array([mean_y])
             y_mean = OnnxAdd(
-                OnnxMul(y_mean_b, var_y, op_version=opv), mean_y, op_version=opv
+                OnnxMul(y_mean_b, var_y.T, op_version=opv), mean_y.T, op_version=opv
             )
 
         y_mean.set_onnx_name_prefix("gpr")
         y_mean_reshaped = OnnxReshapeApi13(
             y_mean,
-            np.array([-1, 1], dtype=np.int64),
+            np.array([-1, mean_y.shape[0]], dtype=np.int64),
             op_version=opv,
             output_names=out[:1],
         )
@@ -192,12 +192,15 @@ def convert_gaussian_process_regressor(
             #     y_var[y_var_negative] = 0.0
             ys0_var = OnnxMax(ys_var, np.array([0], dtype=dtype), op_version=opv)
 
-            if hasattr(op, "_y_train_std") and op._y_train_std != 1:
+            if hasattr(op, "_y_train_std"):
                 # y_var = y_var * self._y_train_std**2
-                ys0_var = OnnxMul(ys0_var, var_y**2, op_version=opv)
+                ys0_var = OnnxMul(
+                    ys0_var, (op._y_train_std**2).astype(dtype), op_version=opv
+                )
 
             # var = np.sqrt(ys0_var)
-            var = OnnxSqrt(ys0_var, output_names=out[1:], op_version=opv)
+            var = OnnxSqrt(ys0_var, op_version=opv)
+            var = OnnxTranspose(var, output_names=out[1:], op_version=opv)
             var.set_onnx_name_prefix("gprv")
             outputs.append(var)
 

tests/test_sklearn_gaussian_process_regressor.py

@@ -24,13 +24,7 @@
     WhiteKernel,
 )
 from sklearn.model_selection import train_test_split
-
-try:
-    # scikit-learn >= 0.22
-    from sklearn.utils._testing import ignore_warnings
-except ImportError:
-    # scikit-learn < 0.22
-    from sklearn.utils.testing import ignore_warnings
+from sklearn.utils._testing import ignore_warnings
 from sklearn.exceptions import ConvergenceWarning
 from skl2onnx.common.data_types import FloatTensorType, DoubleTensorType
 from skl2onnx import to_onnx
@@ -1019,6 +1013,9 @@ def test_gpr_rbf_fitted_return_std_exp_sine_squared_false(self):
             model_onnx,
             verbose=False,
             basename="SklearnGaussianProcessExpSineSquaredStdF-Out0-Dec3",
+            # operator MatMul gets replaced by FusedMatMul but onnxruntime does not check
+            # the availability of the kernel for double.
+            disable_optimisation=True,
         )
         self.check_outputs(
             gp,
@@ -1418,6 +1415,7 @@ def test_x_issue_789_cdist(self):
             pipe.predict(vx1.astype(np.float64)).ravel(), pred[0].ravel()
         )
 
+    @ignore_warnings(category=ConvergenceWarning)
     def test_white_kernel_float(self):
         X, y = make_friedman2(n_samples=500, noise=0, random_state=0)
         tx1, vx1, ty1, vy1 = train_test_split(X, y)
@@ -1436,6 +1434,7 @@ def test_white_kernel_float(self):
             gpr.predict(vx1.astype(np.float32)).ravel(), pred[0].ravel(), rtol=1e-3
         )
 
+    @ignore_warnings(category=ConvergenceWarning)
     def test_white_kernel_double(self):
         X, y = make_friedman2(n_samples=500, noise=0, random_state=0)
         tx1, vx1, ty1, vy1 = train_test_split(X, y)
@@ -1498,11 +1497,29 @@ def test_kernel_white_kernel(self):
         m2 = ker(x, x)
         assert_almost_equal(m2, m1, decimal=5)
 
+    def test_issue_1073_multidimension_process(self):
+        # multioutput gpr
+        n_samples, n_features, n_targets = 1000, 8, 3
+        X, y = make_regression(n_samples, n_features, n_targets=n_targets)
+        tx1, vx1, ty1, vy1 = train_test_split(X, y)
+        model = GaussianProcessRegressor()
+        model.fit(tx1, ty1)
+        initial_type = [("data_in", DoubleTensorType([None, X.shape[1]]))]
+        onx = to_onnx(model, initial_types=initial_type, target_opset=_TARGET_OPSET_)
+        sess = InferenceSession(
+            onx.SerializeToString(), providers=["CPUExecutionProvider"]
+        )
+        pred = sess.run(None, {"data_in": vx1.astype(np.float64)})
+        assert_almost_equal(
+            model.predict(vx1.astype(np.float64)).ravel(), pred[0].ravel()
+        )
+
 
 if __name__ == "__main__":
     # import logging
     # log = logging.getLogger('skl2onnx')
     # log.setLevel(logging.DEBUG)
     # logging.basicConfig(level=logging.DEBUG)
     # TestSklearnGaussianProcessRegressor().test_kernel_white_kernel()
+    # TestSklearnGaussianProcessRegressor().test_issue_1073()
     unittest.main(verbosity=2)