Kernel DRO¶
Kernels enter this package in two different ways. In MMD-DRO, a kernel defines the distance between the empirical distribution and distributions in the ambiguity set. Other examples below use kernels to construct a nonlinear prediction model while retaining a standard DRO ambiguity set.
Choosing an ambiguity-set kernel for MMD-DRO¶
MMD_DRO uses the RBF kernel by default. RBF and Laplacian kernels are characteristic, so their MMD distinguishes probability distributions under standard conditions. The polynomial kernel emphasizes interactions up to a selected degree, while the linear kernel compares first moments only.
For RBF and Laplacian kernels, leaving kernel_gamma=None selects a robust median-distance scale. For a polynomial kernel, the automatic default is gamma = 1 / p, where p is the dimension of the joint sample (X, y); the default degree is 2 and the default offset is 1.0. Small degrees such as 2 or 3 are usually the most stable starting points. Standardize continuous variables before comparing kernels, particularly when using a polynomial kernel.
[1]:
from dro.linear_model.mmd_dro import MMD_DRO
from dro.data.dataloader_classification import classification_basic
feature_dim = 2
X_mmd, y_mmd = classification_basic(
d=feature_dim, num_samples=100, radius=3, visualize=False
)
# Omitting kernel selects the default RBF kernel and median heuristic.
rbf_mmd = MMD_DRO(input_dim=feature_dim, model_type='svm')
# An inhomogeneous quadratic polynomial kernel. gamma=None gives 1 / (d + 1).
polynomial_mmd = MMD_DRO(
input_dim=feature_dim,
model_type='svm',
kernel='polynomial',
kernel_gamma=None,
kernel_degree=2,
kernel_coef0=1.0,
)
polynomial_mmd.update({'eta': 0.1, 'sampling_method': 'hull'})
mmd_params = polynomial_mmd.fit(X_mmd, y_mmd)
(CVXPY) Oct 24 05:35:25 PM: Your problem has 104 variables, 200 constraints, and 0 parameters.
(CVXPY) Oct 24 05:35:25 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Oct 24 05:35:25 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Oct 24 05:35:25 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
(CVXPY) Oct 24 05:35:25 PM: Your problem is compiled with the CPP canonicalization backend.
(CVXPY) Oct 24 05:35:25 PM: Compiling problem (target solver=CLARABEL).
(CVXPY) Oct 24 05:35:25 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> CLARABEL
(CVXPY) Oct 24 05:35:25 PM: Applying reduction Dcp2Cone
(CVXPY) Oct 24 05:35:25 PM: Applying reduction CvxAttr2Constr
(CVXPY) Oct 24 05:35:25 PM: Applying reduction ConeMatrixStuffing
(CVXPY) Oct 24 05:35:25 PM: Applying reduction CLARABEL
(CVXPY) Oct 24 05:35:25 PM: Finished problem compilation (took 3.301e-02 seconds).
(CVXPY) Oct 24 05:35:25 PM: Invoking solver CLARABEL to obtain a solution.
===============================================================================
CVXPY
v1.6.7
===============================================================================
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
-------------------------------------------------------------
Clarabel.rs v0.11.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 305
constraints = 701
nnz(P) = 0
nnz(A) = 21501
cones (total) = 2
: Nonnegative = 1, numel = 600
: SecondOrder = 1, numel = 101
settings:
linear algebra: direct / faer, precision: 64 bit (10 threads)
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-4, max_scale = 1.0e4
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
(CVXPY) Oct 24 05:35:25 PM: Problem status: optimal
(CVXPY) Oct 24 05:35:25 PM: Optimal value: 3.338e-01
(CVXPY) Oct 24 05:35:25 PM: Compilation took 3.301e-02 seconds
0 +1.0810e-01 +1.7921e+02 1.79e+02 8.68e-01 1.39e+01 1.00e+00 1.72e+00 ------
1 +3.1976e+00 +2.1047e+01 5.58e+00 6.93e-02 1.50e+00 5.64e-02 1.91e-01 9.08e-01
2 +1.4128e+00 +4.8852e+00 2.46e+00 2.29e-02 4.51e-01 1.37e-02 7.39e-02 6.73e-01
3 +4.9034e-01 +1.3326e+00 8.42e-01 7.04e-03 1.31e-01 6.80e-03 3.06e-02 9.58e-01
4 +4.1436e-01 +7.7251e-01 3.58e-01 3.08e-03 5.77e-02 2.18e-03 1.46e-02 7.61e-01
5 +3.7566e-01 +5.3947e-01 1.64e-01 1.46e-03 2.72e-02 8.29e-04 7.46e-03 6.29e-01
6 +3.5248e-01 +4.2192e-01 6.94e-02 6.40e-04 1.19e-02 2.54e-04 3.51e-03 7.96e-01
7 +3.5130e-01 +4.1238e-01 6.11e-02 5.61e-04 1.04e-02 2.00e-04 3.08e-03 3.59e-01
8 +3.3569e-01 +3.4889e-01 1.32e-02 1.23e-04 2.29e-03 3.70e-05 7.12e-04 8.67e-01
9 +3.3488e-01 +3.4092e-01 6.04e-03 5.64e-05 1.04e-03 1.42e-05 3.27e-04 7.57e-01
10 +3.3390e-01 +3.3525e-01 1.36e-03 1.27e-05 2.35e-04 3.01e-06 7.39e-05 8.24e-01
11 +3.3385e-01 +3.3477e-01 9.21e-04 8.62e-06 1.60e-04 1.91e-06 5.02e-05 5.19e-01
12 +3.3379e-01 +3.3406e-01 2.72e-04 2.55e-06 4.70e-05 5.51e-07 1.49e-05 7.35e-01
13 +3.3376e-01 +3.3380e-01 3.63e-05 3.41e-07 6.29e-06 6.82e-08 1.99e-06 9.29e-01
14 +3.3376e-01 +3.3377e-01 4.05e-06 3.80e-08 7.01e-07 7.36e-09 2.22e-07 9.27e-01
15 +3.3376e-01 +3.3376e-01 5.26e-08 4.94e-10 9.11e-09 9.56e-11 2.88e-09 9.87e-01
16 +3.3376e-01 +3.3376e-01 6.37e-09 5.97e-11 1.10e-09 1.15e-11 3.48e-10 8.91e-01
---------------------------------------------------------------------------------------------
Terminated with status = Solved
solve time = 98.810498ms
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) Oct 24 05:35:25 PM: Solver (including time spent in interface) took 1.036e-01 seconds
Kernel predictors for Wasserstein and f-divergence DRO¶
Here the kernel constructs a nonlinear predictor; the ambiguity set still follows the selected DRO method.
The cubic polynomial kernel on two input features has ten independent feature directions. We use ten Nyström components to avoid a redundant 500-coefficient optimization problem and keep the example numerically stable with CLARABEL. For this dataset, the reduced representation reconstructs the polynomial kernel to numerical precision. We explicitly choose p=2 and eps=0.01 for the Wasserstein model.
The chi-square, KL, CVaR, and TV examples also use ten-component representations for numerical stability. The sigmoid-kernel CVaR logistic model uses SCS; the other models use CLARABEL. These are reduced kernel representations; tune their component count and robustness parameters for your own data.
[2]:
from dro.linear_model.wasserstein_dro import WassersteinDRO
from dro.data.dataloader_classification import classification_basic
feature_dim = 2
X, y = classification_basic(
d=feature_dim, num_samples=500, radius=3, seed=42, visualize=False
)
kernel_clf_model = WassersteinDRO(
input_dim=feature_dim, model_type='svm', kernel='poly', solver='CLARABEL'
)
kernel_clf_model.update_kernel({
'metric': 'poly', 'kernel_gamma': 1.0, 'n_components': 10,
})
kernel_clf_model.update({'eps': 0.01, 'p': 2})
params = kernel_clf_model.fit(X, y)
print(f'Robust training objective: {kernel_clf_model.robust_obj:.6f}')
accuracy, f1 = kernel_clf_model.score(X, y)
print(f'Training accuracy: {accuracy:.3f}; F1: {f1:.3f}')
Robust training objective: 0.372083
Training accuracy: 0.984; F1: 0.984
[3]:
import numpy as np
from dro.linear_model.chi2_dro import *
from dro.linear_model.kl_dro import *
from dro.linear_model.tv_dro import *
from dro.linear_model.cvar_dro import *
from dro.linear_model.marginal_dro import *
from dro.linear_model.conditional_dro import *
np.random.seed(42)
kernel_clf_model = Chi2DRO(input_dim = feature_dim, model_type = 'svm', kernel = 'poly')
kernel_clf_model.update_kernel({'n_components': 10})
kernel_clf_model.update({'eps': 10})
params = kernel_clf_model.fit(X, y)
kernel_clf_model = KLDRO(input_dim = feature_dim, model_type = 'svm', kernel = 'rbf')
kernel_clf_model.update_kernel({'n_components': 10})
kernel_clf_model.update({'eps': 2})
params = kernel_clf_model.fit(X, y)
kernel_clf_model = CVaRDRO(
input_dim=feature_dim, model_type='logistic', kernel='sigmoid', solver='SCS'
)
kernel_clf_model.update_kernel({'n_components': 10})
kernel_clf_model.update({'alpha': 0.9})
params = kernel_clf_model.fit(X, y)
kernel_clf_model = TVDRO(input_dim = feature_dim, model_type = 'logistic', kernel = 'poly')
kernel_clf_model.update_kernel({'n_components': 10})
kernel_clf_model.update({'eps': 0.01})
params = kernel_clf_model.fit(X, y)
[4]:
kernel_clf_model = MarginalCVaRDRO(input_dim = feature_dim, model_type = 'svm', kernel = 'poly')
kernel_clf_model.update({'alpha': 0.9})
params = kernel_clf_model.fit(X, y)
kernel_clf_model.score(X,y)
[4]:
(0.986, 0.9859994959818554)
[5]:
from sklearn.svm import SVC
from dro.data.dataloader_regression import regression_basic
from dro.data.dataloader_classification import classification_basic
from dro.data.draw_utils import draw_classification
feature_dim = 5
X, y = classification_basic(d = feature_dim, num_samples = 25, radius = 3, visualize = False)
clf = SVC(kernel = 'rbf', C = 100000000)
clf.fit(X, y)
print('coef', clf.dual_coef_, clf.support_vectors_, clf.intercept_)
coef [[-0.25769266 -0.64158585 -3.57086163 -0.78030986 -0.43010956 -0.01490551
0.76823835 0.16369363 1.56160661 1.47227558 0.03362333 1.69602757]] [[ 0.38846926 -0.7028592 1.35277704 0.69878398 -2.12650133]
[ 2.31753572 -0.46290575 1.17834297 1.18731603 -0.94596622]
[ 0.25018034 1.61514873 1.09731754 1.55435329 0.42096142]
[-0.3689567 -0.02826586 -0.84885536 1.28387817 -0.20472226]
[ 0.13204275 -0.69776822 2.16793699 2.95568251 -2.16462365]
[ 1.78316707 -1.07634697 0.80160239 2.94332765 0.57396163]
[-0.00613887 4.01699229 1.16884802 0.82769028 -1.76806615]
[ 0.45436803 2.56600279 0.90566671 -0.6451926 -1.13588996]
[-2.28630569 2.45244185 1.26490427 1.72628874 0.65931799]
[-0.06598713 2.44424396 0.03599602 0.40586944 1.60361199]
[ 0.42349747 1.56956827 2.60746837 -2.13880444 1.43853605]
[ 1.82292115 1.4884194 0.63837633 -0.63730201 0.3482033 ]] [0.47099307]
[6]:
feature_dim = 5
X, y = regression_basic(num_samples = 20, d = feature_dim, noise = 1)
reg_model = WassersteinDRO(input_dim = feature_dim, model_type = 'lad', kernel = 'rbf')
reg_model.update_kernel({'metric': 'rbf'})
reg_model.update({'eps': 1})
reg_model.fit(X, y)
[6]:
{'theta': [28.183840216059043,
-40.587159870262745,
-25.771702054545074,
24.507923542613252,
-20.773407535187975,
-32.971991106392935,
3.7200730037962355,
-71.2804779091182,
70.9566258049175,
-50.902349774153635,
41.142793433932134,
-57.15251168586913,
57.59232483199867,
48.59608462873434,
26.428005719050088,
-0.056328121150369746,
9.225576771876137,
-40.70582945485824,
-30.048318322329195,
19.958741698253867],
'b': array(-37.78934524)}