tidymodels workflow with psvr

This vignette demonstrates the full tidymodels pipeline with psvr: data splitting, preprocessing, hyperparameter tuning by cross-validation, and final model evaluation. We use psvr_rmspe_rbf() (LS-SVR with RMSPE loss, RBF kernel) and tune the regularisation parameter cost ( $\Gamma$ ) against MAPE.

library(psvr)
library(parsnip)
library(rsample)
library(recipes)
library(workflows)
library(tune)
library(dials)
library(yardstick)
library(dplyr)

Data

The synthetic even-function dataset from the package README: $y = 2 + x_1^2 + 0.5\,x_2^2 + \varepsilon$ , $\varepsilon \sim \mathcal{N}(0,\,0.1^2)$ . Targets are strictly positive by construction ( $y > 0$ ).

set.seed(42)
n   <- 200
x1  <- runif(n, -3, 3)
x2  <- runif(n, -3, 3)
y   <- 2 + x1^2 + 0.5 * x2^2 + rnorm(n, sd = 0.1)
dat <- data.frame(y = y, x1 = x1, x2 = x2)

1 — Split

set.seed(1)
split <- initial_split(dat, prop = 0.75)
train <- training(split)
test  <- testing(split)

2 — Preprocessing recipe

Centre and scale all predictors so the RBF kernel operates on a standardised feature space.

rec <- recipe(y ~ x1 + x2, data = train) |>
  step_normalize(all_predictors())

3 — Model spec with `tune()`

Both cost (maps to $\Gamma$ ) and rbf_sigma (the RBF bandwidth $\sigma$ ) are tune() placeholders; the grid search will explore all combinations.

spec <- psvr_rmspe_rbf(cost = tune(), rbf_sigma = tune()) |>
  set_engine("psvr")

4 — Workflow

wf <- workflow() |>
  add_recipe(rec) |>
  add_model(spec)

5 — Tune with 5-fold CV

We search over a 15-point Latin hypercube of cost and rbf_sigma values and evaluate each fold by MAPE. Both search ranges are set from the data, and neither happens automatically — dials cannot finalize either one, so both have to be passed explicitly through param_info:

rbf_sigma_psvr_data() centres the bandwidth range on the median pairwise distance, so it is computed on the baked predictors: the heuristic only means anything on the scale the model is actually fitted on.
cost_psvr_ls_data() widens the cost range. Here cost is $\Gamma$ , and its registered default of $[-2, 10]$ on the log2 scale ( $\Gamma \le 1024$ ) is the $\epsilon$ -SVR range — far too low for LS-SVR, where the optimum scales with var(y) * n. Left at the default the grid tops out at $\Gamma = 1024$ , which on this dataset is an order of magnitude below the value selected once the range is widened — compare the cost column printed below. The search is boundary-trapped: it cannot reach the optimum at all, and nothing warns you, because every candidate it did evaluate was legal. This one cannot be automated even in principle, because tune finalizes parameters from the predictors alone and never passes the outcome to dials::finalize().

The RMSPE LS-SVR only solves an (N+1) × (N+1) linear system — no iterative solver is involved — so 75 fits complete in seconds.

set.seed(2)
folds <- vfold_cv(train, v = 5)

# Data-driven rbf_sigma range centred on median pairwise distance
train_baked      <- rec |> prep() |> bake(new_data = train)
rbf_sigma_custom <- rbf_sigma_psvr_data(train_baked |> select(-y))

wf_params <- extract_parameter_set_dials(wf) |>
  update(
    cost      = cost_psvr_ls_data(train$y),
    rbf_sigma = rbf_sigma_custom
  )

tune_res <- tune_grid(
  wf,
  resamples  = folds,
  grid       = 15,
  param_info = wf_params,
  metrics    = metric_set(yardstick::mape)
)

Cross-validated MAPE for each candidate (lower is better):

collect_metrics(tune_res)[, c("cost", "rbf_sigma", "mean", "std_err")]
#> # A tibble: 15 × 4
#>         cost rbf_sigma  mean std_err
#>        <dbl>     <dbl> <dbl>   <dbl>
#>  1     0.25      4.78  42.0   1.07  
#>  2     0.562     0.478 33.8   1.12  
#>  3     1.26      1.78  36.4   1.13  
#>  4     2.84     12.8   42.0   1.07  
#>  5     6.39      0.178 29.6   1.03  
#>  6    14.4       0.665 10.7   0.244 
#>  7    32.3       2.48  18.2   0.863 
#>  8    72.6       9.23  41.1   1.26  
#>  9   163.        0.248  7.83  0.419 
#> 10   367.        0.923  2.97  0.0925
#> 11   825.        3.44   4.51  0.231 
#> 12  1854.       17.8   40.3   1.37  
#> 13  4169.        0.344  3.78  0.442 
#> 14  9372.        1.28   1.78  0.0907
#> 15 21070.        6.65   2.64  0.221

6 — Select best

best_params <- select_best(tune_res, metric = "mape")
best_params
#> # A tibble: 1 × 3
#>    cost rbf_sigma .config         
#>   <dbl>     <dbl> <chr>           
#> 1 9372.      1.28 pre0_mod14_post0

7 — Final fit and test-set evaluation

last_fit() refits on the full training set with the chosen cost and evaluates once on the held-out test data.

final_wf  <- finalize_workflow(wf, best_params)
final_fit <- last_fit(final_wf, split, metrics = metric_set(yardstick::mape))

collect_metrics(final_fit)
#> # A tibble: 1 × 4
#>   .metric .estimator .estimate .config        
#>   <chr>   <chr>          <dbl> <chr>          
#> 1 mape    standard        2.03 pre0_mod0_post0

Predictions on the test set:

preds <- collect_predictions(final_fit)
head(preds[, c(".row", "y", ".pred")])
#> # A tibble: 6 × 3
#>    .row     y .pred
#>   <int> <dbl> <dbl>
#> 1     3  5.99  5.91
#> 2     4  6.20  6.04
#> 3     5  4.69  4.84
#> 4     6  1.96  2.08
#> 5     8  6.70  6.81
#> 6     9  3.93  3.97

The fitted workflow can also be used directly for new data:

new_obs <- data.frame(x1 = c(0, 1, -2), x2 = c(0, 1, 2))
predict(extract_workflow(final_fit), new_data = new_obs)
#> # A tibble: 3 × 1
#>   .pred
#>   <dbl>
#> 1  2.03
#> 2  3.45
#> 3  8.10

8 — Inspecting the fitted psvr model

The tidymodels layer wraps a psvr_rmspe object (returned by the engine fit wrapper psvr_rmspe_rbf_fit()). Extract it to use print() and coef() directly.

# extract_fit_engine() unwraps the parsnip/workflow layer to the raw psvr
# object -- the same class psvr_rmspe() returns when called directly
engine_fit <- extract_fit_engine(extract_workflow(final_fit))
print(engine_fit)
#> 
#> LS-SVR with RMSPE loss  [psvr_rmspe]
#> 
#>   Kernel:        RBF (sigma = 1.28306)
#>   Gamma:         9371.97
#>   Training obs.: 150

cf <- coef(engine_fit)
# alpha:        N dual variables; weight each training point in
#               f(x) = sum_k alpha_k K(x_k, x) + b
# b:            bias / intercept term
# support_data: all N training inputs (LS-SVR has no sparsity — every training
#               point contributes, so despite the name this is not a subset)
cat(sprintf("b = %.4f  |  alpha range: [%.4f, %.4f]\n",
            cf$b, min(cf$alpha), max(cf$alpha)))
#> b = 14.8943  |  alpha range: [-366.4647, 242.8288]