バイナリ分類のサポートベクターマシン

サポートベクターマシンについて

可分データ
不可分データ
カーネルによる非線形変換

可分データ

データに 2 つのクラスのみ含まれる場合は、サポートベクターマシン (SVM) を使用することができます。SVM では、1 つのクラスのすべてのデータ点を他のクラスのすべてのデータ点から分離する最適な超平面を見つけてデータを分類します。SVM の "最適な" 超平面とは、2 つのクラス間に最大の "マージン" をもつ超平面のことです。マージンとは、内部にデータ点のない、超平面に平行するスラブ平面の最大幅を意味します。

"サポートベクター" は分離超平面に最も近いデータ点です。これらの点はスラブの境界上にあります。次の図は、これらの定義を示しています。+ はタイプ 1 のデータ点、- はタイプ -1 のデータ点を示しています。

数学的形成化: 主問題- この説明は、Hastie、Tibshirani および Friedman [1] と Christianini および Shawe-Taylor [2] に従います。

学習のデータは、点集合 (ベクトル) x_j とそれらのカテゴリ y_j です。次元 d について、x_j ∊ R^d および y_j = ±1 です。超平面の方程式は、次のようになります。

$f (x) = x' β + b = 0$

ここで、β ∊ R^d と b は実数です。

次の問題は、"最適な" 分離超平面 (つまり、判定境界) を定義します。すべてのデータ点 (x_j,y_j) で次のようになるように、||β|| を最小にする β および b を見つけます。

$y_{j} f (x_{j}) \geq 1.$

サポートベクターは、 $y_{j} f (x_{j}) = 1.$ である、境界における x_j です。

数学的な便宜上、通常この問題は $‖ β ‖$ を最小にする同等の問題として与えられます。これは、二次計画法の問題です。最適解 $(\hat{β}, \hat{b})$ により、ベクトル z を次のように分類できます。

$class (z) = sign (z' \hat{β} + \hat{b}) = sign (\hat{f} (z)) .$

$\hat{f} (z)$ は "分類スコア" であり、距離 z は決定境界からのものです。

数学的形成化: 双対問題- この方法は、双対二次計画法の問題を計算上より簡単に解くことができます。双対問題を取得するには、正のラグランジュ乗数 α_j に各制約を乗算し、目的関数から差し引きます。

$L_{P} = \frac{1}{2} β' β - \sum_{j} α_{j} (y_{j} (x_{j}' β + b) - 1),$

ここで、β および b に対する L_P の定常点を求めます。L_P の勾配を 0 に設定すると、以下を得られます。

\begin{matrix} β = \sum_{j} α_{j} y_{j} x_{j} \\ 0 = \sum_{j} α_{j} y_{j} . \end{matrix}

(1)

L_P に代入すると、以下のような双対 L_D を取得できます。

$L_{D} = \sum_{j} α_{j} - \frac{1}{2} \sum_{j} \sum_{k} α_{j} α_{k} y_{j} y_{k} x_{j}' x_{k},$

これを α_j ≥ 0 に対して最大化します。一般的に、多くの α_j は最大 0 です。双対問題での非ゼロの解 α_j は、式 1に示されているように、超平面を定義し、α_jy_jx_j の合計として β が返されます。非ゼロ α_j に対応するデータ点 x_j は "サポートベクター" です。

非ゼロ α_j についての L_D の微分は 0 が最適です。これにより、次が得られます。

$y_{j} f (x_{j}) - 1 = 0.$

特にこの場合は、非ゼロ α_j をもつ j を使用して、解で b の値が返されます。

双対は、標準の二次計画法の問題です。たとえば、Optimization Toolbox™ quadprog (Optimization Toolbox) ソルバーがこのタイプの問題を解きます。

不可分データ

データが分離超平面で使用できないことがあります。この場合は、SVM で "ソフトマージン"、つまり、すべてではありませんが、多くのデータ点を分離する超平面を使用できます。

以下の 2 つのソフトマージン標準定式化があります。どちらの場合も、スラック変数 ξ_j とペナルティパラメーター C の追加が必要です。

