# Now select the features
    X = stations_gdf[feature_columns]
    y = stations_gdf['Temperature']

    # Replace infinite values with NaN
    X = X.replace([np.inf, -np.inf], np.nan)

    # Remove rows with NaN values in X
    valid_indices = ~X.isnull().any(axis=1)
    X = X[valid_indices]
    y = y[valid_indices]
    sample_weights = sample_weights[valid_indices]

    # Define the number of spatial groups (clusters)
    n_clusters = 5  # Adjust based on dataset size and spatial distribution

    # Perform KMeans clustering on the spatial coordinates to define groups for GroupKFold
    kmeans = KMeans(n_clusters=n_clusters, random_state=42)
    stations_gdf_filtered = stations_gdf.loc[valid_indices].copy()
    stations_gdf_filtered['Cluster'] = kmeans.fit_predict(stations_gdf_filtered[['X', 'Y']])

    # Define the groups based on the cluster labels
    groups = stations_gdf_filtered['Cluster']

    # Define GroupKFold cross-validator
    group_kfold = GroupKFold(n_splits=5)

    # Define a custom temporary directory
    custom_tmp_dir = './catboost_tmp'

    # Create the directory if it doesn't exist
    os.makedirs(custom_tmp_dir, exist_ok=True)

    # Define base estimators with updated and more advanced models
    base_estimators = [
        ('rf', RandomForestRegressor(
            random_state=42,
            n_jobs=4,  # Limit CPU usage
            max_depth=25,
            n_estimators=200,
            min_samples_leaf=2,
            min_samples_split=4
        )),
        ('hist_gb', HistGradientBoostingRegressor(
            random_state=42,
            max_iter=200,
            learning_rate=0.05,
            max_depth=8,
            min_samples_leaf=20,
            max_bins=255
        )),
        ('xgb', XGBRegressor(
            objective='reg:squarederror',
            random_state=42,
            tree_method='hist',  # GPU-friendly
            device='cuda',
            sampling_method='gradient_based',
            grow_policy='lossguide',
            max_depth=8,
            subsample=0.8,
            colsample_bytree=0.8,
            learning_rate=0.05,
            n_estimators=200,
            min_child_weight=4,
            gamma=0.2
        )),
        ('catboost', CatBoostRegressor(
            random_seed=42,
            verbose=0,
            depth=8,
            learning_rate=0.05,
            n_estimators=200,
            thread_count=4,
            subsample=0.8,
            rsm=0.8,
            train_dir=custom_tmp_dir
        ))
    ]

    # Updated parameter space for BayesSearchCV focusing on relevant parameters
    param_spaces = {
        # RandomForestRegressor parameters
        'stack__rf__n_estimators': Integer(150, 250),
        'stack__rf__max_depth': Integer(15, 25),
        'stack__rf__min_samples_split': Integer(4, 8),
        'stack__rf__min_samples_leaf': Integer(2, 4),

        # HistGradientBoostingRegressor parameters
        'stack__hist_gb__max_iter': Integer(150, 250),
        'stack__hist_gb__learning_rate': Real(0.03, 0.08, prior='log-uniform'),
        'stack__hist_gb__max_depth': Integer(6, 10),
        'stack__hist_gb__min_samples_leaf': Integer(10, 30),
        'stack__hist_gb__max_bins': Integer(2, 255),

        # XGBRegressor parameters
        'stack__xgb__n_estimators': Integer(150, 250),
        'stack__xgb__learning_rate': Real(0.03, 0.08, prior='log-uniform'),
        'stack__xgb__max_depth': Integer(6, 10),
        'stack__xgb__subsample': Real(0.7, 0.9, prior='uniform'),
        'stack__xgb__colsample_bytree': Real(0.7, 0.9, prior='uniform'),
        'stack__xgb__min_child_weight': Integer(3, 5),
        'stack__xgb__gamma': Real(0.1, 0.3, prior='uniform'),

        # CatBoostRegressor parameters
        'stack__catboost__depth': Integer(6, 10),
        'stack__catboost__learning_rate': Real(0.03, 0.08, prior='log-uniform'),
        'stack__catboost__n_estimators': Integer(150, 250),
        'stack__catboost__subsample': Real(0.7, 0.9, prior='uniform'),
        'stack__catboost__rsm': Real(0.7, 0.9, prior='uniform')
    }

    # Define meta-model
    meta_model = RidgeCV(
        alphas=np.logspace(-3, 3, 7), 
        scoring='r2', 
        store_cv_results=True
    )

    # Define the Stacking Regressor with advanced models
    stacking_regressor = StackingRegressor(
        estimators=base_estimators,
        final_estimator=meta_model,
        passthrough=False,
        n_jobs=4,  # Limited for memory efficiency
        cv=3  # Reduced folds for memory efficiency in final estimator
    )

    # Create a pipeline with StandardScaler and Stacking Regressor
    pipeline = Pipeline([
        ('scaler', StandardScaler()),   # For linear meta-model (Ridge)
        ('stack', stacking_regressor)
    ])

    # Define BayesSearchCV for advanced hyperparameter optimization
    bayes_search = BayesSearchCV(
        estimator=pipeline,
        search_spaces=param_spaces,
        n_iter=40,  # Number of iterations for the search
        cv=group_kfold.split(X, y, groups=groups),
        scoring='r2',
        n_jobs=4,  # Limit for stability and memory usage
        verbose=1,
        random_state=42,
        pre_dispatch='2*n_jobs',  # Memory-efficient dispatch
        error_score='raise'
    )

    # Fit BayesSearchCV with sample weights
    fit_params = {'stack__sample_weight': sample_weights} if sample_weights is not None else {}
    bayes_search.fit(X, y, **fit_params)
    print("BayesSearchCV completed with Spatial GroupKFold and Feature Scaling.")

    # Retrieve the best model
    best_pipeline = bayes_search.best_estimator_
    print(f"Best parameters found: {bayes_search.best_params_}")

