SHOGUN
v2.0.0
|
00001 /*----------------------------------------------------------------------- 00002 * libqp_gsmo.c: implementation of the Generalized SMO algorithm. 00003 * 00004 * DESCRIPTION 00005 * The library provides function which solves the following instance of 00006 * a convex Quadratic Programming task: 00007 * 00008 * min QP(x) := 0.5*x'*H*x + f'*x 00009 * x 00010 * 00011 * s.t. a'*x = b 00012 * LB[i] <= x[i] <= UB[i] for all i=1..n 00013 * 00014 * A precision of the found solution is controlled by the input argument 00015 * TolKKT which defines tightness of the relaxed Karush-Kuhn-Tucker 00016 * stopping conditions. 00017 * 00018 * INPUT ARGUMENTS 00019 * get_col function which returns pointer to the i-th column of H. 00020 * diag_H [float64_t n x 1] vector containing values on the diagonal of H. 00021 * f [float64_t n x 1] vector. 00022 * a [float64_t n x 1] Vector which must not contain zero entries. 00023 * b [float64_t 1 x 1] Scalar. 00024 * LB [float64_t n x 1] Lower bound; -inf is allowed. 00025 * UB [float64_t n x 1] Upper bound; inf is allowed. 00026 * x [float64_t n x 1] solution vector; must be feasible. 00027 * n [uint32_t 1 x 1] dimension of H. 00028 * MaxIter [uint32_t 1 x 1] max number of iterations. 00029 * TolKKT [float64_t 1 x 1] Tightness of KKT stopping conditions. 00030 * print_state print function; if == NULL it is not called. 00031 * 00032 * RETURN VALUE 00033 * structure [libqp_state_T] 00034 * .QP [1x1] Primal objective value. 00035 * .exitflag [1 x 1] Indicates which stopping condition was used: 00036 * -3 ... initial solution vector does not satisfy equality constraint 00037 * -2 ... initial solution vector does not satisfy bounds 00038 * -1 ... not enough memory 00039 * 0 ... Maximal number of iterations reached: nIter >= MaxIter. 00040 * 4 ... Relaxed KKT conditions satisfied. 00041 * .nIter [1x1] Number of iterations. 00042 * 00043 * REFERENCE 00044 * S.S. Keerthi, E.G. Gilbert. Convergence of a generalized SMO algorithm 00045 * for SVM classier design. Technical Report CD-00-01, Control Division, 00046 * Dept. of Mechanical and Production Engineering, National University 00047 * of Singapore, 2000. 00048 * http://citeseer.ist.psu.edu/keerthi00convergence.html 00049 * 00050 * 00051 * Copyright (C) 2006-2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz 00052 * Center for Machine Perception, CTU FEL Prague 00053 * 00054 * This program is free software; you can redistribute it and/or 00055 * modify it under the terms of the GNU General Public 00056 * License as published by the Free Software Foundation; 00057 * Version 3, 29 June 2007 00058 *-------------------------------------------------------------------- */ 00059 00060 #include <math.h> 00061 #include <stdlib.h> 00062 #include <stdio.h> 00063 #include <string.h> 00064 #include <stdint.h> 00065 #include <limits.h> 00066 00067 #include <shogun/lib/common.h> 00068 #include <shogun/lib/external/libqp.h> 00069 00070 namespace shogun 00071 { 00072 00073 libqp_state_T libqp_gsmo_solver(const float64_t* (*get_col)(uint32_t), 00074 float64_t *diag_H, 00075 float64_t *f, 00076 float64_t *a, 00077 float64_t b, 00078 float64_t *LB, 00079 float64_t *UB, 00080 float64_t *x, 00081 uint32_t n, 00082 uint32_t MaxIter, 00083 float64_t TolKKT, 00084 void (*print_state)(libqp_state_T state)) 00085 { 00086 float64_t *col_u; 00087 float64_t *col_v; 00088 float64_t *Nabla; 00089 float64_t minF_up; 00090 float64_t maxF_low; 00091 float64_t tau; 00092 float64_t F_i; 00093 float64_t tau_ub, tau_lb; 00094 uint32_t i, j; 00095 uint32_t u=0, v=0; 00096 libqp_state_T state; 00097 float64_t atx = 0.0; 00098 00099 Nabla = NULL; 00100 00101 /* ------------------------------------------------------------ */ 00102 /* Initialization */ 00103 /* ------------------------------------------------------------ */ 00104 00105 // check bounds of initial guess 00106 for (i=0; i<n; i++) 00107 { 00108 if (x[i]>UB[i]) 00109 { 00110 state.exitflag = -2; 00111 goto cleanup; 00112 } 00113 if (x[i]<LB[i]) 00114 { 00115 state.exitflag = -2; 00116 goto cleanup; 00117 } 00118 } 00119 00120 // check equality constraint 00121 for (i=0; i<n; i++) 00122 atx += a[i]*x[i]; 00123 if (fabs(b-atx)>1e-9) 00124 { 00125 printf("%f \ne %f\n",b,atx); 00126 state.exitflag = -3; 00127 goto cleanup; 00128 } 00129 00130 /* Nabla = H*x + f is gradient*/ 00131 Nabla = (float64_t*)LIBQP_CALLOC(n, sizeof(float64_t)); 00132 if( Nabla == NULL ) 00133 { 00134 state.exitflag=-1; 00135 goto cleanup; 00136 } 00137 00138 /* compute gradient */ 00139 for( i=0; i < n; i++ ) 00140 { 00141 Nabla[i] += f[i]; 00142 if( x[i] != 0 ) { 00143 col_u = (float64_t*)get_col(i); 00144 for( j=0; j < n; j++ ) { 00145 Nabla[j] += col_u[j]*x[i]; 00146 } 00147 } 00148 } 00149 00150 if( print_state != NULL) 00151 { 00152 state.QP = 0; 00153 for(i = 0; i < n; i++ ) 00154 state.QP += 0.5*(x[i]*Nabla[i]+x[i]*f[i]); 00155 00156 print_state( state ); 00157 } 00158 00159 00160 /* ------------------------------------------------------------ */ 00161 /* Main optimization loop */ 00162 /* ------------------------------------------------------------ */ 00163 00164 state.nIter = 0; 00165 state.exitflag = 100; 00166 while( state.exitflag == 100 ) 00167 { 00168 state.nIter ++; 00169 00170 /* find the most violating pair of variables */ 00171 minF_up = LIBQP_PLUS_INF; 00172 maxF_low = -LIBQP_PLUS_INF; 00173 for(i = 0; i < n; i++ ) 00174 { 00175 00176 F_i = Nabla[i]/a[i]; 00177 00178 if(LB[i] < x[i] && x[i] < UB[i]) 00179 { /* i is from I_0 */ 00180 if( minF_up > F_i) { minF_up = F_i; u = i; } 00181 if( maxF_low < F_i) { maxF_low = F_i; v = i; } 00182 } 00183 else if((a[i] > 0 && x[i] == LB[i]) || (a[i] < 0 && x[i] == UB[i])) 00184 { /* i is from I_1 or I_2 */ 00185 if( minF_up > F_i) { minF_up = F_i; u = i; } 00186 } 00187 else if((a[i] > 0 && x[i] == UB[i]) || (a[i] < 0 && x[i] == LB[i])) 00188 { /* i is from I_3 or I_4 */ 00189 if( maxF_low < F_i) { maxF_low = F_i; v = i; } 00190 } 00191 } 00192 00193 /* check KKT conditions */ 00194 if( maxF_low - minF_up <= TolKKT ) 00195 state.exitflag = 4; 00196 else 00197 { 00198 /* SMO update of the most violating pair */ 00199 col_u = (float64_t*)get_col(u); 00200 col_v = (float64_t*)get_col(v); 00201 00202 if( a[u] > 0 ) 00203 { tau_lb = (LB[u]-x[u])*a[u]; tau_ub = (UB[u]-x[u])*a[u]; } 00204 else 00205 { tau_ub = (LB[u]-x[u])*a[u]; tau_lb = (UB[u]-x[u])*a[u]; } 00206 00207 if( a[v] > 0 ) 00208 { tau_lb = LIBQP_MAX(tau_lb,(x[v]-UB[v])*a[v]); tau_ub = LIBQP_MIN(tau_ub,(x[v]-LB[v])*a[v]); } 00209 else 00210 { tau_lb = LIBQP_MAX(tau_lb,(x[v]-LB[v])*a[v]); tau_ub = LIBQP_MIN(tau_ub,(x[v]-UB[v])*a[v]); } 00211 00212 tau = (Nabla[v]/a[v]-Nabla[u]/a[u])/ 00213 (diag_H[u]/(a[u]*a[u]) + diag_H[v]/(a[v]*a[v]) - 2*col_u[v]/(a[u]*a[v])); 00214 00215 tau = LIBQP_MIN(LIBQP_MAX(tau,tau_lb),tau_ub); 00216 00217 x[u] += tau/a[u]; 00218 x[v] -= tau/a[v]; 00219 00220 /* update Nabla */ 00221 for(i = 0; i < n; i++ ) 00222 Nabla[i] += col_u[i]*tau/a[u] - col_v[i]*tau/a[v]; 00223 00224 } 00225 00226 if( state.nIter >= MaxIter ) 00227 state.exitflag = 0; 00228 00229 if( print_state != NULL) 00230 { 00231 state.QP = 0; 00232 for(i = 0; i < n; i++ ) 00233 state.QP += 0.5*(x[i]*Nabla[i]+x[i]*f[i]); 00234 00235 print_state( state ); 00236 } 00237 00238 } 00239 00240 /* compute primal objective value */ 00241 state.QP = 0; 00242 for(i = 0; i < n; i++ ) 00243 state.QP += 0.5*(x[i]*Nabla[i]+x[i]*f[i]); 00244 00245 cleanup: 00246 00247 LIBQP_FREE(Nabla); 00248 00249 return( state ); 00250 } 00251 00252 } /* shogun namespace */ 00253