diff --git a/.github/workflows/array-api.yml b/.github/workflows/array-api.yml
new file mode 100644
index 0000000..700278d
--- /dev/null
+++ b/.github/workflows/array-api.yml
@@ -0,0 +1,34 @@
+name: Array API CI
+
+on:
+  push:
+    branches:
+      - maintenance/**
+  pull_request:
+    branches:
+      - main
+      - maintenance/**
+
+jobs:
+  test:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ubuntu-latest, macos-latest]
+        python-version: ["3.14"]
+
+    steps:
+    - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7.0.0
+      with:
+        python-version: ${{ matrix.python-version }}
+    - name: Install dependencies
+      run: |
+        python -m pip install -v .
+        python -m pip install -v --group test
+        python -m pip install torch array-api-strict array-api-compat
+    - name: test
+      run: |
+        cd $RUNNER_TEMP
+        SCIPY_ARRAY_API=1 python -m pytest --pyargs gfdl --cov=gfdl --cov-report=term-missing
diff --git a/src/gfdl/activations.py b/src/gfdl/activations.py
index e4aae04..e2b7b35 100644
--- a/src/gfdl/activations.py
+++ b/src/gfdl/activations.py
@@ -2,8 +2,7 @@
 Activation functions for Gradient Free Deep Learning estimators.
 """
 
-import numpy as np
-import scipy
+from sklearn.utils._array_api import get_namespace_and_device
 
 
 def relu(z):
@@ -17,10 +16,12 @@ def relu(z):
 
     Returns
     -------
-    numpy.ndarray
+    array
         The output array with only positive values.
     """
-    return np.maximum(0, z)
+    xp, _, device = get_namespace_and_device(z)
+    zero = xp.asarray(0, dtype=z.dtype, device=device)
+    return xp.maximum(zero, z)
 
 
 def tanh(z):
@@ -34,14 +35,15 @@ def tanh(z):
 
     Returns
     -------
-    numpy.ndarray
+    array
         The output array with hyperbolic tangent values.
 
     See Also
     --------
     numpy.tanh : The hyperbolic tangent function.
     """
-    return np.tanh(z)
+    xp, _, _device = get_namespace_and_device(z)
+    return xp.tanh(z)
 
 
 def sigmoid(z):
@@ -55,14 +57,15 @@ def sigmoid(z):
 
     Returns
     -------
-    numpy.ndarray
+    array
         The output array with the function values.
 
     See Also
     --------
     scipy.special.expit : The logistic sigmoid function.
     """
-    return scipy.special.expit(z)
+    xp, _, _ = get_namespace_and_device(z)
+    return xp.exp(-xp.logaddexp(xp.zeros_like(z), -z))
 
 
 def identity(z):
@@ -76,7 +79,7 @@ def identity(z):
 
     Returns
     -------
-    numpy.ndarray
+    array
         The input array is returned unchanged.
     """
     return z
@@ -93,14 +96,17 @@ def softmax(z):
 
     Returns
     -------
-    numpy.ndarray
+    array
         The output array with the function values.
 
     See Also
     --------
     scipy.special.softmax : The softmax function.
     """
-    return scipy.special.softmax(z, axis=-1)
+    xp, _, _ = get_namespace_and_device(z)
+    tmp = z - xp.max(z, axis=1, keepdims=True)
+    exp_tmp = xp.exp(tmp)
+    return exp_tmp / xp.sum(exp_tmp, axis=1, keepdims=True)
 
 
 def softmin(z):
@@ -116,7 +122,7 @@ def softmin(z):
 
     Returns
     -------
-    numpy.ndarray
+    array
         The output array with the function values.
 
     See Also
@@ -137,14 +143,15 @@ def log_sigmoid(z):
 
     Returns
     -------
-    numpy.ndarray
+    array
         The output array with the function values.
 
     See Also
     --------
     scipy.special.log_expit : The logistic sigmoid function.
     """
-    return scipy.special.log_expit(z)
+    xp, _, _ = get_namespace_and_device(z)
+    return -xp.logaddexp(xp.zeros_like(z), -z)
 
 
 def log_softmax(z):
@@ -158,14 +165,18 @@ def log_softmax(z):
 
     Returns
     -------
-    numpy.ndarray
+    array
         The output array with the function values.
 
     See Also
     --------
     scipy.special.log_softmax : The logarithm of softmax function.
     """
-    return scipy.special.log_softmax(z, axis=-1)
+    xp, _, _ = get_namespace_and_device(z)
+    shifted = z - xp.max(z, axis=1, keepdims=True)
+    return shifted - xp.log(
+        xp.sum(xp.exp(shifted), axis=1, keepdims=True)
+    )
 
 
 ACTIVATIONS = {
diff --git a/src/gfdl/model.py b/src/gfdl/model.py
index 2f4f216..b85fb4d 100644
--- a/src/gfdl/model.py
+++ b/src/gfdl/model.py
@@ -13,13 +13,19 @@
     RegressorMixin,
 )
 from sklearn.linear_model import Ridge
-from sklearn.preprocessing import OneHotEncoder
+from sklearn.preprocessing import LabelBinarizer, OneHotEncoder
 from sklearn.utils import column_or_1d
+from sklearn.utils._array_api import (
+    check_same_namespace,
+    get_namespace,
+    get_namespace_and_device,
+    move_to,
+)
 from sklearn.utils.metaestimators import available_if
 from sklearn.utils.multiclass import check_classification_targets, unique_labels
 from sklearn.utils.validation import check_is_fitted, validate_data
 
-from gfdl.activations import resolve_activation
+from gfdl.activations import resolve_activation, softmax
 from gfdl.weights import resolve_weight
 
 
@@ -43,7 +49,19 @@ def __init__(
         self.reg_alpha = reg_alpha
         self.rtol = rtol
 
+    def __sklearn_tags__(self):
+        # Fetch default tags or construct them
+        tags = super().__sklearn_tags__()
+        tags.array_api_support = True
+        return tags
+
     def fit(self, X, Y):
+        xp, _, device = get_namespace_and_device(X)
+        Y = xp.astype(
+            move_to(Y, xp=xp, device=device),
+            X.dtype, copy=False
+        )
+
         # Assumption : X, Y have been pre-processed.
         # X shape: (n_samples, n_features)
         # Y shape: (n_samples, n_classes-1)
@@ -65,22 +83,33 @@ def fit(self, X, Y):
         rng = self.get_generator(self.seed)
 
         self.W_.append(
+            xp.asarray(
             self._weight_mode(
                 self._N, hidden_layer_sizes[0], rng=self.get_generator(self.seed)
-                )
+                ),
+            dtype=X.dtype, device=device
             )
+        )
         self.b_.append(
+            xp.asarray(
             self._weight_mode(1, hidden_layer_sizes[0], rng=rng)
-            .reshape(-1)
+            .reshape(-1), dtype=X.dtype, device=device
+            )
             )
         for i, layer in enumerate(hidden_layer_sizes[1:]):
             # (n_hidden, n_features)
             self.W_.append(
-                self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,)
+                xp.asarray(
+                self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,),
+                dtype=X.dtype, device=device
+                )
                 )
             # (n_hidden,)
             self.b_.append(
-                self._weight_mode(1, layer, rng=rng,).reshape(-1)
+                xp.asarray(
+                self._weight_mode(1, layer, rng=rng,).reshape(-1),
+                dtype=X.dtype, device=device
+                )
                 )
 
         # hypothesis space shape: (n_layers,)
@@ -95,7 +124,7 @@ def fit(self, X, Y):
         # or (n_samples, sum_hidden)
         if self.direct_links:
             Hs.append(X)
-        D = np.hstack(Hs)
+        D = xp.concat(Hs, axis=1)
 
         # beta shape: (sum_hidden+n_features, n_classes-1)
         # or (sum_hidden, n_classes-1)
@@ -103,11 +132,15 @@ def fit(self, X, Y):
         # If reg_alpha is None, use direct solve using
         # MoorePenrose Pseudo-Inverse, otherwise use ridge regularized form.
         if self.reg_alpha is None:
-            self.coeff_ = np.linalg.pinv(D, rtol=self.rtol) @ Y
+            self.coeff_ = xp.linalg.pinv(D, rtol=self.rtol) @ Y
         else:
-            ridge = Ridge(alpha=self.reg_alpha, fit_intercept=False)
+            ridge = Ridge(alpha=self.reg_alpha,
+                          fit_intercept=False,
+                          solver="svd",
+                          )
             ridge.fit(D, Y)
-            self.coeff_ = ridge.coef_.T
+            coef = ridge.coef_
+            self.coeff_ = coef.T if coef.ndim == 2 else coef
         return self
 
     def partial_fit(self, X, y):
@@ -132,6 +165,12 @@ def partial_fit(self, X, y):
         #
         # It's evident that it is indeed possible to update normal equations
         # through summation
+        xp, _, device = get_namespace_and_device(X)
+        y = xp.astype(
+            move_to(y, xp=xp, device=device),
+            X.dtype,
+            copy=False,
+        )
 
         if self.reg_alpha is not None and self.reg_alpha < 0.0:
             raise ValueError("Negative reg_alpha. Expected range : None or [0.0, inf).")
@@ -154,22 +193,28 @@ def partial_fit(self, X, y):
             rng = self.get_generator(self.seed)
 
             self.W_.append(
+                xp.asarray(
                 self._weight_mode(
                     self._N, hidden_layer_sizes[0], rng=self.get_generator(self.seed)
-                    )
+                    ), dtype=X.dtype, device=device)
                 )
             self.b_.append(
+                xp.asarray(
                 self._weight_mode(1, hidden_layer_sizes[0], rng=rng)
-                .reshape(-1)
+                .reshape(-1), dtype=X.dtype, device=device)
                 )
             for i, layer in enumerate(hidden_layer_sizes[1:]):
                 # (n_hidden, n_features)
                 self.W_.append(
-                    self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,)
+                    xp.asarray(
+                    self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,),
+                    dtype=X.dtype, device=device)
                     )
                 # (n_hidden,)
                 self.b_.append(
-                    self._weight_mode(1, layer, rng=rng,).reshape(-1)
+                    xp.asarray(
+                    self._weight_mode(1, layer, rng=rng,).reshape(-1),
+                    dtype=X.dtype, device=device)
                     )
 
         # design matrix shape: (n_samples, sum_hidden+n_features)
@@ -185,11 +230,13 @@ def partial_fit(self, X, y):
         # or (n_samples, sum_hidden)
         if self.direct_links:
             Hs.append(X)
-        D = np.hstack(Hs)
+        D = xp.concat(Hs, axis=1)
 
         if not hasattr(self, "A"):
-            self.A = np.zeros((D.shape[1], D.shape[1]))
-            self.B = np.zeros((D.shape[1], y.shape[1]))
+            self.A = xp.zeros((D.shape[1], D.shape[1]),
+                              dtype=X.dtype, device=device)
+            self.B = xp.zeros((D.shape[1], y.shape[1]),
+                              dtype=X.dtype, device=device)
 
         self.A += D.T @ D
         self.B += D.T @ y
@@ -213,22 +260,33 @@ def partial_fit(self, X, y):
         # = (D.T @ D)^-1 @ D.T @ y
 
         if self.reg_alpha is None:
-            self.coeff_ = np.linalg.pinv(self.A, rtol=self.rtol) @ self.B
+            coef = xp.linalg.pinv(self.A, rtol=self.rtol) @ self.B
+            if coef.ndim == 2 and coef.shape[1] == 1:
+                coef = xp.squeeze(coef, axis=-1)
+            self.coeff_ = coef
         else:
-            # scipy.linalg.solve(self.A + reg_mat, self.B)
+            # xp.linalg.solve(self.A + reg_mat, self.B)
             # is equivalent to
             # ridge = Ridge(
             #   alpha=self.reg_alpha, fit_intercept=False, solver='cholesky'
             #   )
             # ridge.fit(D_all[:current_batch], Y_all[:current_batch])
-            reg_mat = np.eye(self.A.shape[0]) * self.reg_alpha
-            self.coeff_ = scipy.linalg.solve(self.A + reg_mat, self.B)
-        if self.coeff_.ndim == 2 and self.coeff_.shape[1] == 1:
-            self.coeff_ = self.coeff_.ravel()
+            reg_mat = (
+                xp.eye(self.A.shape[0], dtype=self.A.dtype, device=device)
+                * self.reg_alpha
+            )
+
+            coef = xp.linalg.solve(self.A + reg_mat, self.B)
+            if coef.ndim == 2 and coef.shape[1] == 1:
+                coef = xp.squeeze(coef, axis=-1)
+            self.coeff_ = coef
         return self
 
     def predict(self, X):
         check_is_fitted(self)
+        xp, _, _ = get_namespace_and_device(X)
+        check_same_namespace(X, self, attribute="coeff_", method="predict")
+
         Hs = []
         H_prev = X
         for W, b in zip(self.W_, self.b_, strict=False):
@@ -238,7 +296,7 @@ def predict(self, X):
 
         if self.direct_links:
             Hs.append(X)
-        D = np.hstack(Hs)
+        D = xp.concat(Hs, axis=1)
         out = D @ self.coeff_
 
         return out
@@ -411,14 +469,24 @@ def fit(self, X, y):
           Fitted estimator.
         """
         # shape: (n_samples, n_features)
-        X, Y = validate_data(self, X, y)
-        self.classes_ = unique_labels(Y)
+        xp, _, _ = get_namespace_and_device(X)
+        X, Y = validate_data(self, X, y,
+                             dtype=[xp.float64, xp.float32])
 
-        # onehot y
-        # (this is necessary for everything beyond binary classification)
-        self.enc_ = OneHotEncoder(handle_unknown="ignore", sparse_output=False)
-        # shape: (n_samples, n_classes-1)
-        Y = self.enc_.fit_transform(Y.reshape(-1, 1))
+        Y = column_or_1d(Y, warn=True)
+        check_classification_targets(Y)
+
+        self.enc_ = LabelBinarizer(sparse_output=False)
+        Y = self.enc_.fit_transform(Y)
+        self.classes_ = self.enc_.classes_
+
+        # Preserve two columns for binary classification.
+        if self.classes_.shape[0] == 2:
+            xp_y, _ = get_namespace(Y)
+            Y = xp_y.concat((1 - Y, Y), axis=1)
+
+        if self.classes_.shape[0] < 2:
+            raise ValueError("y must contain more than 1 class.")
 
         # call base fit method
         super().fit(X, Y)
@@ -466,7 +534,7 @@ def partial_fit(self, X, y, classes=None):
         ``Ridge()`` as we do in full fit. In ``partial_fit()`` we never have access to
         the full design matrix, we only persist the aggregate quantities
         ``D.T @ D``, ``D.T @ y``. Instead, for the regularized case, we directly solve
-        the system using ``scipy.linalg.solve(A + reg_mat, B)``, which is mathematically
+        the system using ``xp.linalg.solve(A + reg_mat, B)``, which is mathematically
         equivalent to scikit-learn's ``Ridge(solver='cholesky')``.
 
         One other difference between full fit and partial fit arises in the direct
@@ -475,24 +543,39 @@ def partial_fit(self, X, y, classes=None):
         numerical stability and avoid loss of information.
         """
         # shape: (n_samples, n_features)
-        X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_"))
+        xp, _, _ = get_namespace_and_device(X)
+        X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_"),
+                             dtype=[xp.float64, xp.float32])
         Y = column_or_1d(Y, warn=True)
         check_classification_targets(Y)
 
         if not hasattr(self, "classes_"):
             if classes is None:
                 raise TypeError("Classes must not be None for first partial_fit call")
-            self.classes_ = classes
-            self.enc_ = OneHotEncoder(handle_unknown="ignore", sparse_output=False)
-            self.enc_.fit(self.classes_.reshape(-1, 1))
-
-        if not set(unique_labels(Y)) <= set(self.classes_):
+            xp_y, _, device_y_ = get_namespace_and_device(Y)
+            classes = move_to(classes, xp=xp_y, device=device_y_)
+            self.enc_ = LabelBinarizer(sparse_output=False).fit(classes)
+            self.classes_ = self.enc_.classes_
+            if self.classes_.shape[0] < 2:
+                raise ValueError("classes must contain more than 1 class")
+
+        # Check membership without Python sets or converting labels to NumPy.
+        xp_y, _, device_y = get_namespace_and_device(self.classes_)
+        Y = move_to(Y, xp=xp_y, device=device_y)
+        matches = xp_y.expand_dims(Y, axis=1) == xp_y.expand_dims(
+            self.classes_, axis=0
+        )
+        if not xp_y.all(xp_y.any(matches, axis=1)):
             raise ValueError(
-                f"Expected only labels in classes_ = {list(self.classes_)!r}, "
-                f"but got {unique_labels(Y)!r}."
+                f"Expected only labels in classes_ = {self.classes_!r},"
+                f"but got {unique_labels(y)!r}."
             )
         Y = self.enc_.transform(Y.reshape(-1, 1))
 
+        # Preserve two output columns for binary classification.
+        if self.classes_.shape[0] == 2:
+            Y = xp_y.concat((1 - Y, Y), axis=1)
+
         super().partial_fit(X, Y)
         return self
 
@@ -511,9 +594,13 @@ def predict(self, X):
           The predicted classes, with shape (n_samples,) or (n_samples, n_outputs).
         """
         check_is_fitted(self)
-        X = validate_data(self, X, reset=False)
+        check_same_namespace(X, self, attribute="coeff_", method="predict")
+        xp, _, device = get_namespace_and_device(X)
+        X = validate_data(self, X, reset=False,
+                          dtype=[xp.float64, xp.float32])
         out = self.predict_proba(X)
-        y_hat = self.classes_[np.argmax(out, axis=1)]
+        classes = move_to(self.classes_, xp=xp, device=device)
+        y_hat = classes[xp.argmax(out, axis=1)]
         return y_hat
 
     def predict_proba(self, X):
@@ -533,10 +620,12 @@ def predict_proba(self, X):
           have shape (n_samples, n_classes).
         """
         check_is_fitted(self)
-        X = validate_data(self, X, reset=False)
+        check_same_namespace(X, self, attribute="coeff_", method="predict_proba")
+        xp, _ = get_namespace(X)
+        X = validate_data(self, X, reset=False,
+                          dtype=[xp.float64, xp.float32])
         out = super().predict(X)
-        out = np.exp(out - logsumexp(out, axis=1, keepdims=True))
-        return out
+        return softmax(out)
 
 
 class EnsembleGFDL(BaseEstimator):
@@ -876,7 +965,9 @@ def fit(self, X, y):
           The fitted estimator.
         """
         # shape: (n_samples, n_features)
-        X, Y = validate_data(self, X, y)
+        xp, _, _ = get_namespace_and_device(X)
+        X, Y = validate_data(self, X, y,
+                             dtype=[xp.float64, xp.float32])
         self.classes_ = unique_labels(Y)
 
         # onehot y
@@ -940,7 +1031,9 @@ def partial_fit(self, X, y, classes=None):
         numerical stability and avoid loss of information.
         """
         # shape: (n_samples, n_features)
-        X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_"))
+        xp, _, _ = get_namespace_and_device(X)
+        X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_"),
+                             dtype=[xp.float64, xp.float32])
         Y = column_or_1d(Y, warn=True)
         check_classification_targets(Y)
 
@@ -988,7 +1081,9 @@ def predict_proba(self, X):
           have shape (n_samples, n_classes).
         """
         check_is_fitted(self)
-        X = validate_data(self, X, reset=False)
+        xp, _, _ = get_namespace_and_device(X)
+        X = validate_data(self, X, reset=False,
+                          dtype=[xp.float64, xp.float32])
 
         outs = self._forward(X)
         probs = []
@@ -1014,7 +1109,9 @@ def predict(self, X):
           The predicted classes, with shape (n_samples,) or (n_samples, n_outputs).
         """
         check_is_fitted(self)
-        X = validate_data(self, X, reset=False)
+        xp, _, _ = get_namespace_and_device(X)
+        X = validate_data(self, X, reset=False,
+                          dtype=[xp.float64, xp.float32])
 
         if self.voting == "soft":
             P = self.predict_proba(X)
@@ -1193,7 +1290,9 @@ def fit(self, X, y):
         object
           The fitted estimator.
         """
-        X, Y = validate_data(self, X, y, multi_output=True)
+        xp, _, _ = get_namespace_and_device(X)
+        X, Y = validate_data(self, X, y, multi_output=True,
+                             dtype=[xp.float64, xp.float32])
         super().fit(X, Y)
         return self
 
@@ -1241,7 +1340,9 @@ def partial_fit(self, X, y):
         the condition number is squared. This may require a lower ``rtol`` to maintain
         numerical stability and avoid loss of information.
         """
-        X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_"))
+        xp, _, _ = get_namespace_and_device(X)
+        X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_"),
+                             dtype=[xp.float64, xp.float32])
         if Y.ndim == 1:
             Y = Y.reshape(-1, 1)
         super().partial_fit(X, Y)
@@ -1264,5 +1365,8 @@ def predict(self, X):
           (n_samples, n_outputs).
         """
         check_is_fitted(self)
-        X = validate_data(self, X, reset=False)
+        check_same_namespace(X, self, attribute="coeff_", method="predict")
+        xp, _, _ = get_namespace_and_device(X)
+        X = validate_data(self, X, reset=False,
+                          dtype=[xp.float64, xp.float32])
         return super().predict(X)
diff --git a/src/gfdl/tests/test_model.py b/src/gfdl/tests/test_model.py
index 66ce1af..50e4dd4 100644
--- a/src/gfdl/tests/test_model.py
+++ b/src/gfdl/tests/test_model.py
@@ -1,12 +1,26 @@
+
 import numpy as np
 import pytest
 from numpy.testing import assert_allclose
+from sklearn import config_context
 from sklearn.base import clone
 from sklearn.datasets import load_breast_cancer, load_digits, make_classification
+from sklearn.externals import array_api_extra as xpx
 from sklearn.metrics import accuracy_score, roc_auc_score
 from sklearn.model_selection import StratifiedKFold, train_test_split
 from sklearn.preprocessing import OneHotEncoder, StandardScaler
-from sklearn.utils.estimator_checks import parametrize_with_checks
+from sklearn.utils._array_api import (
+    _atol_for_type,
+    move_to,
+    yield_namespace_device_dtype_combinations,
+)
+from sklearn.utils._testing import (
+    _array_api_for_tests,
+    skip_if_array_api_compat_not_configured,
+)
+from sklearn.utils.estimator_checks import (
+    parametrize_with_checks,
+)
 from ucimlrepo import fetch_ucirepo
 
 from gfdl.activations import ACTIVATIONS
@@ -329,11 +343,11 @@ def test_soft_and_hard_can_differ(alpha):
 @pytest.mark.parametrize("Classifier", [GFDLClassifier, EnsembleGFDLClassifier])
 def test_invalid_activation_weight(Classifier):
     X = np.zeros((30, 4))
-    y = np.zeros((30,))
-    invalid_act = Classifier(hidden_layer_sizes=100,
+    y = np.arange(30) % 2
+    invalid_act = Classifier(hidden_layer_sizes=(100,),
                              activation="bogus_activation",
                              weight_scheme="uniform")
-    invalid_weight = Classifier(hidden_layer_sizes=100,
+    invalid_weight = Classifier(hidden_layer_sizes=(100,),
                                 activation="identity",
                                 weight_scheme="bogus_weight")
     # the sklearn estimator API bans input validation in __init__,
@@ -351,7 +365,7 @@ def test_invalid_alpha(Classifier):
     # so we need to call fit() for error handling to kick in:
     # https://scikit-learn.org/stable/developers/develop.html#developing-scikit-learn-estimators
     X = np.zeros((30, 4))
-    y = np.zeros((30,))
+    y = np.arange(30) % 2
     bad_est = Classifier(hidden_layer_sizes=100,
                              activation="identity",
                              weight_scheme="uniform",
@@ -496,7 +510,9 @@ def test_rtol_classifier(reg_alpha, rtol, expected_acc, expected_roc):
     roc_cur = roc_auc_score(y_test, y_hat_cur_proba, multi_class="ovo")
 
     np.testing.assert_allclose(acc_cur, expected_acc)
-    np.testing.assert_allclose(roc_cur, expected_roc)
+    np.testing.assert_allclose(roc_cur, expected_roc,
+                               rtol=1e-4
+                               )
 
 
 @pytest.mark.parametrize("reg_alpha, rtol, expected_acc, expected_roc", [
@@ -822,3 +838,155 @@ def test_gh_85_classifiers(hidden_layer_sizes, classifier):
     clf = classifier(hidden_layer_sizes=hidden_layer_sizes, seed=0)
     with pytest.raises(ValueError, match="must be > 0"):
         clf.fit(X, y)
+
+
+@skip_if_array_api_compat_not_configured
+@pytest.mark.parametrize(
+        "reg_alpha, hidden_layer_sizes, activation, weight_scheme, direct_links",
+        [
+            (2, (10,), "relu", "normal", True),
+            (2, (10, 10), "sigmoid", "zeros", True),
+            (2, (10,), "softmax", "lecun_normal", True),
+            (2, (10,), "tanh", "he_normal", False),
+            (None, (10,), "softmin", "glorot_normal", True),
+            (None, (10, 10), "log_softmax", "uniform", False),
+            (None, (10,), "log_sigmoid", "lecun_uniform", True),
+            (2, (10,), "identity", "glorot_uniform", True),
+            (2, (10,), "identity", "he_uniform", True),
+        ]
+)
+@pytest.mark.parametrize("n_classes", [2, 5])
+@pytest.mark.parametrize(
+    "array_namespace, device_name, dtype_name",
+    yield_namespace_device_dtype_combinations(),
+)
+def test_xp_classifier_fit(hidden_layer_sizes,
+                            direct_links,
+                            activation,
+                            weight_scheme,
+                            reg_alpha,
+                            n_classes,
+                            array_namespace,
+                            device_name,
+                            dtype_name,
+                            ):
+    # Custom make_regression()-based NumPy and array API matching
+    # This test is closely based on `_check_array_api_core` from sklearn.
+    # Here, we have better control of `atol`, `rtol`, and simulating data.
+    # https://github.com/scikit-learn/scikit-learn/pull/33215#issuecomment-5953082117.
+    xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name)
+
+    atol = _atol_for_type(dtype_name) * 10
+    estimator = GFDLClassifier
+
+    X, y = make_classification(n_samples=500, n_features=20,
+        n_informative=5, n_classes=n_classes, random_state=0,
+    )
+    X = X.astype(dtype_name, copy=False)
+    X_xp = xp.asarray(X, device=device)
+
+    # NumPy .fit
+    est = estimator(
+        reg_alpha=reg_alpha,
+        hidden_layer_sizes=hidden_layer_sizes,
+        activation=activation,
+        weight_scheme=weight_scheme,
+        direct_links=direct_links,
+        seed=0,)
+    with config_context(array_api_dispatch=False):
+        est.fit(X, y)
+
+    # xp (array API) .fit
+    est_xp = clone(est)
+    with config_context(array_api_dispatch=True):
+        est_xp.fit(X_xp, y)
+
+    # Check fitted attributes are close
+    array_attributes = {
+        key: value for key, value in vars(est).items() if isinstance(value, np.ndarray)
+    }
+    for attribute_name, attribute_value in array_attributes.items():
+        est_xp_attr = getattr(est_xp, attribute_name)
+        est_xp_attr_np = move_to(est_xp_attr, xp=np, device="cpu")
+        xpx.testing.assert_close(
+            attribute_value, est_xp_attr_np,
+            err_msg=f"{attribute_name} not the same",
+            atol=atol,
+        )
+
+    # Check method outputs are close
+    methods = [
+        "predict",
+        "predict_proba",
+    ]
+    for method_name in methods:
+        with config_context(array_api_dispatch=False):
+            result = getattr(est, method_name)(X)
+        with config_context(array_api_dispatch=True):
+            result_xp = getattr(est_xp, method_name)(X_xp)
+        result_xp_np = move_to(result_xp, xp=np, device="cpu")
+        xpx.testing.assert_close(
+            result, result_xp_np,
+            err_msg=f"{method_name}() did not the return the same result",
+            atol=atol,
+        )
+
+
+@skip_if_array_api_compat_not_configured
+@pytest.mark.parametrize("reg_alpha", [None, 2])
+@pytest.mark.parametrize("hidden_layer_sizes", [(10,), (10, 10)])
+@pytest.mark.parametrize("n_classes", [2, 5])
+@pytest.mark.parametrize(
+    "array_namespace, device_name, dtype_name",
+    yield_namespace_device_dtype_combinations(),
+)
+def test_xp_classifier_partial_fit(hidden_layer_sizes,
+                                    reg_alpha,
+                                    n_classes,
+                                    array_namespace,
+                                    device_name,
+                                    dtype_name,
+                                    ):
+
+    atol = _atol_for_type(dtype_name) * 10
+    batch = 2
+
+    xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name)
+    X_np, y_np = make_classification(
+        n_samples=6, n_features=20, n_informative=5, random_state=0,
+        n_classes=n_classes
+    )
+    X_np = X_np.astype(dtype_name, copy=False)
+    X_xp = xp.asarray(X_np, device=device)
+
+    params = {
+        "reg_alpha": reg_alpha,
+        "hidden_layer_sizes": hidden_layer_sizes,
+        "seed": 0
+    }
+    estimator = GFDLClassifier
+
+    predictions = []
+    for X, dispatch in [(X_np, False), (X_xp, True)]:
+        with config_context(array_api_dispatch=dispatch):
+            model = estimator(**params)
+            for start in range(0, X.shape[0], batch):
+                model.partial_fit(
+                    X[start:start + batch, ...], y_np[start:start + batch],
+                    classes=np.unique(y_np)
+                )
+            predictions.extend(
+                (move_to(model.predict(X), xp=np, device="cpu"),
+                move_to(model.predict_proba(X), xp=np, device="cpu"))
+            )
+
+    xpx.testing.assert_close(
+        predictions[0], predictions[2],
+        err_msg="predict() did not the return the same result",
+        atol=atol,
+    )
+    xpx.testing.assert_close(
+        predictions[1], predictions[3],
+        err_msg="predict() did not the return the same result",
+        atol=atol,
+    )
diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py
index c0cab55..a0b6c89 100644
--- a/src/gfdl/tests/test_regression.py
+++ b/src/gfdl/tests/test_regression.py
@@ -1,11 +1,26 @@
+
+
 import numpy as np
 import pytest
+from sklearn import config_context
 from sklearn.base import clone
 from sklearn.datasets import fetch_openml, make_regression
+from sklearn.externals import array_api_extra as xpx
 from sklearn.metrics import r2_score
 from sklearn.model_selection import train_test_split
 from sklearn.preprocessing import StandardScaler
-from sklearn.utils.estimator_checks import parametrize_with_checks
+from sklearn.utils._array_api import (
+    _atol_for_type,
+    move_to,
+    yield_namespace_device_dtype_combinations,
+)
+from sklearn.utils._testing import (
+    _array_api_for_tests,
+    skip_if_array_api_compat_not_configured,
+)
+from sklearn.utils.estimator_checks import (
+    parametrize_with_checks,
+)
 
 from gfdl.activations import ACTIVATIONS
 from gfdl.model import GFDLRegressor
@@ -257,3 +272,142 @@ def test_preserve_class_inputs():
     for k, v in actual.items():
         assert v == expected[k]
         assert isinstance(v, type(expected[k]))
+
+
+@skip_if_array_api_compat_not_configured
+@pytest.mark.parametrize(
+        "reg_alpha, hidden_layer_sizes, activation, weight_scheme, direct_links",
+        [
+            (2, (10,), "relu", "normal", True),
+            (2, (10, 10), "sigmoid", "zeros", True),
+            (2, (10,), "softmax", "lecun_normal", True),
+            (2, (10,), "tanh", "he_normal", False),
+            (None, (10,), "softmin", "glorot_normal", True),
+            (None, (10, 10), "log_softmax", "uniform", False),
+            (None, (10,), "log_sigmoid", "lecun_uniform", True),
+            (2, (10,), "identity", "glorot_uniform", True),
+            (2, (10,), "identity", "he_uniform", True),
+        ]
+)
+@pytest.mark.parametrize(
+    "array_namespace, device_name, dtype_name",
+    yield_namespace_device_dtype_combinations(),
+)
+def test_xp_regressor_fit(hidden_layer_sizes,
+                            direct_links,
+                            activation,
+                            weight_scheme,
+                            reg_alpha,
+                            array_namespace,
+                            device_name,
+                            dtype_name,
+                            ):
+    # Custom make_regression()-based NumPy and array API matching
+    # This test is closely based on `_check_array_api_core` from sklearn.
+    # Here, we have better control of `atol`, `rtol`, and simulating data.
+    # https://github.com/scikit-learn/scikit-learn/pull/33215#issuecomment-5953082117.
+    xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name)
+
+    atol = _atol_for_type(dtype_name) * 10
+    estimator = GFDLRegressor
+
+    X, y = make_regression(
+        n_samples=500, n_features=20, n_informative=5, n_targets=1, random_state=0,
+    )
+    X = X.astype(dtype_name, copy=False)
+    X_xp = xp.asarray(X, device=device)
+
+    # NumPy .fit
+    est = estimator(
+        reg_alpha=reg_alpha,
+        hidden_layer_sizes=hidden_layer_sizes,
+        activation=activation,
+        weight_scheme=weight_scheme,
+        direct_links=direct_links,
+        seed=0,)
+    with config_context(array_api_dispatch=False):
+        est.fit(X, y)
+
+    # xp (array API) .fit
+    est_xp = clone(est)
+    with config_context(array_api_dispatch=True):
+        est_xp.fit(X_xp, y)
+
+    # Check fitted attributes are close
+    array_attributes = {
+        key: value for key, value in vars(est).items() if isinstance(value, np.ndarray)
+    }
+    for attribute_name, attribute_value in array_attributes.items():
+        est_xp_attr = getattr(est_xp, attribute_name)
+        est_xp_attr_np = move_to(est_xp_attr, xp=np, device="cpu")
+        xpx.testing.assert_close(
+            attribute_value, est_xp_attr_np,
+            err_msg=f"{attribute_name} not the same",
+            atol=atol,
+        )
+
+    # Check method outputs are close
+    methods = [
+        "predict",
+    ]
+    for method_name in methods:
+        with config_context(array_api_dispatch=False):
+            result = getattr(est, method_name)(X)
+        with config_context(array_api_dispatch=True):
+            result_xp = getattr(est_xp, method_name)(X_xp)
+        result_xp_np = move_to(result_xp, xp=np, device="cpu")
+        xpx.testing.assert_close(
+            result, result_xp_np,
+            err_msg=f"{method_name}() did not the return the same result",
+            atol=atol,
+        )
+
+
+@skip_if_array_api_compat_not_configured
+@pytest.mark.parametrize("reg_alpha", [None, 2])
+@pytest.mark.parametrize("hidden_layer_sizes", [(10,), (10, 10)])
+@pytest.mark.parametrize(
+    "array_namespace, device_name, dtype_name",
+    yield_namespace_device_dtype_combinations(),
+)
+def test_xp_regressor_partial_fit(hidden_layer_sizes,
+                                    reg_alpha,
+                                    array_namespace,
+                                    device_name,
+                                    dtype_name,
+                                    ):
+
+    atol = _atol_for_type(dtype_name) * 10
+    batch = 2
+
+    xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name)
+    X_np, y_np = make_regression(
+        n_samples=6, n_features=20, n_informative=5, random_state=0,
+    )
+    X_np = X_np.astype(dtype_name, copy=False)
+    X_xp = xp.asarray(X_np, device=device)
+
+    params = {
+        "reg_alpha": reg_alpha,
+        "hidden_layer_sizes": hidden_layer_sizes,
+        "seed": 0
+    }
+    estimator = GFDLRegressor
+
+    predictions = []
+    for X, dispatch in [(X_np, False), (X_xp, True)]:
+        with config_context(array_api_dispatch=dispatch):
+            model = estimator(**params)
+            for start in range(0, X.shape[0], batch):
+                model.partial_fit(
+                    X[start:start + batch, ...], y_np[start:start + batch],
+                )
+            predictions.append(
+                move_to(model.predict(X), xp=np, device="cpu")
+            )
+
+    xpx.testing.assert_close(
+        predictions[0], predictions[1],
+        err_msg="predict() did not the return the same result",
+        atol=atol,
+    )