Taiyo Mineo
Easy Support Vector Machine
Diff: SVM.cpp
- Revision:
- 5:792afbb0bcf3
- Parent:
- 4:01a20b89db32
--- a/SVM.cpp Mon Feb 16 07:53:14 2015 +0000
+++ b/SVM.cpp Wed Feb 18 15:01:12 2015 +0000
@@ -2,252 +2,255 @@
SVM::SVM(int dim_sample, int n_sample, float* sample_data, int* sample_label)
{
- this->dim_signal = dim_sample;
- this->n_sample = n_sample;
+ this->dim_signal = dim_sample;
+ this->n_sample = n_sample;
- // 各配列(ベクトル)のアロケート
- alpha = new float[n_sample];
- // grammat = new float[n_sample * n_sample];
- label = new int[n_sample];
- sample = new float[dim_signal * n_sample];
- sample_max = new float[dim_signal];
- sample_min = new float[dim_signal];
+ // 各配列(ベクトル)のアロケート
+ alpha = new float[n_sample];
+ // grammat = new float[n_sample * n_sample];
+ label = new int[n_sample];
+ sample = new float[dim_signal * n_sample];
+ sample_max = new float[dim_signal];
+ sample_min = new float[dim_signal];
+
+ // サンプルのコピー
+ memcpy(this->sample, sample_data,
+ sizeof(float) * dim_signal * n_sample);
+ memcpy(this->label, sample_label,
+ sizeof(int) * n_sample);
- // サンプルのコピー
- memcpy(this->sample, sample_data,
- sizeof(float) * dim_signal * n_sample);
- memcpy(this->label, sample_label,
- sizeof(int) * n_sample);
+ // 正規化のための最大最小値
+ memset(sample_max, 0, sizeof(float) * dim_sample);
+ memset(sample_min, 0, sizeof(float) * dim_sample);
+ for (int i = 0; i < dim_signal; i++) {
+ float value;
+ sample_min[i] = FLT_MAX;
+ for (int j = 0; j < n_sample; j++) {
+ value = MATRIX_AT(this->sample, dim_signal, j, i);
+ if ( value > sample_max[i] ) {
+ sample_max[i] = value;
+ } else if ( value < sample_min[i] ) {
+ //printf("min[%d] : %f -> ", i, sample_min[i]);
+ sample_min[i] = value;
+ //printf("min[%d] : %f \dim_signal", i, value);
+ }
+ }
+ }
- // 正規化のための最大最小値
- memset(sample_max, 0, sizeof(float) * dim_sample);
- memset(sample_min, 0, sizeof(float) * dim_sample);
- for (int i = 0; i < dim_signal; i++) {
- float value;
- sample_min[i] = FLT_MAX;
- for (int j = 0; j < n_sample; j++) {
- value = MATRIX_AT(this->sample, dim_signal, j, i);
- if ( value > sample_max[i] ) {
- sample_max[i] = value;
- } else if ( value < sample_min[i] ) {
- //printf("min[%d] : %f -> ", i, sample_min[i]);
- sample_min[i] = value;
- //printf("min[%d] : %f \dim_signal", i, value);
+ // 信号の正規化 : 死ぬほど大事
+ for (int i = 0; i < dim_signal; i++) {
+ float max,min;
+ max = sample_max[i];
+ min = sample_min[i];
+ for (int j = 0; j < n_sample; j++) {
+ // printf("[%d,%d] %f -> ", i, j, MATRIX_AT(this->sample, dim_signal, j, i));
+ MATRIX_AT(this->sample, dim_signal, j, i) = ( MATRIX_AT(this->sample, dim_signal, j, i) - min ) / (max - min);
+ // printf("%f \r\n", MATRIX_AT(this->sample, dim_signal, j, i));
+ }
+ }
+
+ /* // グラム行列の計算 : メモリの制約上,廃止
+ for (int i = 0; i < n_sample; i++) {
+ for (int j = i; j < n_sample; j++) {
+ MATRIX_AT(grammat,n_sample,i,j) = kernel_function(&(MATRIX_AT(this->sample,dim_signal,i,0)), &(MATRIX_AT(this->sample,dim_signal,j,0)), dim_signal);
+ // グラム行列は対称
+ if ( i != j ) {
+ MATRIX_AT(grammat,n_sample,j,i) = MATRIX_AT(grammat,n_sample,i,j);
+ }
}
}
- }
-
- // 信号の正規化 : 死ぬほど大事
- for (int i = 0; i < dim_signal; i++) {
- float max,min;
- max = sample_max[i];
- min = sample_min[i];
- for (int j = 0; j < n_sample; j++) {
- // printf("[%d,%d] %f -> ", i, j, MATRIX_AT(this->sample, dim_signal, j, i));
- MATRIX_AT(this->sample, dim_signal, j, i) = ( MATRIX_AT(this->sample, dim_signal, j, i) - min ) / (max - min);
- // printf("%f \r\n", MATRIX_AT(this->sample, dim_signal, j, i));
- }
- }
+ */
- /* // グラム行列の計算 : メモリの制約上,廃止
- for (int i = 0; i < n_sample; i++) {
- for (int j = i; j < n_sample; j++) {
- MATRIX_AT(grammat,n_sample,i,j) = kernel_function(&(MATRIX_AT(this->sample,dim_signal,i,0)), &(MATRIX_AT(this->sample,dim_signal,j,0)), dim_signal);
- // グラム行列は対称
- if ( i != j ) {
- MATRIX_AT(grammat,n_sample,j,i) = MATRIX_AT(grammat,n_sample,i,j);
- }
- }
- }
- */
+ // 学習関連の設定. 例によって経験則
+ this->maxIteration = 5000;
+ this->epsilon = float(0.00001);
+ this->eta = float(0.05);
+ this->learn_alpha = float(0.8) * this->eta;
+ this->status = SVM_NOT_LEARN;
- // 学習関連の設定. 例によって経験則
- this->maxIteration = 5000;
- this->epsilon = float(0.00001);
- this->eta = float(0.05);
- this->learn_alpha = float(0.8) * this->eta;
- this->status = SVM_NOT_LEARN;
+ // ソフトマージンの係数. 両方ともFLT_MAXとすることでハードマージンと一致.
+ // また, 設定するときはどちらか一方のみにすること.
+ C1 = FLT_MAX;
+ C2 = 5;
- // ソフトマージンの係数. 両方ともFLT_MAXとすることでハードマージンと一致.
- // また, 設定するときはどちらか一方のみにすること.
- C1 = FLT_MAX;
- C2 = 5;
-
- srand((unsigned int)time(NULL));
-
+ srand((unsigned int)time(NULL));
}
// 楽園追放
-SVM::~SVM(void)
+SVM::~SVM(void)
{
- delete [] alpha; delete [] label;
+ delete [] alpha;
+ delete [] label;
delete [] sample;
- delete [] sample_max; delete [] sample_min;
+ delete [] sample_max;
+ delete [] sample_min;
}
// 再急勾配法（ｻｰｾﾝwww）による学習
int SVM::learning(void)
{
- int iteration; // 学習繰り返しカウント
-
- float* diff_alpha; // 双対問題の勾配値
- float* pre_diff_alpha; // 双対問題の前回の勾配値（慣性項に用いる）
- float* pre_alpha; // 前回の双対係数ベクトル（収束判定に用いる）
- register float diff_sum; // 勾配計算用の小計
- register float kernel_val; // カーネル関数とC2を含めた項
-
- //float plus_sum, minus_sum; // 正例と負例の係数和
+ int iteration; // 学習繰り返しカウント
- // 配列（ベクトル）のアロケート
- diff_alpha = new float[n_sample];
- pre_diff_alpha = new float[n_sample];
- pre_alpha = new float[n_sample];
-
- status = SVM_NOT_LEARN; // 学習を未完了に
- iteration = 0; // 繰り返し回数を初期化
+ float* diff_alpha; // 双対問題の勾配値
+ float* pre_diff_alpha; // 双対問題の前回の勾配値（慣性項に用いる）
+ float* pre_alpha; // 前回の双対係数ベクトル（収束判定に用いる）
+ register float diff_sum; // 勾配計算用の小計
+ register float kernel_val; // カーネル関数とC2を含めた項
- // 双対係数の初期化.乱択
- for (int i = 0; i < n_sample; i++ ) {
- // 欠損データの係数は0にして使用しない
- if ( label[i] == 0 ) {
- alpha[i] = 0;
- continue;
- }
- alpha[i] = uniform_rand(1.0) + 1.0;
- }
+ //float plus_sum, minus_sum; // 正例と負例の係数和
- // 学習ループ
- while ( iteration < maxIteration ) {
+ // 配列（ベクトル）のアロケート
+ diff_alpha = new float[n_sample];
+ pre_diff_alpha = new float[n_sample];
+ pre_alpha = new float[n_sample];
- printf("ite: %d diff_norm : %f alpha_dist : %f \r\n", iteration, two_norm(diff_alpha, n_sample), vec_dist(alpha, pre_alpha, n_sample));
- // 前回の更新値の記録
- memcpy(pre_alpha, alpha, sizeof(float) * n_sample);
- if ( iteration >= 1 ) {
- memcpy(pre_diff_alpha, diff_alpha, sizeof(float) * n_sample);
- } else {
- // 初回は0埋めで初期化
- memset(diff_alpha, 0, sizeof(float) * n_sample);
- memset(pre_diff_alpha, 0, sizeof(float) * n_sample);
- }
+ status = SVM_NOT_LEARN; // 学習を未完了に
+ iteration = 0; // 繰り返し回数を初期化
- // 勾配値の計算
- for (int i=0; i < n_sample; i++) {
- diff_sum = 0;
- for (int j=0; j < n_sample;j++) {
- // C2を踏まえたカーネル関数値
- kernel_val = kernel_function(&(MATRIX_AT(sample,dim_signal,i,0)), &(MATRIX_AT(sample,dim_signal,j,0)), dim_signal);
- // kernel_val = MATRIX_AT(grammat,n_sample,i,j); // via Gram matrix
- if (i == j) {
- kernel_val += (1/C2);
+ // 双対係数の初期化.乱択
+ for (int i = 0; i < n_sample; i++ ) {
+ // 欠損データの係数は0にして使用しない
+ if ( label[i] == 0 ) {
+ alpha[i] = 0;
+ continue;
}
- diff_sum += alpha[j] * label[j] * kernel_val;
- }
- diff_sum *= label[i];
- diff_alpha[i] = 1 - diff_sum;
+ alpha[i] = uniform_rand(1.0) + 1.0;
}
- // 双対変数の更新
- for (int i=0; i < n_sample; i++) {
- if ( label[i] == 0 ) {
- continue;
- }
- //printf("alpha[%d] : %f -> ", i, alpha[i]);
- alpha[i] = pre_alpha[i]
- + eta * diff_alpha[i]
- + learn_alpha * pre_diff_alpha[i];
- //printf("%f \dim_signal", alpha[i]);
+ // 学習ループ
+ while ( iteration < maxIteration ) {
+
+ printf("ite: %d diff_norm : %f alpha_dist : %f \r\n", iteration, two_norm(diff_alpha, n_sample), vec_dist(alpha, pre_alpha, n_sample));
+ // 前回の更新値の記録
+ memcpy(pre_alpha, alpha, sizeof(float) * n_sample);
+ if ( iteration >= 1 ) {
+ memcpy(pre_diff_alpha, diff_alpha, sizeof(float) * n_sample);
+ } else {
+ // 初回は0埋めで初期化
+ memset(diff_alpha, 0, sizeof(float) * n_sample);
+ memset(pre_diff_alpha, 0, sizeof(float) * n_sample);
+ }
+
+ // 勾配値の計算
+ for (int i=0; i < n_sample; i++) {
+ diff_sum = 0;
+ for (int j=0; j < n_sample; j++) {
+ // C2を踏まえたカーネル関数値
+ kernel_val = kernel_function(&(MATRIX_AT(sample,dim_signal,i,0)), &(MATRIX_AT(sample,dim_signal,j,0)), dim_signal);
+ // kernel_val = MATRIX_AT(grammat,n_sample,i,j); // via Gram matrix
+ if (i == j) {
+ kernel_val += (1/C2);
+ }
+ diff_sum += alpha[j] * label[j] * kernel_val;
+ }
+ diff_sum *= label[i];
+ diff_alpha[i] = 1 - diff_sum;
+ }
+
+ // 双対変数の更新
+ for (int i=0; i < n_sample; i++) {
+ if ( label[i] == 0 ) {
+ continue;
+ }
+ //printf("alpha[%d] : %f -> ", i, alpha[i]);
+ alpha[i] = pre_alpha[i]
+ + eta * diff_alpha[i]
+ + learn_alpha * pre_diff_alpha[i];
+ //printf("%f \dim_signal", alpha[i]);
- // 非数/無限チェック
- if ( isnan(alpha[i]) || isinf(alpha[i]) ) {
- fprintf(stderr, "Detected NaN or Inf Dual-Coffience : pre_alhpa[%d]=%f -> alpha[%d]=%f", i, pre_alpha[i], i, alpha[i]);
- return SVM_DETECT_BAD_VAL;
- }
+ // 非数/無限チェック
+ if ( isnan(alpha[i]) || isinf(alpha[i]) ) {
+ fprintf(stderr, "Detected NaN or Inf Dual-Coffience : pre_alhpa[%d]=%f -> alpha[%d]=%f", i, pre_alpha[i], i, alpha[i]);
+ return SVM_DETECT_BAD_VAL;
+ }
+
+ }
- }
+ // 係数の制約条件1:正例と負例の双対係数和を等しくする.
+ // 手法:標本平均に寄せる
+ float norm_sum = 0;
+ for (int i = 0; i < n_sample; i++ ) {
+ norm_sum += (label[i] * alpha[i]);
+ }
+ norm_sum /= n_sample;
- // 係数の制約条件1:正例と負例の双対係数和を等しくする.
- // 手法:標本平均に寄せる
- float norm_sum = 0;
- for (int i = 0; i < n_sample; i++ ) {
- norm_sum += (label[i] * alpha[i]);
- }
- norm_sum /= n_sample;
+ for (int i = 0; i < n_sample; i++ ) {
+ if ( label[i] == 0 ) {
+ continue;
+ }
+ alpha[i] -= (norm_sum / label[i]);
+ }
- for (int i = 0; i < n_sample; i++ ) {
- if ( label[i] == 0 ) {
- continue;
- }
- alpha[i] -= (norm_sum / label[i]);
+ // 係数の制約条件2:双対係数は非負
+ for (int i = 0; i < n_sample; i++ ) {
+ if ( alpha[i] < 0 ) {
+ alpha[i] = 0;
+ } else if ( alpha[i] > C1 ) {
+ // C1を踏まえると,係数の上限はC1となる.
+ alpha[i] = C1;
+ }
+ }
+
+ // 収束判定 : 凸計画問題なので,収束時は大域最適が
+ // 保証されている.
+ if ( (vec_dist(alpha, pre_alpha, n_sample) < epsilon)
+ || (two_norm(diff_alpha, n_sample) < epsilon) ) {
+ // 学習の正常完了
+ status = SVM_LEARN_SUCCESS;
+ break;
+ }
+
+ // 学習繰り返し回数のインクリメント
+ iteration++;
}
- // 係数の制約条件2:双対係数は非負
- for (int i = 0; i < n_sample; i++ ) {
- if ( alpha[i] < 0 ) {
- alpha[i] = 0;
- } else if ( alpha[i] > C1 ) {
- // C1を踏まえると,係数の上限はC1となる.
- alpha[i] = C1;
- }
+ if (iteration >= maxIteration) {
+ fprintf(stderr, "Learning is not convergenced. (iteration count > maxIteration) \r\n");
+ status = SVM_NOT_CONVERGENCED;
+ } else if ( status != SVM_LEARN_SUCCESS ) {
+ status = SVM_NOT_LEARN;
}
- // 収束判定 : 凸計画問題なので,収束時は大域最適が
- // 保証されている.
- if ( (vec_dist(alpha, pre_alpha, n_sample) < epsilon)
- || (two_norm(diff_alpha, n_sample) < epsilon) ) {
- // 学習の正常完了
- status = SVM_LEARN_SUCCESS;
- break;
- }
-
- // 学習繰り返し回数のインクリメント
- iteration++;
- }
+ // 領域開放
+ delete [] diff_alpha;
+ delete [] pre_diff_alpha;
+ delete [] pre_alpha;
- if (iteration >= maxIteration) {
- fprintf(stderr, "Learning is not convergenced. (iteration count > maxIteration) \r\n");
- status = SVM_NOT_CONVERGENCED;
- } else if ( status != SVM_LEARN_SUCCESS ) {
- status = SVM_NOT_LEARN;
- }
-
- // 領域開放
- delete [] diff_alpha;
- delete [] pre_diff_alpha;
- delete [] pre_alpha;
-
- return status;
+ return status;
}
// 未知データのネットワーク値を計算
float SVM::predict_net(float* data)
{
- // 学習の終了を確認
- if (status != SVM_LEARN_SUCCESS && status != SVM_SET_ALPHA) {
- fprintf(stderr, "Learning is not completed yet.");
- //exit(1);
- return SVM_NOT_LEARN;
- }
+ // 学習の終了を確認
+ if (status != SVM_LEARN_SUCCESS && status != SVM_SET_ALPHA) {
+ fprintf(stderr, "Learning is not completed yet.");
+ //exit(1);
+ return SVM_NOT_LEARN;
+ }
- float* norm_data = new float[dim_signal];
+ float* norm_data = new float[dim_signal];
+
+ // 信号の正規化
+ for (int i = 0; i < dim_signal; i++) {
+ norm_data[i] = ( data[i] - sample_min[i] ) / ( sample_max[i] - sample_min[i] );
+ }
- // 信号の正規化
- for (int i = 0; i < dim_signal; i++) {
- norm_data[i] = ( data[i] - sample_min[i] ) / ( sample_max[i] - sample_min[i] );
- }
+ // ネットワーク値の計算
+ float net = 0;
+ for (int l=0; l < n_sample; l++) {
+ // **係数が正に相当するサンプルはサポートベクトル**
+ if(alpha[l] > 0) {
+ net += label[l] * alpha[l]
+ * kernel_function(&(MATRIX_AT(sample,dim_signal,l,0)), norm_data, dim_signal);
+ }
+ }
+
+ delete [] norm_data;
- // ネットワーク値の計算
- float net = 0;
- for (int l=0; l < n_sample; l++) {
- // **係数が正に相当するサンプルはサポートベクトル**
- if(alpha[l] > 0) {
- net += label[l] * alpha[l]
- * kernel_function(&(MATRIX_AT(sample,dim_signal,l,0)), norm_data, dim_signal);
- }
- }
-
- return net;
+ return net;
}
@@ -258,9 +261,9 @@
float* optimal_w = new float[dim_signal]; // 最適時の係数(not 双対係数)
float sigmoid_param; // シグモイド関数の温度パラメタ
float norm_w; // 係数の2乗ノルム
-
+
net = SVM::predict_net(data);
-
+
// 最適時の係数を計算
for (int n = 0; n < dim_signal; n++ ) {
optimal_w[n] = 0;
@@ -270,11 +273,9 @@
}
norm_w = two_norm(optimal_w, dim_signal);
sigmoid_param = 1 / ( norm_w * logf( (1 - epsilon) / epsilon ) );
-
+
probability = sigmoid_func(net/sigmoid_param);
-
- delete [] optimal_w;
-
+
// 打ち切り:誤差epsilon以内ならば, 1 or 0に打ち切る.
if ( probability > (1 - epsilon) ) {
return float(1);
@@ -282,23 +283,27 @@
return float(0);
}
+ delete [] optimal_w;
+
return probability;
-
+
}
// 未知データの識別
int SVM::predict_label(float* data)
{
- return (predict_net(data) >= 0) ? 1 : (-1);
+ return (predict_net(data) >= 0) ? 1 : (-1);
}
// 双対係数のゲッター
-float* SVM::get_alpha(void) {
+float* SVM::get_alpha(void)
+{
return (float *)alpha;
}
// 双対係数のセッター
-void SVM::set_alpha(float* alpha_data, int nsample) {
+void SVM::set_alpha(float* alpha_data, int nsample)
+{
if ( nsample != n_sample ) {
fprintf( stderr, " set_alpha : number of sample isn't match with arg. n_samle= %d, arg= %d \r\n", n_sample, nsample);
return;