L¹ ノルム問題は次のとおりです。
$\min_{β, b, ξ} (\frac{1}{2} β' β + C \sum_{j} ξ_{j})$
条件
$\begin{matrix} y_{j} f (x_{j}) \geq 1 - ξ_{j} \\ ξ_{j} \geq 0. \end{matrix}$
L¹ ノルムは、ξ_j をその二乗ではなくスラック変数として使用することを表します。fitcsvm の 3 つのソルバーオプション (SMO、ISDA、L1QP) を使用すると、L¹ ノルム問題が最小化されます。
L² ノルム問題は次のとおりです。
$\min_{β, b, ξ} (\frac{1}{2} β' β + C \sum_{j} ξ_{j}^{2})$
これも、同じ制約の影響を受けます。

これらの定式化では、C を増やすとスラック変数 ξ_j に重みが追加されることがわかります。つまり、最適化によりクラス間の分離がより厳格になります。同様に、C を 0 まで下げると、誤分類の重要度が低下します。

数学的形成化: 双対問題- 計算を簡単にするために、このソフトマージン定式化に対する L¹ 双対問題を考えてみます。ラグランジュ乗数 μ_j を使用する場合、L¹ ノルム問題に対して最小化を行う関数は次のようになります。

$L_{P} = \frac{1}{2} β' β + C \sum_{j} ξ_{j} - \sum_{j} α_{j} (y_{i} f (x_{j}) - (1 - ξ_{j})) - \sum_{j} μ_{j} ξ_{j},$

ここで、β、b および正の ξ_j に対する L_P の定常点を求めます。L_P の勾配を 0 に設定すると、以下を得られます。

$\begin{matrix} β = \sum_{j} α_{j} y_{j} x_{j} \\ \sum_{j} α_{j} y_{j} = 0 \\ α_{j} = C - μ_{j} \\ α_{j}, μ_{j}, ξ_{j} \geq 0. \end{matrix}$

これらの方程式は、次の双対形式に直結します。

$\max_{α} \sum_{j} α_{j} - \frac{1}{2} \sum_{j} \sum_{k} α_{j} α_{k} y_{j} y_{k} x_{j}' x_{k}$

以下の制約に従います。

$\begin{array}{l} \sum_{j} y_{j} α_{j} = 0 \\ 0 \leq α_{j} \leq C . \end{array}$

最終的な不等式のセット 0 ≤ α_j ≤ C は、C が "ボックス制約" と呼ばれる理由を示します。C はラグランジュ乗数 α_j の許容値を "ボックス" の内部、つまり境界領域内に維持します。

b の勾配方程式は、サポートベクターに対応する、非ゼロ α_j のセットによって解 b を返します。

L² ノルム問題の双対は類似した方法で作成および解くことができます。詳細は、Christianini および Shawe-Taylor [2] (第 6 章) を参照してください。

fitcsvm の実装- 両方の双対ソフトマージン問題は、二次計画法の問題です。内部では、fitcsvm には問題を解くためのさまざまなアルゴリズムがあります。

1 クラスまたはバイナリ分類の場合、データで予想される外れ値の小数を設定していないと (OutlierFraction を参照)、既定のソルバーは逐次最小最適化 (SMO) になります。SMO を使用すると、一連の 2 点最小化により 1 ノルム問題が最小化されます。最適化時に、SMO では線形制約 $\sum_{i} α_{i} y_{i} = 0,$ が考慮され、モデルにバイアス項が明示的に含まれます。SMO は比較的高速です。SMO についての詳細は、[3]を参照してください。
バイナリ分類の場合、データで予想される外れ値の小数を設定すると、既定のソルバーは反復単一データアルゴリズム (ISDA) になります。SMO と同様、ISDA は 1 ノルム問題を解決します。SMO と異なり、ISDA は一連の 1 点最小化により最小化されますが、線形制約を考慮しません。また、モデルにはバイアス項が明示的に含まれません。ISDA についての詳細は、[4]を参照してください。
1 クラスまたはバイナリ分類の場合、Optimization Toolbox のライセンスがあると、quadprog (Optimization Toolbox) を使用して 1 ノルム問題を解くことができます。quadprog はメモリを大量に使用しますが、高い精度で二次計画法を解きます。詳細は、二次計画法の定義 (Optimization Toolbox)を参照してください。

カーネルによる非線形変換

一部のバイナリ分類問題には、有用な分離基準として使用できる単純な超平面がない場合もあります。これらの問題に対応して、SVM 分離超平面のほぼすべての単純化を保持する数学的アプローチのバリアントが存在します。

このアプローチでは、以下の再生カーネル理論の結果を使用します。

次のプロパティを使用する関数 G(x₁,x₂) のクラスがあります。線形空間 S と、次のように x を S にマッピングする関数 φ があります。
G(x₁,x₂) = <φ(x₁),φ(x₂)>.
ドット積は空間 S で生じます。
この関数のクラスには以下が含まれます。
- 多項式: 正の整数 p について、
  G(x₁,x₂) = (1 + x₁′x₂)^p.
- 放射基底関数 (ガウス):
  G(x₁,x₂) = exp(–∥x₁–x₂)∥²)
- 多層パーセプトロンまたはシグモイド (ニューラルネットワーク): 正の数値 p₁ および負の数値 p₂ について、
  G(x₁,x₂) = tanh(p₁x₁′x₂ + p₂)
  メモ
  p₁ および p₂ のすべてのセットが有効な再生カーネルを生成するわけではありません。
  fitcsvm ではシグモイドカーネルはサポートされません。代わりに、名前と値のペアの引数 'KernelFunction' を使用して、シグモイドカーネルを定義および指定できます。詳細については、カスタムカーネルを使用した SVM 分類器の学習を参照してください。

カーネルを使用する数学的アプローチは、超平面の計算方法に依存します。超平面分類のすべての計算では、ドット積しか使用しません。したがって、非線形カーネルでは同一の計算と解アルゴリズムを使用し、非線形の分類器を得ることができます。結果の分類器は、特定の空間 S の超平面になりますが、空間 S を識別または確認する必要はありません。

サポートベクターマシンの使用

教師あり学習モデルの場合と同様に、最初にサポートベクターマシンを学習し、その後分類器を交差検証します。この学習したマシンを使用して新しいデータを分類 (予測) します。また、満足する予測精度を得るには、さまざまな SVM カーネル関数を使用します。カーネル関数のパラメーターは調整しなければなりません。

SVM 分類器の学習
SVM 分類器による新しいデータの分類
SVM 分類器の調整

SVM 分類器の学習

fitcsvm を使用して、SVM 分類器を学習し、オプションで交差検証します。最も一般的な構文は、次のとおりです。

SVMModel = fitcsvm(X,Y,'KernelFunction','rbf',...
    'Standardize',true,'ClassNames',{'negClass','posClass'});

入力は以下になります。

X — 予測子データの行列。各行は 1 件の観測で、各行は 1 件の予測子です。
Y — 各行が X の該当行の値に対応しているクラスラベルの配列。Y は categorical 配列、文字配列、string 配列、logical ベクトル、数値ベクトル、または文字ベクトルの cell 配列です。
KernelFunction — 2 クラス学習の場合、既定値は 'linear' です。この値が適用されると、データは超平面で分割されます。値 'gaussian' (または 'rbf') は 1 クラス学習の場合の既定値であり、ガウス (または放射基底関数) カーネルを使用するよう指定します。SVM 分類器の学習を成功させるための重要な手順として、適切なカーネル関数を選択します。
Standardize — 分類器の学習の前に予測子を標準化する必要があるかどうかを示すフラグ。
ClassNames — 陰性クラスと陽性クラスを区別するか、データに含まれるクラスを指定します。陰性クラスは最初の要素 (または文字配列の行)、たとえば 'negClass' であり、陽性クラスは 2 番目の要素 (または文字配列の行)、たとえば 'posClass' です。ClassNames は Y と同じデータ型でなくてはなりません。クラス名を指定することをお勧めします。特に、異なる分類器の性能を比較する場合に有効です。

結果として、学習されたモデル (SVMModel) には、SVM アルゴリズムの最適化されたパラメーターが含まれ、新しいデータを分類できるようになります。

学習の管理に使用できる名前と値のペアについての詳細は、fitcsvm のリファレンスページを参照してください。

SVM 分類器による新しいデータの分類

predict を使用して新しいデータを分類します。学習された SVM 分類器 (SVMModel) を使用して新しいデータを分類するには、次の構文を使用します。

[label,score] = predict(SVMModel,newX);

結果として生成されるベクトル label は、X の各行の分類を表します。score はソフトスコアの n 行 2 列の行列です。各行は X の行 (新しい観測) に対応しています。最初の列には陰性のクラスに分類される観測のスコアが含まれ、2 番目の列には陽性のクラスに分類される観測のスコアが含まれます。

スコアではなく事後確率を推定する場合は、学習済み SVM 分類器 (SVMModel) を最初に fitPosterior に渡します。これは、スコアから事後確率への変換関数をスコアにあてはめます。構文は、次のようになります。

ScoreSVMModel = fitPosterior(SVMModel,X,Y);

分類器 ScoreSVMModel の ScoreTransform プロパティには、最適変換関数が含まれています。ScoreSVMModel を predict に渡します。出力引数 score はスコアは返しませんが、陰性のクラス (score の 1 列目) または陽性のクラス (score の 2 列目) に分類される観測の事後確率が格納されます。

SVM 分類器の調整

交差検証損失を最小化するパラメーター値を求めるには、fitcsvm の名前と値のペアの引数 'OptimizeHyperparameters' を使用します。使用可能なパラメーターは、'BoxConstraint'、'KernelFunction'、'KernelScale'、'PolynomialOrder' および 'Standardize' です。たとえば、ベイズ最適化の使用による SVM 分類器のあてはめの最適化を参照してください。または、bayesopt を使用した交差検証済み SVM 分類器の最適化で示されているように、関数 bayesopt を使用できます。関数 bayesopt では、より柔軟に最適化をカスタマイズできます。関数 bayesopt は、関数 fitcsvm では最適化できないパラメーターを含む任意のパラメーターの最適化に使用できます。

以下の方法に従って、分類器のパラメーターを手動で調整することもできます。

データを fitcsvm に渡し、名前と値のペアの引数 'KernelScale','auto' を設定します。学習された SVM モデルの名前が SVMModel であるとします。ヒューリスティック手法に従ってカーネルスケールが選択されます。ヒューリスティック手法では副標本抽出を使用します。そのため、結果を再現するには、分類器の学習の前に rng を使用して乱数シードを設定します。
分類器を crossval に渡して、分類器を交差検証します。既定では、10 分割交差検証が実行されます。
交差検証された SVM モデルを kfoldLoss に渡し、分類誤差を推定して保持します。
SVM 分類器に再学習させますが、'KernelScale' および 'BoxConstraint' 名前と値のペアの引数は調整します。
- BoxConstraint — 方法の 1 つに、ボックス制約パラメーターの等比数列を試すというものがあります。たとえば、1e-5 から 1e5 まで 10 倍ずつ大きくした 11 個の値があるとします。BoxConstraint を大きくすると、サポートベクターの数が少なくなる可能性がありますが、学習時間が長くなる可能性もあります。
- KernelScale — 方法の 1 つに、元のカーネルスケールでスケーリングされた RBF sigma パラメーターの等比数列を試すというものがあります。以下のように行います。
  1. 元のカーネルスケール (たとえば ks など) をドット表記を使用して取得します。ks = SVMModel.KernelParameters.Scale
  2. 元の係数を新しいカーネルスケール係数として使用します。たとえば、ks を 1e-5 ～ 1e5 の 11 個の値 (係数 10 ずつ増加) で乗算します。

分類誤差が最小になるモデルを選択します。精度を高めるためにパラメーターの調整がさらに必要となる場合があります。初期パラメーターから開始し、今回は係数 1.2 を使用して、別の交差検証ステップを実行します。

ガウスカーネルを使用した SVM 分類器の学習

ライブスクリプトを開く

この例では、ガウスカーネル関数を使用する非線形分類器の生成方法を示します。まず、2 次元の単位円板内の点から構成されるクラスと、半径 1 から半径 2 までの環形内の点から構成される別のクラスを生成します。次に、ガウス放射基底関数カーネルを使用したデータに基づいて分類器を生成します。モデルは円対称であるため、既定の線形分類器は明らかにこの問題には不適切です。ボックス制約パラメーターを Inf に設定し、厳格な分類にします。つまり、誤分類された学習点がなくなります。他のカーネル関数は、厳格な分類を提供できない場合があるため、この厳格なボックス制約では機能しない可能性があります。rbf 分類器でクラスを分離できる場合でも、結果が過学習になる可能性があります。

単位円板に一様分布する 100 個の点を生成します。これを行うため、一様確率変数の平方根により半径 r を生成し、(0, $2 π$ ) 内で一様に角度 t を生成して、(r cos( t ), r sin( t )) に点を配置します。

rng(1); % For reproducibility
r = sqrt(rand(100,1)); % Radius
t = 2*pi*rand(100,1);  % Angle
data1 = [r.*cos(t), r.*sin(t)]; % Points

環帯に一様分布する 100 個の点を生成します。この場合も半径は平方根に比例しますが、今回は 1 ～ 4 の一様分布の平方根になります。

r2 = sqrt(3*rand(100,1)+1); % Radius
t2 = 2*pi*rand(100,1);      % Angle
data2 = [r2.*cos(t2), r2.*sin(t2)]; % points

点をプロットし、比較のため半径 1 と 2 の円をプロットします。

figure;
plot(data1(:,1),data1(:,2),'r.','MarkerSize',15)
hold on
plot(data2(:,1),data2(:,2),'b.','MarkerSize',15)
ezpolar(@(x)1);ezpolar(@(x)2);
axis equal
hold off

データを 1 つの行列に挿入し、分類のベクトルを作成します。

data3 = [data1;data2];
theclass = ones(200,1);
theclass(1:100) = -1;

KernelFunction を 'rbf' に、BoxConstraint を Inf に設定して SVM 分類器に学習させます。判定境界をプロットして、サポートベクトルのフラグを設定します。

%Train the SVM Classifier
cl = fitcsvm(data3,theclass,'KernelFunction','rbf',...
    'BoxConstraint',Inf,'ClassNames',[-1,1]);

% Predict scores over the grid
d = 0.02;
[x1Grid,x2Grid] = meshgrid(min(data3(:,1)):d:max(data3(:,1)),...
    min(data3(:,2)):d:max(data3(:,2)));
xGrid = [x1Grid(:),x2Grid(:)];
[~,scores] = predict(cl,xGrid);

% Plot the data and the decision boundary
figure;
h(1:2) = gscatter(data3(:,1),data3(:,2),theclass,'rb','.');
hold on
ezpolar(@(x)1);
h(3) = plot(data3(cl.IsSupportVector,1),data3(cl.IsSupportVector,2),'ko');
contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'k');
legend(h,{'-1','+1','Support Vectors'});
axis equal
hold off

fitcsvm によって半径 1 の円に近い分類器が生成されます。学習データが無作為であるために差異が発生しています。

既定のパラメーターで学習させると、さらに円状に近い分類境界になりますが、一部の学習データが誤分類されます。また、BoxConstraint の既定値は 1 です。そのため、サポートベクトルの数はさらに多くなります。

cl2 = fitcsvm(data3,theclass,'KernelFunction','rbf');
[~,scores2] = predict(cl2,xGrid);

figure;
h(1:2) = gscatter(data3(:,1),data3(:,2),theclass,'rb','.');
hold on
ezpolar(@(x)1);
h(3) = plot(data3(cl2.IsSupportVector,1),data3(cl2.IsSupportVector,2),'ko');
contour(x1Grid,x2Grid,reshape(scores2(:,2),size(x1Grid)),[0 0],'k');
legend(h,{'-1','+1','Support Vectors'});
axis equal
hold off

カスタムカーネルを使用した SVM 分類器の学習

スクリプトを開く

この例では、シグモイドカーネルなどのカスタムカーネル関数を使用して SVM 分類器を学習させる方法と、カスタムカーネル関数のパラメーターを調整する方法を説明します。

単位円の内部に無作為な点集合を生成します。第 1 象限と第 3 象限の点を陽性のクラスとしてラベル付けし、第 2 象限と第 4 象限の点を陰性のクラスとしてラベル付けします。

rng(1);  % For reproducibility
n = 100; % Number of points per quadrant

r1 = sqrt(rand(2*n,1));                     % Random radii
t1 = [pi/2*rand(n,1); (pi/2*rand(n,1)+pi)]; % Random angles for Q1 and Q3
X1 = [r1.*cos(t1) r1.*sin(t1)];             % Polar-to-Cartesian conversion

r2 = sqrt(rand(2*n,1));
t2 = [pi/2*rand(n,1)+pi/2; (pi/2*rand(n,1)-pi/2)]; % Random angles for Q2 and Q4
X2 = [r2.*cos(t2) r2.*sin(t2)];

X = [X1; X2];        % Predictors
Y = ones(4*n,1);
Y(2*n + 1:end) = -1; % Labels

データをプロットする。

figure;
gscatter(X(:,1),X(:,2),Y);
title('Scatter Diagram of Simulated Data')

特徴空間の 2 つの行列を入力として受け入れ、シグモイドカーネルを使用してグラム行列に変換する関数を作成します。

function G = mysigmoid(U,V)
% Sigmoid kernel function with slope gamma and intercept c
gamma = 1;
c = -1;
G = tanh(gamma*U*V' + c);
end

このコードを mysigmoid という名前のファイルとして MATLAB® パス上に保存します。

シグモイドカーネル関数を使用して SVM 分類器に学習させます。データの標準化をお勧めします。

Mdl1 = fitcsvm(X,Y,'KernelFunction','mysigmoid','Standardize',true);

Mdl1 は、推定されたパラメーターを含む ClassificationSVM 分類器です。

データをプロットし、サポートベクトルと判定境界を特定します。

    % Compute the scores over a grid
d = 0.02; % Step size of the grid
[x1Grid,x2Grid] = meshgrid(min(X(:,1)):d:max(X(:,1)),...
    min(X(:,2)):d:max(X(:,2)));
xGrid = [x1Grid(:),x2Grid(:)];        % The grid
[~,scores1] = predict(Mdl1,xGrid); % The scores

figure;
h(1:2) = gscatter(X(:,1),X(:,2),Y);
hold on
h(3) = plot(X(Mdl1.IsSupportVector,1),...
    X(Mdl1.IsSupportVector,2),'ko','MarkerSize',10);
    % Support vectors
contour(x1Grid,x2Grid,reshape(scores1(:,2),size(x1Grid)),[0 0],'k');
    % Decision boundary
title('Scatter Diagram with the Decision Boundary')
legend({'-1','1','Support Vectors'},'Location','Best');
hold off

カーネルのパラメーターを調整すると、判定境界の形状を改善することができます。これによりサンプル内の誤分類率も低下する場合もありますが、最初にサンプル外の誤分類率を特定する必要があります。

10 分割交差検証を使用して、標本外の誤分類率を特定します。

CVMdl1 = crossval(Mdl1);
misclass1 = kfoldLoss(CVMdl1);
misclass1

misclass1 =

    0.1350

サンプル外の誤分類率は 13.5% です。

別のシグモイド関数を作成します。ただし、gamma = 0.5; に設定します。

function G = mysigmoid2(U,V)
% Sigmoid kernel function with slope gamma and intercept c
gamma = 0.5;
c = -1;
G = tanh(gamma*U*V' + c);
end

このコードを mysigmoid2 という名前のファイルとして MATLAB® パス上に保存します。

調整したシグモイドカーネルを使用して、別の SVM 分類器に学習をさせます。データと決定領域をプロットし、サンプル外の誤分類率を特定します。

Mdl2 = fitcsvm(X,Y,'KernelFunction','mysigmoid2','Standardize',true);
[~,scores2] = predict(Mdl2,xGrid);

figure;
h(1:2) = gscatter(X(:,1),X(:,2),Y);
hold on
h(3) = plot(X(Mdl2.IsSupportVector,1),...
    X(Mdl2.IsSupportVector,2),'ko','MarkerSize',10);
title('Scatter Diagram with the Decision Boundary')
contour(x1Grid,x2Grid,reshape(scores2(:,2),size(x1Grid)),[0 0],'k');
legend({'-1','1','Support Vectors'},'Location','Best');
hold off

CVMdl2 = crossval(Mdl2);
misclass2 = kfoldLoss(CVMdl2);
misclass2

misclass2 =

    0.0450

シグモイド勾配調整を行うと、新しい判定境界によりサンプル内近似の精度は高くなりなますが、交差検証率は 66% 以上縮小します。

ベイズ最適化の使用による SVM 分類器のあてはめの最適化

ライブスクリプトを開く

この例では、関数 fitcsvm および名前と値のペア OptimizeHyperparameters を使用して SVM 分類器を最適化する方法を示します。この分類は、混合ガウスモデルによる点の位置に作用します。モデルの詳細については、The Elements of Statistical Learning, Hastie, Tibshirani, and Friedman (2009) の 17 ページを参照してください。このモデルでは、平均 (1,0) および単位分散をもつ 2 次元の独立した正規分布になっている 10 個の基底点をはじめに "green" クラスについて生成します。また、平均 (0,1) と単位分散による 2 次元の独立した正規として分布される "red" クラスにも、10 個の基底点が生成されます。クラス (green と red) ごとに、次のように 100 個の無作為な点を生成します。

適切な色の基底点 m を一様にランダムに選択します。
平均 m と分散 I/5 (I は 2 行 2 列の単位行列) をもつ 2 次元正規分布を使用して、独立した無作為な点を生成します。最適化のアドバンテージをより明確に示すため、この例では I/50 という分散を使用します。

点と分類器の生成

クラスごとに 10 個の基底点を生成します。

rng default % For reproducibility
grnpop = mvnrnd([1,0],eye(2),10);
redpop = mvnrnd([0,1],eye(2),10);

基底点を表示します。

plot(grnpop(:,1),grnpop(:,2),'go')
hold on
plot(redpop(:,1),redpop(:,2),'ro')
hold off

赤の基底点の一部が緑の基底点の近くにあるため、位置のみによるデータ点の分類は難しいかもしれません。

各クラスについて 100 個ずつのデータ点を生成します。

redpts = zeros(100,2);grnpts = redpts;
for i = 1:100
    grnpts(i,:) = mvnrnd(grnpop(randi(10),:),eye(2)*0.02);
    redpts(i,:) = mvnrnd(redpop(randi(10),:),eye(2)*0.02);
end

データ点を表示します。

figure
plot(grnpts(:,1),grnpts(:,2),'go')
hold on
plot(redpts(:,1),redpts(:,2),'ro')
hold off

分類用のデータの準備

データを 1 つの行列に格納し、各点のクラスにラベルを付けるベクトル grp を作成します。

cdata = [grnpts;redpts];
grp = ones(200,1);
% Green label 1, red label -1
grp(101:200) = -1;

交差検証の準備

交差検証用の分割を設定します。これにより、各ステップで最適化に使用される学習セットと検定セットが決まります。

c = cvpartition(200,'KFold',10);

近似の最適化

適切な近似、つまり交差検証損失が小さい近似を求めるため、ベイズ最適化を使用するようにオプションを設定します。すべての最適化で同じ交差検証分割 c を使用します。

再現性を得るために、'expected-improvement-plus' の獲得関数を使用します。

opts = struct('Optimizer','bayesopt','ShowPlots',true,'CVPartition',c,...
    'AcquisitionFunctionName','expected-improvement-plus');
svmmod = fitcsvm(cdata,grp,'KernelFunction','rbf',...
    'OptimizeHyperparameters','auto','HyperparameterOptimizationOptions',opts)

|=====================================================================================================|
| Iter | Eval   | Objective   | Objective   | BestSoFar   | BestSoFar   | BoxConstraint|  KernelScale |
|      | result |             | runtime     | (observed)  | (estim.)    |              |              |
|=====================================================================================================|
|    1 | Best   |       0.345 |     0.35552 |       0.345 |       0.345 |      0.00474 |       306.44 |

|    2 | Best   |       0.115 |     0.26672 |       0.115 |     0.12678 |       430.31 |       1.4864 |

|    3 | Accept |        0.52 |     0.21117 |       0.115 |      0.1152 |     0.028415 |     0.014369 |

|    4 | Accept |        0.61 |     0.44425 |       0.115 |     0.11504 |       133.94 |    0.0031427 |

|    5 | Accept |        0.34 |     0.46272 |       0.115 |     0.11504 |     0.010993 |       5.7742 |

|    6 | Best   |       0.085 |     0.20733 |       0.085 |    0.085039 |       885.63 |      0.68403 |

|    7 | Accept |       0.105 |     0.18607 |       0.085 |    0.085428 |       0.3057 |      0.58118 |

|    8 | Accept |        0.21 |     0.22047 |       0.085 |     0.09566 |      0.16044 |      0.91824 |

|    9 | Accept |       0.085 |     0.20821 |       0.085 |     0.08725 |       972.19 |      0.46259 |

|   10 | Accept |         0.1 |     0.36196 |       0.085 |    0.090952 |       990.29 |        0.491 |

|   11 | Best   |        0.08 |     0.34245 |        0.08 |    0.079362 |       2.5195 |        0.291 |

|   12 | Accept |        0.09 |     0.31556 |        0.08 |     0.08402 |       14.338 |      0.44386 |

|   13 | Accept |         0.1 |     0.20414 |        0.08 |     0.08508 |    0.0022577 |      0.23803 |

|   14 | Accept |        0.11 |     0.26065 |        0.08 |    0.087378 |       0.2115 |      0.32109 |

|   15 | Best   |        0.07 |     0.21537 |        0.07 |    0.081507 |        910.2 |      0.25218 |

|   16 | Best   |       0.065 |     0.21929 |       0.065 |    0.072457 |       953.22 |      0.26253 |

|   17 | Accept |       0.075 |     0.30938 |       0.065 |    0.072554 |       998.74 |      0.23087 |

|   18 | Accept |       0.295 |     0.23381 |       0.065 |    0.072647 |       996.18 |       44.626 |

|   19 | Accept |        0.07 |     0.45571 |       0.065 |     0.06946 |       985.37 |      0.27389 |

|   20 | Accept |       0.165 |     0.30732 |       0.065 |    0.071622 |     0.065103 |      0.13679 |

|=====================================================================================================|
| Iter | Eval   | Objective   | Objective   | BestSoFar   | BestSoFar   | BoxConstraint|  KernelScale |
|      | result |             | runtime     | (observed)  | (estim.)    |              |              |
|=====================================================================================================|
|   21 | Accept |       0.345 |     0.22239 |       0.065 |    0.071764 |        971.7 |       999.01 |

|   22 | Accept |        0.61 |     0.21176 |       0.065 |    0.071967 |    0.0010168 |    0.0010005 |

|   23 | Accept |       0.345 |     0.30206 |       0.065 |    0.071959 |    0.0010674 |       999.18 |

|   24 | Accept |        0.35 |     0.20608 |       0.065 |    0.071863 |    0.0010003 |       40.628 |

|   25 | Accept |        0.24 |     0.43093 |       0.065 |    0.072124 |       996.55 |       10.423 |

|   26 | Accept |        0.61 |     0.45597 |       0.065 |    0.072068 |       958.64 |    0.0010026 |

|   27 | Accept |        0.47 |     0.21528 |       0.065 |     0.07218 |       993.69 |     0.029723 |

|   28 | Accept |         0.3 |     0.19438 |       0.065 |    0.072291 |       993.15 |       170.01 |

|   29 | Accept |        0.16 |     0.41633 |       0.065 |    0.072104 |       992.81 |       3.8594 |

|   30 | Accept |       0.365 |     0.20951 |       0.065 |    0.072112 |    0.0010017 |     0.044287 |

__________________________________________________________
Optimization completed.
MaxObjectiveEvaluations of 30 reached.
Total function evaluations: 30
Total elapsed time: 84.1762 seconds
Total objective function evaluation time: 8.6528

Best observed feasible point:
    BoxConstraint    KernelScale
    _____________    ___________

       953.22          0.26253  

Observed objective function value = 0.065
Estimated objective function value = 0.073726
Function evaluation time = 0.21929

Best estimated feasible point (according to models):
    BoxConstraint    KernelScale
    _____________    ___________

       985.37          0.27389  

Estimated objective function value = 0.072112
Estimated function evaluation time = 0.27508

svmmod = 
  ClassificationSVM
                         ResponseName: 'Y'
                CategoricalPredictors: []
                           ClassNames: [-1 1]
                       ScoreTransform: 'none'
                      NumObservations: 200
    HyperparameterOptimizationResults: [1x1 BayesianOptimization]
                                Alpha: [77x1 double]
                                 Bias: -0.2352
                     KernelParameters: [1x1 struct]
                       BoxConstraints: [200x1 double]
                      ConvergenceInfo: [1x1 struct]
                      IsSupportVector: [200x1 logical]
                               Solver: 'SMO'


  Properties, Methods

最適化したモデルの損失を求めます。

lossnew = kfoldLoss(fitcsvm(cdata,grp,'CVPartition',c,'KernelFunction','rbf',...
    'BoxConstraint',svmmod.HyperparameterOptimizationResults.XAtMinObjective.BoxConstraint,...
    'KernelScale',svmmod.HyperparameterOptimizationResults.XAtMinObjective.KernelScale))

lossnew = 0.0650

この損失は、[観測された目的関数値] 下の最適化の出力で表示される損失と同じです。

最適化された分類器を可視化します。

d = 0.02;
[x1Grid,x2Grid] = meshgrid(min(cdata(:,1)):d:max(cdata(:,1)),...
    min(cdata(:,2)):d:max(cdata(:,2)));
xGrid = [x1Grid(:),x2Grid(:)];
[~,scores] = predict(svmmod,xGrid);
figure;
h = nan(3,1); % Preallocation
h(1:2) = gscatter(cdata(:,1),cdata(:,2),grp,'rg','+*');
hold on
h(3) = plot(cdata(svmmod.IsSupportVector,1),...
    cdata(svmmod.IsSupportVector,2),'ko');
contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'k');
legend(h,{'-1','+1','Support Vectors'},'Location','Southeast');
axis equal
hold off

SVM 分類モデルの事後確率領域のプロット

ライブスクリプトを開く

この例では、観測値のグリッドで SVM モデルの事後確率を予測して事後確率をグリッドにプロットする方法を示します。事後確率をプロットすると、判定境界が明らかになります。

フィッシャーのアヤメのデータセットを読み込みます。花弁の長さと幅を使用して分類器に学習させ、データから virginica 種を削除します。

load fisheriris
classKeep = ~strcmp(species,'virginica');
X = meas(classKeep,3:4);
y = species(classKeep);

データを使用して SVM 分類器に学習させます。クラスの順序を指定することをお勧めします。

SVMModel = fitcsvm(X,y,'ClassNames',{'setosa','versicolor'});

最適なスコア変換関数を推定します。

rng(1); % For reproducibility
[SVMModel,ScoreParameters] = fitPosterior(SVMModel);

Warning: Classes are perfectly separated. The optimal score-to-posterior transformation is a step function.

ScoreParameters

ScoreParameters = struct with fields:
                        Type: 'step'
                  LowerBound: -0.8431
                  UpperBound: 0.6897
    PositiveClassProbability: 0.5000

クラスは可分であるため、最適なスコア変換関数はステップ関数です。ScoreParameters の LowerBound と UpperBound の各フィールドは、クラス分離超平面 (マージン) 内の観測値に対応するスコアの範囲の上限と下限を表します。学習の観測値はマージン内には入りません。新しいスコアがこの範囲内にある場合、対応する観測値に陽性のクラスの事後確率 (ScoreParameters の PositiveClassProbability フィールドの値) が割り当てられます。

観測された予測子領域の値のグリッドを定義します。グリッド内の各インスタンスの事後確率を予測します。

xMax = max(X);
xMin = min(X);
d = 0.01;
[x1Grid,x2Grid] = meshgrid(xMin(1):d:xMax(1),xMin(2):d:xMax(2));

[~,PosteriorRegion] = predict(SVMModel,[x1Grid(:),x2Grid(:)]);

学習データの陽性のクラスの事後確率領域をプロットします。

figure;
contourf(x1Grid,x2Grid,...
        reshape(PosteriorRegion(:,2),size(x1Grid,1),size(x1Grid,2)));
h = colorbar;
h.Label.String = 'P({\it{versicolor}})';
h.YLabel.FontSize = 16;
caxis([0 1]);
colormap jet;

hold on
gscatter(X(:,1),X(:,2),y,'mc','.x',[15,10]);
sv = X(SVMModel.IsSupportVector,:);
plot(sv(:,1),sv(:,2),'yo','MarkerSize',15,'LineWidth',2);
axis tight
hold off

2 クラス学習では、クラスが可分な場合、陽性クラスの事後確率 0 を観測値に割り当てる領域、1 を割り当てる領域、および陽性クラスの事前確率を割り当てる領域という 3 つの領域があります。

線形サポートベクターマシンによるイメージの分析

ライブスクリプトを開く

この例では、線形 SVM バイナリ学習器から構成される誤り訂正出力符号 (ECOC) モデルを学習させることにより、イメージのどの象限に形状があるかを判別する方法を示します。また、サポートベクター、そのラベル、および推定した $α$ 係数を格納する ECOC モデルのディスク領域の消費状況も示します。

データセットの作成

50 x 50 のイメージに半径が 5 の円を無作為に配置します。5000 個のイメージを作成します。円がある象限を示すラベルを各イメージについて作成します。第 1 象限は右上、第 2 象限は左上、第 3 象限は左下、第 4 象限は右下です。予測子は各ピクセルの明度です。

d = 50;  % Height and width of the images in pixels
n = 5e4; % Sample size

X = zeros(n,d^2); % Predictor matrix preallocation 
Y = zeros(n,1);   % Label preallocation
theta = 0:(1/d):(2*pi);
r = 5;            % Circle radius
rng(1);           % For reproducibility

for j = 1:n
    figmat = zeros(d);                       % Empty image
    c = datasample((r + 1):(d - r - 1),2);   % Random circle center
    x = r*cos(theta) + c(1);                 % Make the circle 
    y = r*sin(theta) + c(2);               
    idx = sub2ind([d d],round(y),round(x));  % Convert to linear indexing
    figmat(idx) = 1;                         % Draw the circle
    X(j,:) = figmat(:);                % Store the data
    Y(j) = (c(2) >= floor(d/2)) + 2*(c(2) < floor(d/2)) + ...
        (c(1) < floor(d/2)) + ...
        2*((c(1) >= floor(d/2)) & (c(2) < floor(d/2))); % Determine the quadrant
end

観測値をプロットします。

figure
imagesc(figmat)
h = gca;
h.YDir = 'normal';
title(sprintf('Quadrant %d',Y(end)))

ECOC モデルを学習させます。

25% のホールドアウト標本を使用し、学習標本とホールドアウト標本のインデックスを指定します。

p = 0.25;
CVP = cvpartition(Y,'Holdout',p); % Cross-validation data partition
isIdx = training(CVP);            % Training sample indices
oosIdx = test(CVP);               % Test sample indices

バイナリ学習器のサポートベクターを格納するように指定する SVM テンプレートを作成します。このテンプレートと学習データを fitcecoc に渡してモデルを学習させます。学習標本の分類誤差を決定します。

t = templateSVM('SaveSupportVectors',true);
MdlSV = fitcecoc(X(isIdx,:),Y(isIdx),'Learners',t);
isLoss = resubLoss(MdlSV)

isLoss = 0

MdlSV は、学習させたマルチクラスモデル ClassificationECOC です。これには、各バイナリ学習器の学習データとサポートベクターが格納されています。イメージ分析のようにデータセットが大きい場合、このモデルは大量にメモリを消費する可能性があります。

この ECOC モデルによるディスク領域の消費量を調べます。

infoMdlSV = whos('MdlSV');
mbMdlSV = infoMdlSV.bytes/1.049e6

mbMdlSV = 763.6151

このモデルでは 763.6 MB を消費しています。

モデルの効率の向上

標本外性能を評価できます。また、サポートベクター、関連パラメーターおよび学習データが含まれていない圧縮したモデルで過適合になっているかどうかも評価できます。

学習させた ECOC モデルからサポートベクターおよび関連パラメーターを破棄します。次に、compact を使用して、結果のモデルから学習データを破棄します。

Mdl = discardSupportVectors(MdlSV);
CMdl = compact(Mdl);
info = whos('Mdl','CMdl');
[bytesCMdl,bytesMdl] = info.bytes;
memReduction = 1 - [bytesMdl bytesCMdl]/infoMdlSV.bytes

memReduction = 1×2

    0.0626    0.9996

この場合、サポートベクターを破棄するとメモリ消費量が約 6% 削減されます。サポートベクターを圧縮および破棄すると、サイズが約 99.96% 小さくなります。

サポートベクターの代替の管理方法として、大きいボックス制約 (100 など) を指定して学習中のサポートベクター数を削減します。使用するサポートベクター数が少ない SVM モデルが望ましく、消費メモリも少なくなりますが、ボックス制約の値を大きくすると学習時間が長くなる傾向があります。

MdlSV と Mdl をワークスペースから削除します。

clear Mdl MdlSV

ホールドアウト標本性能の評価

ホールドアウト標本の分類誤差を計算します。ホールドアウト標本の予測のサンプルをプロットします。

oosLoss = loss(CMdl,X(oosIdx,:),Y(oosIdx))

oosLoss = 0

yHat = predict(CMdl,X(oosIdx,:));
nVec = 1:size(X,1);
oosIdx = nVec(oosIdx);

figure;
for j = 1:9
    subplot(3,3,j)
    imagesc(reshape(X(oosIdx(j),:),[d d]))
    h = gca;
    h.YDir = 'normal';
    title(sprintf('Quadrant: %d',yHat(j)))
end
text(-1.33*d,4.5*d + 1,'Predictions','FontSize',17)

このモデルでは、すべてのホールドアウト標本の観測値が正しく分類されています。

参考

bayesopt | fitcsvm | kfoldLoss

参照

[1] Hastie, T., R. Tibshirani, and J. Friedman. The Elements of Statistical Learning, second edition. New York: Springer, 2008.

[2] Christianini, N., and J. Shawe-Taylor. An Introduction to Support Vector Machines and Other Kernel-Based Learning Methods. Cambridge, UK: Cambridge University Press, 2000.

[3] Fan, R.-E., P.-H. Chen, and C.-J. Lin. “Working set selection using second order information for training support vector machines.” Journal of Machine Learning Research, Vol 6, 2005, pp. 1889–1918.

[4] Kecman V., T. -M. Huang, and M. Vogt. “Iterative Single Data Algorithm for Training Kernel Machines from Huge Data Sets: Theory and Performance.” In Support Vector Machines: Theory and Applications. Edited by Lipo Wang, 255–274. Berlin: Springer-Verlag, 2005.