Returns this:[14:51:16] WARNING: C:\buildkite-agent\builds\buildkite-windows-cpu-autoscaling-group-i-0015a694724fa8361-1\xgboost\xgboost-ci-windows\src\common\error_msg.cc:58: Falling back to prediction using DMatrix due to mismatched devices. This might lead to higher memory usage and slower performance. XGBoost is running on: cuda:0, while the input data is on: cpu.Potential solutions:- Use a data structure that matches the device ordinal in the booster.- Set the device for booster before call to inplace_predict.Fix that please. I own a GTX 1650 GPU. Note that's I'm using hist in combination with device cuda because that's the implementation starting v2.0.0 instead of gpu_hist.

Here's my XGBoost version documentation on GPU support:
XGBoost GPU Support
This page contains information about GPU algorithms supported in XGBoost.

Note

CUDA 11.0, Compute Capability 5.0 required (See this list to look up compute capability of your GPU card.)

CUDA Accelerated Tree Construction Algorithms
Most of the algorithms in XGBoost including training, prediction and evaluation can be accelerated with CUDA-capable GPUs.

Usage
To enable GPU acceleration, specify the device parameter as cuda. In addition, the device ordinal (which GPU to use if you have multiple devices in the same node) can be specified using the cuda:<ordinal> syntax, where <ordinal> is an integer that represents the device ordinal. XGBoost defaults to 0 (the first device reported by CUDA runtime).

The GPU algorithms currently work with CLI, Python, R, and JVM packages. See Installation Guide for details.

Python example
params = dict()
params["device"] = "cuda"
params["tree_method"] = "hist"
Xy = xgboost.QuantileDMatrix(X, y)
xgboost.train(params, Xy)
With the Scikit-Learn interface
XGBRegressor(tree_method="hist", device="cuda")
GPU-Accelerated SHAP values
XGBoost makes use of GPUTreeShap as a backend for computing shap values when the GPU is used.

booster.set_param({"device": "cuda:0"})
shap_values = booster.predict(dtrain, pred_contribs=True)
shap_interaction_values = model.predict(dtrain, pred_interactions=True)
See Use GPU to speedup SHAP value computation for a worked example.

Multi-node Multi-GPU Training
XGBoost supports fully distributed GPU training using Dask, Spark and PySpark. For getting started with Dask see our tutorial Distributed XGBoost with Dask and worked examples XGBoost Dask Feature Walkthrough, also Python documentation Dask API for complete reference. For usage with Spark using Scala see XGBoost4J-Spark-GPU Tutorial (version 1.6.1+). Lastly for distributed GPU training with PySpark, see Distributed XGBoost with PySpark.

Memory usage
The following are some guidelines on the device memory usage of the hist tree method on GPU.

Memory inside xgboost training is generally allocated for two reasons - storing the dataset and working memory.

The dataset itself is stored on device in a compressed ELLPACK format. The ELLPACK format is a type of sparse matrix that stores elements with a constant row stride. This format is convenient for parallel computation when compared to CSR because the row index of each element is known directly from its address in memory. The disadvantage of the ELLPACK format is that it becomes less memory efficient if the maximum row length is significantly more than the average row length. Elements are quantised and stored as integers. These integers are compressed to a minimum bit length. Depending on the number of features, we usually don’t need the full range of a 32 bit integer to store elements and so compress this down. The compressed, quantised ELLPACK format will commonly use 1/4 the space of a CSR matrix stored in floating point.

Working memory is allocated inside the algorithm proportional to the number of rows to keep track of gradients, tree positions and other per row statistics. Memory is allocated for histogram bins proportional to the number of bins, number of features and nodes in the tree. For performance reasons we keep histograms in memory from previous nodes in the tree, when a certain threshold of memory usage is passed we stop doing this to conserve memory at some performance loss.

If you are getting out-of-memory errors on a big dataset, try the or xgboost.QuantileDMatrix or external memory version. Note that when external memory is used for GPU hist, it’s best to employ gradient based sampling as well. Last but not least, inplace_predict can be preferred over predict when data is already on GPU. Both QuantileDMatrix and inplace_predict are automatically enabled if you are using the scikit-learn interface.

CPU-GPU Interoperability
The model can be used on any device regardless of the one used to train it. For instance, a model trained using GPU can still work on a CPU-only machine and vice versa. For more information about model serialization, see Introduction to Model IO.

Developer notes
The application may be profiled with annotations by specifying USE_NTVX to cmake. Regions covered by the ‘Monitor’ class in CUDA code will automatically appear in the nsight profiler when verbosity is set to 3.

References
Mitchell R, Frank E. (2017) Accelerating the XGBoost algorithm using GPU computing. PeerJ Computer Science 3:e127 https://doi.org/10.7717/peerj-cs.127

NVIDIA Parallel Forall: Gradient Boosting, Decision Trees and XGBoost with CUDA

Out-of-Core GPU Gradient Boosting

Contributors
Many thanks to the following contributors (alphabetical order):

Andrey Adinets

Jiaming Yuan

Jonathan C. McKinney

Matthew Jones

Philip Cho

Rong Ou

Rory Mitchell

Shankara Rao Thejaswi Nanditale

Sriram Chandramouli

Vinay Deshpande

Please report bugs to the XGBoost issues list: https://github.com/dmlc/xgboost/issues. For general questions please visit our user form: https://discuss.xgboost.ai/.