00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070 #include <math.h>
00071 #include <stdlib.h>
00072 #include <stdio.h>
00073 #include <string.h>
00074 #include <stdint.h>
00075 #include <limits.h>
00076
00077 #include "classifier/svm/libocas_common.h"
00078 #include "classifier/svm/qpssvmlib.h"
00079
00080 namespace shogun
00081 {
00082
00083
00084
00085
00086
00087
00088 int8_t qpssvm_solver(const void* (*get_col)(uint32_t),
00089 float64_t *diag_H,
00090 float64_t *f,
00091 float64_t b,
00092 uint16_t *I,
00093 float64_t *x,
00094 uint32_t n,
00095 uint32_t tmax,
00096 float64_t tolabs,
00097 float64_t tolrel,
00098 float64_t *QP,
00099 float64_t *QD,
00100 uint32_t verb)
00101 {
00102 float64_t *x_nequ;
00103 float64_t *d;
00104 float64_t *col_u, *col_v;
00105 float64_t LB;
00106 float64_t UB;
00107 float64_t tmp;
00108 float64_t improv;
00109 float64_t tmp_num;
00110 float64_t tmp_den=0;
00111 float64_t tau=0;
00112 float64_t delta;
00113 float64_t yu;
00114 uint32_t *inx;
00115 uint32_t *nk;
00116 uint32_t m;
00117 uint32_t t;
00118 uint32_t u=0;
00119 uint32_t v=0;
00120 uint32_t k;
00121 uint32_t i, j;
00122 int8_t exitflag;
00123
00124
00125
00126
00127
00128
00129 x_nequ=NULL;
00130 inx=NULL;
00131 nk=NULL;
00132 d=NULL;
00133
00134
00135 for( i=0, m=0; i < n; i++ ) m = CMath::max(m,(uint32_t) I[i]);
00136
00137
00138 x_nequ = (float64_t*) OCAS_CALLOC(m, sizeof(float64_t));
00139 if( x_nequ == NULL )
00140 {
00141 exitflag=-2;
00142 goto cleanup;
00143 }
00144
00145
00146 inx = (uint32_t*) OCAS_CALLOC(m*n, sizeof(uint32_t));
00147 if( inx == NULL )
00148 {
00149 exitflag=-2;
00150 goto cleanup;
00151 }
00152
00153 nk = (uint32_t*) OCAS_CALLOC(m, sizeof(uint32_t));
00154 if( nk == NULL )
00155 {
00156 exitflag=-2;
00157 goto cleanup;
00158 }
00159
00160 for( i=0; i < m; i++ ) x_nequ[i] = b;
00161 for( i=0; i < n; i++ ) {
00162 k = I[i]-1;
00163 x_nequ[k] -= x[i];
00164 inx[INDEX2(nk[k],k,n)] = i;
00165 nk[k]++;
00166 }
00167
00168
00169 d = (float64_t*) OCAS_CALLOC(n, sizeof(float64_t));
00170 if( d == NULL )
00171 {
00172 exitflag=-2;
00173 goto cleanup;
00174 }
00175
00176
00177 for( i=0; i < n; i++ ) {
00178 if( x[i] > 0 ) {
00179 col_u = (float64_t*)get_col(i);
00180 for( j=0; j < n; j++ ) {
00181 d[j] += col_u[j]*x[i];
00182 }
00183 }
00184 }
00185 for( i=0; i < n; i++ ) d[i] += f[i];
00186
00187
00188
00189 for( i=0, UB = 0, LB=0; i < n; i++) {
00190 UB += x[i]*(f[i]+d[i]);
00191 LB += x[i]*(f[i]-d[i]);
00192 }
00193 UB = 0.5*UB;
00194 LB = 0.5*LB;
00195
00196
00197
00198
00199
00200
00201
00202
00203 for( i=0; i < m; i++ ) {
00204 for( j=0, tmp = OCAS_PLUS_INF; j < nk[i]; j++ ) {
00205 tmp = CMath::min(tmp, d[inx[INDEX2(j,i,n)]]);
00206 }
00207 if( tmp < 0) LB += b*tmp;
00208 }
00209
00210 exitflag = 0;
00211 t = 0;
00212
00213
00214 while( (exitflag == 0) && (t < tmax))
00215 {
00216 t++;
00217
00218 exitflag = 1;
00219 for( k=0; k < m; k++ )
00220 {
00221
00222
00223
00224
00225
00226 for( j=0, tmp = OCAS_PLUS_INF, delta = 0; j < nk[k]; j++ ) {
00227 i = inx[INDEX2(j,k,n)];
00228 delta += x[i]*d[i];
00229 if( tmp > d[i] ) {
00230 tmp = d[i];
00231 u = i;
00232 }
00233 }
00234
00235
00236 if( d[u] < 0) yu = b; else yu = 0;
00237
00238
00239 delta -= yu*d[u];
00240
00241 if( delta > tolabs/m && delta > tolrel*CMath::abs(UB)/m)
00242 {
00243 exitflag = 0;
00244
00245 if( yu > 0 )
00246 {
00247 col_u = (float64_t*)get_col(u);
00248
00249 improv = -OCAS_PLUS_INF;
00250 for( j=0; j < nk[k]; j++ ) {
00251 i = inx[INDEX2(j,k,n)];
00252
00253
00254
00255 if(x[i] > 0) {
00256
00257 tmp_num = x[i]*(d[i] - d[u]);
00258 tmp_den = x[i]*x[i]*(diag_H[u] - 2*col_u[i] + diag_H[i]);
00259 if( tmp_den > 0 ) {
00260 if( tmp_num < tmp_den ) {
00261 tmp = tmp_num*tmp_num / tmp_den;
00262 } else {
00263 tmp = tmp_num - 0.5 * tmp_den;
00264 }
00265 }
00266 if( tmp > improv ) {
00267 improv = tmp;
00268 tau = CMath::min(1.0,tmp_num/tmp_den);
00269 v = i;
00270 }
00271 }
00272 }
00273
00274 tmp_num = -x_nequ[k]*d[u];
00275 if( tmp_num > 0 ) {
00276 tmp_den = x_nequ[k]*x_nequ[k]*diag_H[u];
00277 if( tmp_den > 0 ) {
00278 if( tmp_num < tmp_den ) {
00279 tmp = tmp_num*tmp_num / tmp_den;
00280 } else {
00281 tmp = tmp_num - 0.5 * tmp_den;
00282 }
00283 }
00284 } else {
00285 tmp = -OCAS_PLUS_INF;
00286 }
00287
00288 if( tmp > improv ) {
00289 tau = CMath::min(1.0,tmp_num/tmp_den);
00290 for( i = 0; i < n; i++ ) {
00291 d[i] += x_nequ[k]*tau*col_u[i];
00292 }
00293 x[u] += tau*x_nequ[k];
00294 x_nequ[k] -= tau*x_nequ[k];
00295
00296 } else {
00297
00298
00299 col_v = (float64_t*)get_col(v);
00300 for( i = 0; i < n; i++ ) {
00301 d[i] += x[v]*tau*(col_u[i]-col_v[i]);
00302 }
00303
00304 x[u] += tau*x[v];
00305 x[v] -= tau*x[v];
00306 }
00307 }
00308 else
00309 {
00310 improv = -OCAS_PLUS_INF;
00311 for( j=0; j < nk[k]; j++ ) {
00312 i = inx[INDEX2(j,k,n)];
00313
00314
00315
00316 if( x[i] > 0 && d[i] > 0) {
00317
00318 tmp_num = x[i]*d[i];
00319 tmp_den = x[i]*x[i]*diag_H[i];
00320 if( tmp_den > 0 ) {
00321 if( tmp_num < tmp_den ) {
00322 tmp = tmp_num*tmp_num / tmp_den;
00323 } else {
00324 tmp = tmp_num - 0.5 * tmp_den;
00325 }
00326 }
00327 if( tmp > improv ) {
00328 improv = tmp;
00329 tau = CMath::min(1.0,tmp_num/tmp_den);
00330 v = i;
00331 }
00332 }
00333 }
00334
00335
00336 col_v = (float64_t*)get_col(v);
00337 for( i = 0; i < n; i++ ) {
00338 d[i] -= x[v]*tau*col_v[i];
00339 }
00340
00341 x_nequ[k] += tau*x[v];
00342 x[v] -= tau*x[v];
00343 }
00344
00345 UB = UB - improv;
00346 }
00347
00348 }
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358 for( i=0, UB = 0, LB=0; i < n; i++) {
00359 UB += x[i]*(f[i]+d[i]);
00360 LB += x[i]*(f[i]-d[i]);
00361 }
00362 UB = 0.5*UB;
00363 LB = 0.5*LB;
00364
00365 for( k=0; k < m; k++ ) {
00366 for( j=0,tmp = OCAS_PLUS_INF; j < nk[k]; j++ ) {
00367 i = inx[INDEX2(j,k,n)];
00368
00369 tmp = CMath::min(tmp, d[i]);
00370 }
00371 if( tmp < 0) LB += b*tmp;
00372 }
00373
00374 if( verb > 0 && (exitflag > 0 || (t % verb)==0 ))
00375 {
00376 float64_t gap=(UB!=0) ? (UB-LB)/CMath::abs(UB) : 0;
00377 SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(tolrel), 6);
00378 }
00379
00380 }
00381
00382
00383 if( UB-LB < tolabs ) exitflag = 1;
00384 else if(UB-LB < CMath::abs(UB)*tolrel ) exitflag = 2;
00385 else exitflag = 0;
00386
00387
00388
00389
00390 *QP = UB;
00391 *QD = LB;
00392
00393
00394
00395
00396 cleanup:
00397 OCAS_FREE( d );
00398 OCAS_FREE( inx );
00399 OCAS_FREE( nk );
00400 OCAS_FREE( x_nequ );
00401
00402 return( exitflag );
00403
00404 }
00405 }