?? qpcsolver.c
字號:
aHa = tmp_aHa2;
ac = tmp_ac2;
/*
Ha = Ha + gamma2*alpha(v)*(H(:,u)-H(:,v));
alpha = alpha*(1-gamma2)+gamma2*alpha2;
*/
tmp = alpha[v]*gamma2;
for( i = 0; i < dim; i++ ) {
Ha[i] = Ha[i] + tmp*(col_u[i] - col_v[i]);
}
alpha[u] = alpha[u] + tmp;
alpha[v] = alpha[v] - tmp;
break;
case 3:
aHa = tmp_aHa3;
ac = tmp_ac3;
/*
Ha = gamma3*Ha + H(:,u)*(1-gamma3+gamma3*alpha(v)) - ...
gamma3*alpha(v)*H(:,v);
alpha = alpha1*(1-gamma3)+gamma3*alpha2;
*/
tmp = alpha[v]*gamma3;
for( i = 0; i < dim; i++ ) {
Ha[i] = gamma3*Ha[i] + col_u[i]*(1-gamma3+tmp) - tmp*col_v[i];
alpha[i] = alpha[i]*gamma3;
}
alpha[u] = alpha[u] + 1 - gamma3 + tmp;
alpha[v] = alpha[v] - tmp;
break;
}
}
UB = 0.5*aHa + ac;
min_beta = PLUS_INF;
max_beta = MINUS_INF;
for( i = 0; i < dim; i++ )
{
beta = Ha[i]+ vector_c[i];
if( alpha[i] !=0 && max_beta < beta )
{
v = i;
max_beta = beta;
}
if( beta < min_beta )
{
u = i;
min_beta = beta;
}
}
LB = min_beta - 0.5*aHa;
/* Stopping conditions */
if( UB-LB <= tolabs ) exitflag = 1;
else if( UB-LB <= ABS(UB)*tolrel ) exitflag = 2;
else if(t >= tmax) exitflag = 0;
/* print info */
if(verb && (t % verb) == 0 ) {
mexPrintf("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
t, UB, LB, UB-LB,(UB-LB)/UB);
}
/* Store selected values */
if( t < History_size ) {
History[INDEX(0,t,2)] = LB;
History[INDEX(1,t,2)] = UB;
}
else {
tmp_ptr = mxCalloc((History_size+HISTORY_BUF)*2,sizeof(double));
if( tmp_ptr == NULL ) mexErrMsgTxt("Not enough memory.");
for( i = 0; i < History_size; i++ ) {
tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
}
tmp_ptr[INDEX(0,t,2)] = LB;
tmp_ptr[INDEX(1,t,2)] = UB;
History_size += HISTORY_BUF;
mxFree( History );
History = tmp_ptr;
}
}
/* print info about last iteration*/
if(verb && (t % verb) ) {
mexPrintf("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
t, UB, LB, UB-LB,(UB-LB)/UB);
}
/*------------------------------------------------------- */
/* Set outputs */
/*------------------------------------------------------- */
(*ptr_t) = t;
(*ptr_History) = History;
/* Free memory */
mxFree( Ha );
return( exitflag );
}
/* --------------------------------------------------------------
QPC solver based on the Kowalczyk's algorithm
(Maximal Margin Perceptron).
Usage: exitflag = qpc_kowalczyk( &get_col, diag_H, vector_c, dim,
tmax, tolabs, tolrel, &alpha, &t, &History );
-------------------------------------------------------------- */
int qpc_kowalczyk(const void* (*get_col)(long,long),
double *diag_H,
double *vector_c,
long dim,
long tmax,
double tolabs,
double tolrel,
double *alpha,
long *ptr_t,
double **ptr_History,
long verb)
{
double LB;
double UB;
double aHa, ac;
double tmp, tmp1;
double tmp_aHa, tmp_ac, tmp_UB, tmp_gamma;
double min_beta, beta;
double x10, x20, x11, x12, x22;
double delta;
double gamma;
double *History;
double *Ha;
double *tmp_ptr;
double *col_inx;
long i;
long inx;
long t;
long History_size;
int exitflag;
/* ------------------------------------------------------------ */
/* Initialization */
/* ------------------------------------------------------------ */
Ha = mxCalloc(dim, sizeof(double));
if( Ha == NULL ) mexErrMsgTxt("Not enough memory.");
History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
History = mxCalloc(History_size*2,sizeof(double));
if( History == NULL ) mexErrMsgTxt("Not enough memory.");
/* inx = argmin(0.5*diag_H + vector_c ); */
for( tmp1 = PLUS_INF, i = 0; i < dim; i++ ) {
tmp = 0.5*diag_H[i] + vector_c[i];
if( tmp1 > tmp) {
tmp1 = tmp;
inx = i;
}
}
col_inx = (double*)get_col(inx,-1);
for( min_beta = PLUS_INF, i = 0; i < dim; i++ )
{
alpha[i] = 0;
Ha[i] = col_inx[i];
beta = Ha[i] + vector_c[i];
if( beta < min_beta ) {
min_beta = beta;
}
}
alpha[inx] = 1;
aHa = diag_H[inx];
ac = vector_c[inx];
UB = 0.5*aHa + ac;
LB = min_beta - 0.5*aHa;
t = 0;
History[INDEX(0,0,2)] = LB;
History[INDEX(1,0,2)] = UB;
if( verb ) {
mexPrintf("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
UB, LB, UB-LB,(UB-LB)/UB);
}
/* Stopping conditions */
if( UB-LB <= tolabs ) exitflag = 1;
else if(UB-LB <= ABS(UB)*tolrel ) exitflag = 2;
else exitflag = -1;
/* ------------------------------------------------------------ */
/* Main optimization loop */
/* ------------------------------------------------------------ */
while( exitflag == -1 )
{
t++;
x11 = aHa;
x10 = ac;
/* searches for rule which yileds the biggest improvement */
for( i = 0; i < dim; i++ )
{
delta = Ha[i] + vector_c[i] - aHa - ac;
tmp_UB = PLUS_INF;
if( delta < 0 )
{
/* Kozinec rule */
x12 = Ha[i];
x20 = vector_c[i];
x22 = diag_H[i];
tmp_gamma = (x11 - x12 + x10 - x20)/(x11 - 2*x12 + x22);
tmp_gamma = MIN(1,tmp_gamma);
tmp_aHa = (1-tmp_gamma)*(1-tmp_gamma)*x11+2*(1-tmp_gamma)*tmp_gamma*x12
+tmp_gamma*tmp_gamma*x22;
tmp_ac = (1-tmp_gamma)*x10+tmp_gamma*x20;
tmp_UB = 0.5*tmp_aHa + tmp_ac;
}
else if( delta > 0 && alpha[i] < 1 && alpha[i] > 0)
{
x12 = (x11-alpha[i]*Ha[i])/(1-alpha[i]);
x22 = (x11-2*alpha[i]*Ha[i]+alpha[i]*alpha[i]*diag_H[i])/
((1-alpha[i])*(1-alpha[i]));
x20 = (x10-alpha[i]*vector_c[i])/(1-alpha[i]);
tmp_gamma = (x11 - x12 + x10 - x20)/(x11 - 2*x12 + x22);
tmp_gamma = MIN(1,tmp_gamma);
tmp_aHa = (1-tmp_gamma)*(1-tmp_gamma)*x11+2*(1-tmp_gamma)*tmp_gamma*x12
+ tmp_gamma*tmp_gamma*x22;
tmp_ac = (1-tmp_gamma)*x10+tmp_gamma*x20;
tmp_UB = 0.5*tmp_aHa + tmp_ac;
}
if( tmp_UB < UB )
{
UB = tmp_UB;
gamma = tmp_gamma;
aHa = tmp_aHa;
ac = tmp_ac;
inx = i;
}
}
col_inx = (double*)get_col(inx,-1);
/* Use the update with biggest improvement */
delta = Ha[inx] + vector_c[inx] - x11 - x10;
if( delta < 0 )
{
/* Kozinec rule */
for(i = 0; i < dim; i++ ) {
Ha[i] = Ha[i]*(1-gamma) + gamma*col_inx[i];
alpha[i] = alpha[i]*(1-gamma);
}
alpha[inx] = alpha[inx] + gamma;
}
else
{
/* Inverse Kozinec rule */
tmp = gamma*alpha[inx];
tmp1 = 1-alpha[inx];
for(i = 0; i < dim; i++ ) {
Ha[i] = (Ha[i]*(tmp1+tmp) - tmp*col_inx[i])/tmp1;
alpha[i] = alpha[i]*(1-gamma) + gamma*alpha[i]/tmp1;
}
alpha[inx] = alpha[inx] - tmp/tmp1;
}
min_beta = PLUS_INF;
for( i = 0; i < dim; i++ )
{
beta = Ha[i]+ vector_c[i];
if( beta < min_beta ) { min_beta = beta; }
}
LB = min_beta - 0.5*aHa;
/* Stopping conditions */
if( UB-LB <= tolabs ) exitflag = 1;
else if( UB-LB <= ABS(UB)*tolrel ) exitflag = 2;
else if(t >= tmax) exitflag = 0;
/* print info */
if(verb && (t % verb) == 0 ) {
mexPrintf("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
t, UB, LB, UB-LB,(UB-LB)/UB);
}
/* Store selected values */
if( t < History_size ) {
History[INDEX(0,t,2)] = LB;
History[INDEX(1,t,2)] = UB;
}
else {
tmp_ptr = mxCalloc((History_size+HISTORY_BUF)*2,sizeof(double));
if( tmp_ptr == NULL ) mexErrMsgTxt("Not enough memory.");
for( i = 0; i < History_size; i++ ) {
tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
}
tmp_ptr[INDEX(0,t,2)] = LB;
tmp_ptr[INDEX(1,t,2)] = UB;
History_size += HISTORY_BUF;
mxFree( History );
History = tmp_ptr;
}
}
/* print info about last iteration*/
if(verb && (t % verb) ) {
mexPrintf("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
t, UB, LB, UB-LB,(UB-LB)/UB);
}
/*------------------------------------------------------- */
/* Set outputs */
/*------------------------------------------------------- */
(*ptr_t) = t;
(*ptr_History) = History;
/* Free memory */
mxFree( Ha );
return( exitflag );
}
/* --------------------------------------------------------------
QPC solver based on the Kozinec algorithm.
Usage: exitflag = qpc_kozinec( &get_col, diag_H, vector_c, dim,
tmax, tolabs, tolrel, &alpha, &t, &History );
-------------------------------------------------------------- */
int qpc_kozinec(const void* (*get_col)(long,long),
double *diag_H,
double *vector_c,
long dim,
long tmax,
double tolabs,
double tolrel,
double *alpha,
long *ptr_t,
double **ptr_History,
long verb)
{
double LB;
double UB;
double aHa, ac;
double tmp, tmp1;
double min_beta, beta;
double lambda;
double *History;
double *Ha;
double *tmp_ptr;
double *col_u;
long u, inx;
long new_u;
long i;
long t;
long History_size;
int exitflag;
/* ------------------------------------------------------------ */
/* Initialization */
/* ------------------------------------------------------------ */
Ha = mxCalloc(dim, sizeof(double));
if( Ha == NULL ) mexErrMsgTxt("Not enough memory.");
History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
History = mxCalloc(History_size*2,sizeof(double));
if( History == NULL ) mexErrMsgTxt("Not enough memory.");
/* inx = argmin(0.5*diag_H + vector_c ); */
for( tmp1 = PLUS_INF, i = 0; i < dim; i++ ) {
tmp = 0.5*diag_H[i] + vector_c[i];
if( tmp1 > tmp) {
tmp1 = tmp;
inx = i;
}
}
col_u = (double*)get_col(inx,-1);
for( min_beta = PLUS_INF, i = 0; i < dim; i++ )
{
alpha[i] = 0;
Ha[i] = col_u[i];
beta = Ha[i] + vector_c[i];
if( beta < min_beta ) {
min_beta = beta;
u = i;
}
}
alpha[inx] = 1;
aHa = diag_H[inx];
ac = vector_c[inx];
UB = 0.5*aHa + ac;
LB = min_beta - 0.5*aHa;
t = 0;
History[INDEX(0,0,2)] = LB;
History[INDEX(1,0,2)] = UB;
if( verb ) {
mexPrintf("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
UB, LB, UB-LB,(UB-LB)/UB);
}
/* Stopping conditions */
if( UB-LB <= tolabs ) exitflag = 1;
else if(UB-LB <= ABS(UB)*tolrel ) exitflag = 2;
else exitflag = -1;
/* ------------------------------------------------------------ */
/* Main optimization loop */
/* ------------------------------------------------------------ */
while( exitflag == -1 )
{
t++;
col_u = (double*)get_col(u,-1);
/* Adaptation rule and update */
lambda = (aHa - Ha[u] + ac - vector_c[u])/(aHa - 2*Ha[u] + diag_H[u]);
lambda = MIN(1,lambda);
aHa = aHa*(1-lambda)*(1-lambda) + 2*lambda*(1-lambda)*Ha[u]
+ diag_H[u]*lambda*lambda;
ac = ac*(1-lambda) + lambda*vector_c[u];
min_beta = PLUS_INF;
for( i = 0; i < dim; i++ ) {
alpha[i] = alpha[i]*(1-lambda);
Ha[i] = Ha[i]*(1-lambda) + lambda*col_u[i];
beta = Ha[i] + vector_c[i];
if( min_beta > beta ) {
min_beta = beta;
new_u = i;
}
}
alpha[u] = alpha[u] + lambda;
UB = 0.5*aHa + ac;
LB = min_beta - 0.5*aHa;
u = new_u;
/* Stopping conditions */
if( UB-LB <= tolabs ) exitflag = 1;
else if( UB-LB <= ABS(UB)*tolrel ) exitflag = 2;
else if(t >= tmax) exitflag = 0;
/* print info about last iteration*/
if(verb && (t % verb) == 0 ) {
mexPrintf("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
t, UB, LB, UB-LB,(UB-LB)/UB);
}
/* Store selected values */
if( t < History_size ) {
History[INDEX(0,t,2)] = LB;
History[INDEX(1,t,2)] = UB;
}
else {
tmp_ptr = mxCalloc((History_size+HISTORY_BUF)*2,sizeof(double));
if( tmp_ptr == NULL ) mexErrMsgTxt("Not enough memory.");
for( i = 0; i < History_size; i++ ) {
tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
}
tmp_ptr[INDEX(0,t,2)] = LB;
tmp_ptr[INDEX(1,t,2)] = UB;
History_size += HISTORY_BUF;
mxFree( History );
History = tmp_ptr;
}
}
/* print info about last iteration*/
if(verb && (t % verb) ) {
mexPrintf("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
t, UB, LB, UB-LB,(UB-LB)/UB);
}
/*------------------------------------------------------- */
/* Set outputs */
/*------------------------------------------------------- */
(*ptr_t) = t;
(*ptr_History) = History;
/* Free memory */
mxFree( Ha );
return( exitflag );
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -