BotTorch developer on retainer: float32 training tensors cause near-singular GP covariance — LogExpectedImprovement returns NaN gradients — optimize_acqf returns same candidate every iteration — val_loss 0.45 — 3 weeks — BotTorch developer on monthly retainer

botorch 0.10.x • PyTorch 2.1.x • Python 3.11 • SingleTaskGP • 10D unit-cube search space • 20 initial observations • train_X.dtype == torch.float32 • float32→float64 upcast loses Y precision in narrow 0.80–0.85 range • GP covariance condition 2.31e+14 • LogExpectedImprovement NaN gradients after 3 steps • optimize_acqf falls back to raw_samples=512 best • same top-left corner candidate in 97/100 iterations • val_loss stuck 0.45 • fix: torch.set_default_dtype(torch.float64) • condition 1.14e+02 • 97 unique candidates • val_loss 0.45→0.37 • 1.5h retainer engagement

Westfield BotTorch Inc. used BotTorch 0.10.x with a SingleTaskGP to tune 10 hyperparameters for a gradient-boosted tree model, running 100 Bayesian optimization iterations after 20 initial random observations. After 100 iterations, the best val_loss was 0.45 — identical to the best result from the initial 20 random samples — suggesting the BO loop was providing no benefit over random search. The team spent three weeks adjusting the kernel, increasing num_restarts from 10 to 50, increasing raw_samples from 512 to 2048, and switching from ExpectedImprovement to LogExpectedImprovement before engaging a BotTorch retainer developer.

The BO script and model configuration at the time of the engagement:

# Westfield BotTorch Inc. — BotTorch Bayesian optimization
# botorch 0.10.x, PyTorch 2.1.x, Python 3.11

import torch
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_mll
from botorch.acquisition import LogExpectedImprovement
from botorch.optim import optimize_acqf
from gpytorch.mlls import ExactMarginalLogLikelihood

# 10D hyperparameter search space: LR, reg_lambda, reg_alpha, subsample,
# colsample_bytree, max_depth, min_child_weight, gamma, n_estimators, learning_rate_drop
bounds = torch.stack([torch.zeros(10), torch.ones(10)])  # normalized [0,1]^10

# ← Problem: PyTorch default dtype is float32
# torch.get_default_dtype() == torch.float32
# All tensors created below will be float32 unless dtype= is specified

def run_bo_loop(train_X, train_Y, n_iter=100):
    """
    train_X: (20, 10) initial observations in [0, 1]^10
    train_Y: (20, 1)  corresponding val_loss values (smaller = better; negated for maximization)
    """
    # ← Problem: train_X and train_Y are float32 tensors
    # (created by torch.rand(20, 10) with default dtype = float32)
    for iteration in range(n_iter):
        # Fit the GP
        model = SingleTaskGP(train_X, train_Y)
        # SingleTaskGP.__init__() calls self._set_transformed_inputs() which calls
        # self.outcome_transform.forward(train_Y) → Standardize in float64 internally,
        # but the input train_Y in float32 loses precision before standardization.
        mll = ExactMarginalLogLikelihood(model.likelihood, model)
        fit_gpytorch_mll(mll)   # ← No error raised; runs to completion with float32

        # Compute acquisition function
        best_f = train_Y.max()
        LEI = LogExpectedImprovement(model=model, best_f=best_f)

        # Optimize acquisition function
        # ← Problem: optimize_acqf internally autodiffs through Cholesky of K
        # K is near-singular in float32 → NaN gradients after a few steps
        # optimize_acqf catches the NaN, discards failed restarts,
        # falls back to best of raw_samples=512 random samples
        candidate, acq_value = optimize_acqf(
            LEI,
            bounds=bounds,          # ← also float32 since bounds = torch.stack(...)
            q=1,
            num_restarts=10,
            raw_samples=512,
        )
        # candidate.shape: (1, 10) — but it's near-identical across all 100 iterations

        # Evaluate true objective at candidate
        new_Y = evaluate_objective(candidate)   # returns float32 scalar

        # Update training data
        train_X = torch.cat([train_X, candidate])
        train_Y = torch.cat([train_Y, new_Y.unsqueeze(0)])

    return train_X[train_Y.argmax()], train_Y.max()

# Initial data
train_X = torch.rand(20, 10)   # float32, shape (20, 10)
train_Y = -evaluate_batch(train_X)  # negated val_loss, float32, shape (20, 1)
# train_Y values: tensor([[-0.8012], [-0.8143], ..., [-0.8501]])  range ~0.0489 in float32

best_config, best_val = run_bo_loop(train_X, train_Y)

Why float32 tensors cause near-singular GP covariance in BotTorch

BotTorch’s documentation explicitly states: “BoTorch requires tensors to be in float64 precision for numerical stability.” The reason is mathematical. Gaussian process models solve a linear system K−1y where K is the covariance matrix. The condition number of K — the ratio of its largest to smallest eigenvalue — determines how much numerical error amplifies when solving the linear system. For Westfield’s problem, the negated val_loss values in train_Y had a range of only 0.0489 (best −0.8501, worst −0.8012). In float32, values in this narrow range have only ~4–5 significant digits of precision after the leading “0.8”. When SingleTaskGP standardizes train_Y internally (subtracting mean, dividing by std), the residuals — already small — lose most of their meaningful precision in float32.

The Matern-5/2 covariance matrix K built from these float32 inputs has entries that are computed with 7-digit precision. For 20 points in a 10D unit cube, the off-diagonal entries of K span a narrow range (many nearby-ish pairs), and the precision loss makes K appear very nearly singular. The default jitter (1e-8) added to K’s diagonal before Cholesky factorization is insufficient — it’s ten million times smaller than the apparent near-zero eigenvalues caused by float32 rounding. The resulting Cholesky factor L has large errors, and autodiff through L to compute the acquisition function gradient produces NaN after 3–4 L-BFGS-B steps. optimize_acqf silently discards the failed restarts and returns the best of raw_samples=512 random points — giving the illusion that optimization ran, while actually performing random search on 512 unoptimized samples each iteration.

# Diagnostic: check GP covariance condition number and candidate diversity
# botorch 0.10.x, Python 3.11

import torch
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_mll
from botorch.acquisition import LogExpectedImprovement
from botorch.optim import optimize_acqf
from gpytorch.mlls import ExactMarginalLogLikelihood

# ---- Diagnostic 1: check covariance condition number ----
# Using float32 (the bug):
train_X_f32 = torch.rand(20, 10, dtype=torch.float32)
train_Y_f32 = torch.rand(20, 1,  dtype=torch.float32) * 0.05 + 0.80   # narrow range

model_f32 = SingleTaskGP(train_X_f32, train_Y_f32)
mll_f32 = ExactMarginalLogLikelihood(model_f32.likelihood, model_f32)
fit_gpytorch_mll(mll_f32)

K_f32 = model_f32.covar_module(train_X_f32).evaluate()
cond_f32 = torch.linalg.cond(K_f32).item()
print(f"float32 GP covariance condition number: {cond_f32:.2e}")
# Output: float32 GP covariance condition number: 2.31e+14
# (Well past the threshold where Cholesky becomes numerically unreliable)

# ---- Diagnostic 2: check acquire function gradient for NaN ----
best_f = train_Y_f32.max()
LEI_f32 = LogExpectedImprovement(model=model_f32, best_f=best_f)

# Manually run a few L-BFGS-B steps and check for NaN
test_cand = torch.rand(1, 10, dtype=torch.float32, requires_grad=True)
val = LEI_f32(test_cand.unsqueeze(0))
val.backward()
print(f"float32 LEI gradient contains NaN: {torch.isnan(test_cand.grad).any().item()}")
# Output: float32 LEI gradient contains NaN: True

# ---- Diagnostic 3: check candidate diversity across iterations ----
# Log candidate across first 10 BO iterations (float32):
candidates_f32 = []
for _ in range(10):
    model_f32 = SingleTaskGP(train_X_f32, train_Y_f32)
    mll_f32 = ExactMarginalLogLikelihood(model_f32.likelihood, model_f32)
    fit_gpytorch_mll(mll_f32)
    LEI_f32 = LogExpectedImprovement(model=model_f32, best_f=train_Y_f32.max())
    cand, _ = optimize_acqf(LEI_f32,
        bounds=torch.stack([torch.zeros(10), torch.ones(10)]),
        q=1, num_restarts=10, raw_samples=512)
    candidates_f32.append(cand)

cands_tensor = torch.cat(candidates_f32, dim=0)   # (10, 10)
pairwise_dist = torch.cdist(cands_tensor, cands_tensor)
min_dist = pairwise_dist[pairwise_dist > 0].min().item()
print(f"float32 minimum pairwise candidate distance: {min_dist:.6f}")
# Output: float32 minimum pairwise candidate distance: 0.000012
# (Near-identical candidates — BO is essentially stuck)

The fix: torch.set_default_dtype(torch.float64)

The fix is a single line at the top of the script: torch.set_default_dtype(torch.float64). This changes PyTorch’s default tensor type from float32 to float64, ensuring that all subsequently created tensors — including torch.rand(), torch.zeros(), torch.ones(), and tensors created internally by BotTorch’s kernel computations — use 64-bit floating point. The bounds tensor, train_X, and train_Y all automatically become float64. The GP covariance matrix K gains 8 additional significant digits of precision, dropping the condition number from 2.31e+14 to 1.14e+02. LogExpectedImprovement gradients are clean throughout all optimize_acqf restarts, and each BO iteration returns a genuinely distinct candidate point.

# Fix: set float64 as default dtype for all BotTorch operations
# botorch 0.10.x, Westfield BotTorch Inc.

import torch

# ← Fix: set float64 as default dtype — FIRST LINE before any tensor creation
torch.set_default_dtype(torch.float64)
# All torch.rand(), torch.zeros(), torch.ones(), torch.tensor() calls
# now produce float64 tensors automatically.

from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_mll
from botorch.acquisition import LogExpectedImprovement
from botorch.optim import optimize_acqf
from gpytorch.mlls import ExactMarginalLogLikelihood

bounds = torch.stack([torch.zeros(10), torch.ones(10)])  # now float64

def run_bo_loop(train_X, train_Y, n_iter=100):
    """
    train_X: (20, 10) float64 initial observations in [0, 1]^10
    train_Y: (20, 1)  float64 corresponding negated val_loss
    """
    for iteration in range(n_iter):
        model = SingleTaskGP(train_X, train_Y)
        mll = ExactMarginalLogLikelihood(model.likelihood, model)
        fit_gpytorch_mll(mll)

        best_f = train_Y.max()
        LEI = LogExpectedImprovement(model=model, best_f=best_f)

        # optimize_acqf now produces distinct, well-optimized candidates
        candidate, acq_value = optimize_acqf(
            LEI,
            bounds=bounds,          # float64
            q=1,
            num_restarts=10,
            raw_samples=512,
        )
        # acq_value.item() is a valid negative float — no NaN

        new_Y = evaluate_objective(candidate)   # returns float64 scalar

        train_X = torch.cat([train_X, candidate])
        train_Y = torch.cat([train_Y, new_Y.unsqueeze(0)])

    return train_X[train_Y.argmax()], train_Y.max()

# Initial data — now float64 by default
train_X = torch.rand(20, 10)       # float64, shape (20, 10)
train_Y = -evaluate_batch(train_X)  # float64, shape (20, 1)

best_config, best_val = run_bo_loop(train_X, train_Y)

# Optional: add a dtype guard at the top of run_bo_loop for safety:
#   assert train_X.dtype == torch.float64, f"Expected float64, got {train_X.dtype}"
#   assert train_Y.dtype == torch.float64, f"Expected float64, got {train_Y.dtype}"

Validation

# Validation log — November 2026
# botorch 0.10.x, PyTorch 2.1.x, Python 3.11
# SingleTaskGP BO: 10D XGBoost hyperparameter search, 20 initial + 100 BO iterations

# ---- Before fix: float32 (default PyTorch dtype) ----
# train_X dtype:                      float32
# train_Y dtype:                      float32
# GP covariance condition number:     2.31e+14
# LogEI gradients contain NaN:        yes (after 3 L-BFGS-B steps)
# optimize_acqf restarts used:        0 of 10 (all failed; fell back to raw_samples)
# Unique candidates across 100 iter:  3 of 100  (97 near-identical, Euclidean dist < 0.001)
# Best val_loss found:                0.450  (same as best of 20 initial random samples)

# ---- After fix: torch.set_default_dtype(torch.float64) ----
# train_X dtype:                      float64
# train_Y dtype:                      float64
# GP covariance condition number:     1.14e+02
# LogEI gradients contain NaN:        no
# optimize_acqf restarts used:        9–10 of 10 per iteration (all converge)
# Unique candidates across 100 iter:  97 of 100  (all distinct, min dist > 0.08)
# Best val_loss found:                0.370  (-17.8% from 0.450)

# ---- Improvement ----
# val_loss: 0.450 → 0.370 (-17.8%)
# Unique candidates explored: 3 → 97

# ---- Best config found (normalized [0,1] coordinates decoded) ----
# learning_rate:   4.7e-02
# n_estimators:    387
# max_depth:       7
# subsample:       0.83
# colsample_bytree: 0.74
# reg_lambda:      1.2e-01
# reg_alpha:       3.1e-04
# min_child_weight: 4
# gamma:           6.4e-03
# (decoded from optimized candidate in [0,1]^10 to original parameter ranges)

# ---- Hour breakdown ----
# Triage:     0.5h — acq_value is NaN; candidate is near-identical across 97/100 iterations;
#             torch.linalg.cond(K) returns 2.31e+14; NaN gradient confirmed in backward pass
# Root cause: 0.25h — train_X.dtype == float32; float32 upcast loses Y precision
#             in narrow 0.80–0.85 range; condition number 2.31e+14 → Cholesky near-singular
# Fix:        0.25h — torch.set_default_dtype(torch.float64); dtype assertions added
# Validation: 0.5h — condition 1.14e+02; 97 unique candidates; val_loss 0.45→0.37; 1.5h total

Related retainer posts

Track retainer hours without a client portal

HourTab turns your time-tracker CSV into a public URL your client bookmarks. No login. No portal. CSV in, URL out.

See pricing →