In this tutorial, we implement NVIDIA cuML as a GPU-accelerated machine learning framework and build a practical workflow that demonstrates how RAPIDS can accelerate familiar data science and machine learning tasks. We begin by configuring the GPU environment and examining cuml.accel, which lets us accelerate existing scikit-learn workloads with minimal code changes, before moving to the native cuML API for direct CuPy and cuDF interoperability. We then benchmark CPU and GPU implementations of PCA, K-Means, nearest-neighbor search, logistic regression, random forests, and DBSCAN, while using synchronized timing to obtain meaningful performance measurements. We also build GPU-based manifold-learning and clustering pipelines with UMAP, t-SNE, HDBSCAN, and trustworthiness metrics; explore high-throughput forest inference with FIL; validate GPU-generated SHAP explanations; perform hyperparameter optimization with scikit-learn meta-estimators; and finally serialize trained models while examining portability between GPU and CPU environments.
import os
import sys
import time
import json
import shutil
import warnings
import subprocess
import importlib
import traceback
warnings.filterwarnings("ignore")
QUICK = False
SEED = 42
SCALE = 0.25 if QUICK else 1.0
N_MAIN = int(200_000 * SCALE)
D_MAIN = 64
N_RF = int(50_000 * SCALE)
D_RF = 32
N_NN_INDEX = int(50_000 * SCALE)
N_NN_QUERY = int(5_000 * SCALE)
N_DBSCAN = int(20_000 * SCALE)
N_MANIFOLD = int(60_000 * SCALE)
N_ACCEL = int(80_000 * SCALE)
RESULTS = []
NOTES = []
def banner(title):
line = "=" * 78
print(f"\n{line}\n {title}\n{line}", flush=True)
def section(title, fn, *args, **kwargs):
banner(title)
t0 = time.perf_counter()
try:
fn(*args, **kwargs)
except Exception:
print(f"[!] Section skipped due to an error:\n{traceback.format_exc()}")
print(f"[section wall time: {time.perf_counter() - t0:.1f}s]", flush=True)
def bootstrap():
if shutil.which("nvidia-smi") is None:
raise SystemExit(
"No NVIDIA GPU found. In Colab: Runtime > Change runtime type > GPU."
)
print(subprocess.run(
["nvidia-smi",
"--query-gpu=name,memory.total,compute_cap,driver_version",
"--format=csv"],
capture_output=True, text=True).stdout)
try:
import cuml
print("cuML already available — skipping install.")
except ImportError:
print("Installing RAPIDS cuML (this takes ~1-3 minutes)...")
pin = ""
try:
import cudf
major_minor = ".".join(cudf.__version__.split("+")[0].split(".")[:2])
pin = f"=={major_minor}.*"
print(f" Pinning to the preinstalled cuDF line: cuml-cu12{pin}")
except Exception:
print(" cuDF not found; installing the latest stable cuml-cu12.")
cmd = [sys.executable, "-m", "pip", "install", "-q",
"--extra-index-url=https://pypi.nvidia.com", f"cuml-cu12{pin}"]
print("$ " + " ".join(cmd))
rc = subprocess.run(cmd).returncode
if rc != 0:
raise SystemExit(
"pip install failed. Alternative that always works on Colab:\n"
" !git clone https://github.com/rapidsai/rapidsai-csp-utils.git\n"
" !python rapidsai-csp-utils/colab/pip-install.py"
)
importlib.invalidate_caches()
import cuml
import cupy
print(f"cuml {cuml.__version__}")
print(f"cupy {cupy.__version__}")
try:
import cudf
print(f"cudf {cudf.__version__}")
except Exception:
pass
import sklearn
print(f"sklearn {sklearn.__version__} (cuML requires scikit-learn >= 1.6)")
bootstrap()
import numpy as np
import cupy as cp
import cuml
import matplotlib.pyplot as plt
from cuml.datasets import make_classification as gpu_make_classification
from cuml.datasets import make_blobs as gpu_make_blobs
rng = np.random.RandomState(SEED)
cp.random.seed(SEED)
class Timer:
def __init__(self, label, sync=True):
self.label = label
self.sync = sync
def __enter__(self):
if self.sync:
cp.cuda.runtime.deviceSynchronize()
self.t0 = time.perf_counter()
return self
def __exit__(self, *exc):
if self.sync:
cp.cuda.runtime.deviceSynchronize()
self.dt = time.perf_counter() - self.t0
print(f" {self.label:<44s} {self.dt:8.3f}s")
return False
def to_numpy(a):
if isinstance(a, cp.ndarray):
return cp.asnumpy(a)
if hasattr(a, "to_numpy"):
return a.to_numpy()
return np.asarray(a)
def record(task, cpu_s, gpu_s):
RESULTS.append((task, cpu_s, gpu_s))
if cpu_s and gpu_s:
print(f" -> {task}: {cpu_s / gpu_s:.1f}x speedup\n")
ACCEL_SCRIPT = f'''
import time
import numpy as np
from sklearn.datasets import make_blobs
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn.neighbors import NearestNeighbors
from sklearn.linear_model import Ridge
X, y = make_blobs(n_samples={N_ACCEL}, n_features=32, centers=12, random_state=0)
X = X.astype("float32"); y = y.astype("float32")
t0 = time.perf_counter()
PCA(n_components=8).fit_transform(X)
KMeans(n_clusters=12, n_init=1, random_state=0).fit(X)
NearestNeighbors(n_neighbors=8).fit(X[:{N_ACCEL // 2}]).kneighbors(X[:5000])
Ridge(alpha=1.0).fit(X, y)
Ridge(alpha=1.0, positive=True).fit(X[:5000], y[:5000])
print("MODELTIME %.3f" % (time.perf_counter() - t0))
'''
def demo_accel():
path = "/content/_accel_demo.py" if os.path.isdir("/content") else "_accel_demo.py"
with open(path, "w") as f:
f.write(ACCEL_SCRIPT)
def run(cmd, label):
print(f"\n$ {' '.join(cmd[1:])}")
t0 = time.perf_counter()
p = subprocess.run(cmd, capture_output=True, text=True)
wall = time.perf_counter() - t0
out = p.stdout + p.stderr
model_s = None
for line in out.splitlines():
if line.startswith("MODELTIME"):
model_s = float(line.split()[1])
print(out.strip()[:4000])
print(f"[{label}] model time = {model_s}s | process wall = {wall:.1f}s")
return model_s
cpu_s = run([sys.executable, path], "stock sklearn")
cmd = [sys.executable, "-m", "cuml.accel", "--profile", path]
gpu_s = run(cmd, "cuml.accel")
if gpu_s is None:
gpu_s = run([sys.executable, "-m", "cuml.accel", path], "cuml.accel")
record("cuml.accel (sklearn script, unmodified)", cpu_s, gpu_s)
NOTES.append(
"cuml.accel needed ZERO source changes; the profile table above shows "
"which calls ran on GPU and why Ridge(positive=True) fell back to CPU."
)
We configure the tutorial environment, define dataset sizes and benchmarking utilities, and verify that an NVIDIA GPU is available. We install and initialize RAPIDS cuML when necessary, set up CuPy and reproducibility controls, and create synchronized timing and result-tracking helpers. We also demonstrate cuml.accel by running an unmodified scikit-learn workload and comparing its CPU execution with GPU-accelerated execution.
def demo_native_api():
from cuml.preprocessing import StandardScaler
from cuml.model_selection import train_test_split
X, y = gpu_make_blobs(n_samples=50_000, n_features=8, centers=5,
random_state=SEED, dtype=np.float32)
print(f"cuml.datasets output lives on device: {type(X).__module__}, "
f"shape={X.shape}, dtype={X.dtype}")
try:
import cudf
df = cudf.DataFrame(X, columns=[f"f{i}" for i in range(X.shape[1])])
back = df.values
ptr_a = X.__cuda_array_interface__["data"][0]
ptr_b = back.__cuda_array_interface__["data"][0]
print(f"CuPy ptr = {hex(ptr_a)}")
print(f"cuDF->CuPy= {hex(ptr_b)}")
print("Same device pointer (true zero-copy)? ", ptr_a == ptr_b)
print("Note: a column-major DataFrame round trip may re-pack; what "
"matters is that no host (CPU) round trip ever happens.")
scaled = StandardScaler().fit_transform(df)
print(f"StandardScaler(cuDF) -> {type(scaled).__name__}")
except Exception as e:
print(f"cuDF interop skipped: {e}")
from cuml.decomposition import PCA
pca = PCA(n_components=3).fit(X)
print(f"default (mirrors input) -> {type(pca.transform(X)).__name__}")
with cuml.using_output_type("numpy"):
print(f"inside using_output_type() -> {type(pca.transform(X)).__name__}")
print(f"after the context manager -> {type(pca.transform(X)).__name__}")
NOTES.append(
"Keep output_type as CuPy/cuDF inside a pipeline; converting to NumPy "
"on every step forces a device->host copy and eats the speedup."
)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.2, random_state=SEED)
print(f"train_test_split -> {Xtr.shape} / {Xte.shape}, still on device: "
f"{isinstance(Xtr, cp.ndarray)}")
We work directly with the native cuML API and explore how GPU-resident data moves between CuPy, cuDF, and cuML components. We inspect device pointers to understand zero-copy interoperability and use cuML output-type controls to manage whether results remain on the GPU or return as NumPy arrays. We also perform a GPU-native train-test split so that our data remains on the device throughout the workflow.
def demo_benchmarks():
from sklearn.decomposition import PCA as skPCA
from sklearn.cluster import KMeans as skKMeans, DBSCAN as skDBSCAN
from sklearn.neighbors import NearestNeighbors as skNN
from sklearn.linear_model import LogisticRegression as skLR
from sklearn.ensemble import RandomForestClassifier as skRF
from cuml.decomposition import PCA as cuPCA
from cuml.cluster import KMeans as cuKMeans, DBSCAN as cuDBSCAN
from cuml.neighbors import NearestNeighbors as cuNN
from cuml.linear_model import LogisticRegression as cuLR
from cuml.ensemble import RandomForestClassifier as cuRF
print(f"Generating {N_MAIN:,} x {D_MAIN} on the GPU...")
Xg, yg = gpu_make_classification(n_samples=N_MAIN, n_features=D_MAIN,
n_informative=32, n_classes=4,
random_state=SEED)
Xg = Xg.astype(cp.float32)
yg = yg.astype(cp.int32)
Xc, yc = cp.asnumpy(Xg), cp.asnumpy(yg)
print(f" device array: {Xg.nbytes / 1e6:.0f} MB\n")
print("PCA (n_components=16)")
with Timer("sklearn", sync=False) as t:
skPCA(n_components=16, random_state=SEED).fit_transform(Xc)
cpu = t.dt
with Timer("cuML") as t:
cuPCA(n_components=16, random_state=SEED).fit_transform(Xg)
record("PCA", cpu, t.dt)
print("KMeans (k=16)")
with Timer("sklearn", sync=False) as t:
skKMeans(n_clusters=16, n_init=1, max_iter=100,
random_state=SEED).fit(Xc)
cpu = t.dt
with Timer("cuML") as t:
cuKMeans(n_clusters=16, n_init=1, max_iter=100,
random_state=SEED).fit(Xg)
record("KMeans", cpu, t.dt)
print(f"NearestNeighbors k=16 ({N_NN_INDEX:,} index / {N_NN_QUERY:,} query)")
idx_g, q_g = Xg[:N_NN_INDEX], Xg[N_NN_INDEX:N_NN_INDEX + N_NN_QUERY]
idx_c, q_c = cp.asnumpy(idx_g), cp.asnumpy(q_g)
with Timer("sklearn (brute)", sync=False) as t:
skNN(n_neighbors=16, algorithm="brute", n_jobs=-1).fit(idx_c).kneighbors(q_c)
cpu = t.dt
with Timer("cuML") as t:
d_gpu, i_gpu = cuNN(n_neighbors=16).fit(idx_g).kneighbors(q_g)
record("NearestNeighbors", cpu, t.dt)
print("LogisticRegression (multinomial, lbfgs/QN)")
with Timer("sklearn", sync=False) as t:
sk_lr = skLR(max_iter=200, n_jobs=-1).fit(Xc, yc)
cpu = t.dt
with Timer("cuML") as t:
cu_lr = cuLR(max_iter=200).fit(Xg, yg)
record("LogisticRegression", cpu, t.dt)
print(f" accuracy sklearn={sk_lr.score(Xc, yc):.4f} "
f"cuML={float((cu_lr.predict(Xg) == yg).mean()):.4f} "
"(different solvers, so small differences are expected)\n")
print(f"RandomForestClassifier (100 trees, depth 12, {N_RF:,} x {D_RF})")
Xr_g, yr_g = gpu_make_classification(n_samples=N_RF, n_features=D_RF,
n_informative=16, n_classes=2,
random_state=SEED)
Xr_g = Xr_g.astype(cp.float32)
yr_g = yr_g.astype(cp.int32)
Xr_c, yr_c = cp.asnumpy(Xr_g), cp.asnumpy(yr_g)
with Timer("sklearn", sync=False) as t:
skRF(n_estimators=100, max_depth=12, n_jobs=-1,
random_state=SEED).fit(Xr_c, yr_c)
cpu = t.dt
with Timer("cuML") as t:
cu_rf = cuRF(n_estimators=100, max_depth=12, n_bins=128,
n_streams=4, random_state=SEED).fit(Xr_g, yr_g)
record("RandomForest (fit)", cpu, t.dt)
globals()["_RF_ARTIFACTS"] = (cu_rf, Xr_g, yr_g, Xr_c, yr_c)
print(f"DBSCAN ({N_DBSCAN:,} x 8)")
Xd_g, _ = gpu_make_blobs(n_samples=N_DBSCAN, n_features=8, centers=6,
cluster_std=0.6, random_state=SEED,
dtype=np.float32)
Xd_c = cp.asnumpy(Xd_g)
with Timer("sklearn", sync=False) as t:
lab_c = skDBSCAN(eps=0.9, min_samples=8, n_jobs=-1).fit_predict(Xd_c)
cpu = t.dt
with Timer("cuML") as t:
lab_g = cuDBSCAN(eps=0.9, min_samples=8).fit_predict(Xd_g)
record("DBSCAN", cpu, t.dt)
print(f" clusters found: sklearn={len(set(lab_c.tolist())) - 1}, "
f"cuML={len(set(cp.asnumpy(lab_g).tolist())) - 1}\n")
We benchmark scikit-learn and cuML implementations of PCA, K-Means, nearest neighbors, logistic regression, random forests, and DBSCAN. We generate datasets on the GPU, synchronize CUDA operations for fair timing, and record the speedup each accelerated algorithm achieves. We also compare model behavior and retain the trained cuML random forest so that we can reuse it later in the tutorial.
def demo_manifold():
from cuml.manifold import UMAP, TSNE
from cuml.metrics import trustworthiness
X, y = gpu_make_blobs(n_samples=N_MANIFOLD, n_features=48, centers=8,
cluster_std=1.6, random_state=SEED, dtype=np.float32)
print(f"data: {X.shape}")
embeddings = {}
for n_neighbors, min_dist in [(15, 0.1), (50, 0.0)]:
key = f"UMAP(n_neighbors={n_neighbors}, min_dist={min_dist})"
with Timer(key) as t:
emb = UMAP(n_neighbors=n_neighbors, min_dist=min_dist,
n_components=2, random_state=SEED).fit_transform(X)
sub = slice(0, min(5000, X.shape[0]))
tw = trustworthiness(X[sub], emb[sub], n_neighbors=10)
print(f" trustworthiness = {tw:.4f}")
embeddings[key] = (emb, t.dt, tw)
with Timer("TSNE(method='fft')") as t:
tsne_emb = TSNE(n_components=2, perplexity=30,
random_state=SEED).fit_transform(X)
embeddings["TSNE"] = (tsne_emb, t.dt, float("nan"))
best_key = max([k for k in embeddings if k.startswith("UMAP")],
key=lambda k: embeddings[k][2])
emb = embeddings[best_key][0]
print(f"\nClustering the '{best_key}' embedding with GPU HDBSCAN")
try:
from cuml.cluster import HDBSCAN
with Timer("HDBSCAN") as t:
hdb = HDBSCAN(min_cluster_size=max(int(50 * SCALE), 5),
min_samples=10, prediction_data=True).fit(emb)
labels = cp.asarray(hdb.labels_)
n_clusters = int(labels.max()) + 1
noise = float((labels == -1).mean())
print(f" clusters={n_clusters} noise fraction={noise:.3f}")
try:
from cuml.metrics.cluster import adjusted_rand_score
print(f" adjusted Rand index vs ground truth: "
f"{adjusted_rand_score(y, labels):.4f}")
except Exception as e:
print(f" ARI skipped: {e}")
try:
from cuml.cluster.hdbscan import all_points_membership_vectors
mv = all_points_membership_vectors(hdb)
print(f" soft-cluster membership matrix: {tuple(mv.shape)}")
except Exception as e:
print(f" soft clustering skipped: {e}")
except Exception as e:
print(f" HDBSCAN step skipped: {e}")
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
keys = list(embeddings)[:3]
for ax, k in zip(axes, keys):
e = to_numpy(embeddings[k][0])
c = to_numpy(y)
ax.scatter(e[:, 0], e[:, 1], c=c, s=1.5, cmap="tab10", alpha=0.6)
ax.set_title(f"{k}\n{embeddings[k][1]:.2f}s", fontsize=9)
ax.set_xticks([]); ax.set_yticks([])
plt.suptitle("GPU manifold learning (colored by ground-truth cluster)")
plt.tight_layout(); plt.show()
We build an unsupervised GPU pipeline using UMAP and t-SNE to reduc