| void Solver::Solve | ( | int | l, | |
| const Kernel & | Q, | |||
| const double * | b_, | |||
| const schar * | y_, | |||
| double * | alpha_, | |||
| double | Cp, | |||
| double | Cn, | |||
| double | eps, | |||
| SolutionInfo * | si, | |||
| int | shrinking | |||
| ) |
Definition at line 441 of file svm.cpp.
References active_set, active_size, alpha, alpha_status, b, calculate_rho(), clone(), do_shrinking(), G, G_bar, get_C(), Kernel::get_Q(), info(), info_flush(), is_lower_bound(), is_upper_bound(), max, min, Solver::SolutionInfo::obj, Q, reconstruct_gradient(), Solver::SolutionInfo::rho, select_working_set(), unshrinked, update_alpha_status(), Solver::SolutionInfo::upper_bound_n, Solver::SolutionInfo::upper_bound_p, and y.
Referenced by Solver_NU::Solve(), solve_c_svc(), solve_epsilon_svr(), and solve_one_class().
00444 { 00445 this->l = l; 00446 this->Q = &Q; 00447 clone(b, b_,l); 00448 clone(y, y_,l); 00449 clone(alpha,alpha_,l); 00450 this->Cp = Cp; 00451 this->Cn = Cn; 00452 this->eps = eps; 00453 unshrinked = false; 00454 00455 // initialize alpha_status 00456 { 00457 alpha_status = new char[l]; 00458 for(int i=0;i<l;i++) 00459 update_alpha_status(i); 00460 } 00461 00462 // initialize active set (for shrinking) 00463 { 00464 active_set = new int[l]; 00465 for(int i=0;i<l;i++) 00466 active_set[i] = i; 00467 active_size = l; 00468 } 00469 00470 // initialize gradient 00471 { 00472 G = new double[l]; 00473 G_bar = new double[l]; 00474 int i; 00475 for(i=0;i<l;i++) 00476 { 00477 G[i] = b[i]; 00478 G_bar[i] = 0; 00479 } 00480 for(i=0;i<l;i++) 00481 if(!is_lower_bound(i)) 00482 { 00483 Qfloat *Q_i = Q.get_Q(i,l); 00484 double alpha_i = alpha[i]; 00485 int j; 00486 for(j=0;j<l;j++) 00487 G[j] += alpha_i*Q_i[j]; 00488 if(is_upper_bound(i)) 00489 for(j=0;j<l;j++) 00490 G_bar[j] += get_C(i) * Q_i[j]; 00491 } 00492 } 00493 00494 // optimization step 00495 00496 int iter = 0; 00497 int counter = min(l,1000)+1; 00498 00499 while(1) 00500 { 00501 // show progress and do shrinking 00502 00503 if(--counter == 0) 00504 { 00505 counter = min(l,1000); 00506 if(shrinking) do_shrinking(); 00507 info("."); info_flush(); 00508 } 00509 00510 int i,j; 00511 if(select_working_set(i,j)!=0) 00512 { 00513 // reconstruct the whole gradient 00514 reconstruct_gradient(); 00515 // reset active set size and check 00516 active_size = l; 00517 info("*"); info_flush(); 00518 if(select_working_set(i,j)!=0) 00519 break; 00520 else 00521 counter = 1; // do shrinking next iteration 00522 } 00523 00524 ++iter; 00525 00526 // update alpha[i] and alpha[j], handle bounds carefully 00527 00528 const Qfloat *Q_i = Q.get_Q(i,active_size); 00529 const Qfloat *Q_j = Q.get_Q(j,active_size); 00530 00531 double C_i = get_C(i); 00532 double C_j = get_C(j); 00533 00534 double old_alpha_i = alpha[i]; 00535 double old_alpha_j = alpha[j]; 00536 00537 if(y[i]!=y[j]) 00538 { 00539 double delta = (-G[i]-G[j])/std::max(Q_i[i]+Q_j[j]+2*Q_i[j],(Qfloat)0); 00540 double diff = alpha[i] - alpha[j]; 00541 alpha[i] += delta; 00542 alpha[j] += delta; 00543 00544 if(diff > 0) 00545 { 00546 if(alpha[j] < 0) 00547 { 00548 alpha[j] = 0; 00549 alpha[i] = diff; 00550 } 00551 } 00552 else 00553 { 00554 if(alpha[i] < 0) 00555 { 00556 alpha[i] = 0; 00557 alpha[j] = -diff; 00558 } 00559 } 00560 if(diff > C_i - C_j) 00561 { 00562 if(alpha[i] > C_i) 00563 { 00564 alpha[i] = C_i; 00565 alpha[j] = C_i - diff; 00566 } 00567 } 00568 else 00569 { 00570 if(alpha[j] > C_j) 00571 { 00572 alpha[j] = C_j; 00573 alpha[i] = C_j + diff; 00574 } 00575 } 00576 } 00577 else 00578 { 00579 double delta = (G[i]-G[j])/std::max(Q_i[i]+Q_j[j]-2*Q_i[j],(Qfloat)0); 00580 double sum = alpha[i] + alpha[j]; 00581 alpha[i] -= delta; 00582 alpha[j] += delta; 00583 if(sum > C_i) 00584 { 00585 if(alpha[i] > C_i) 00586 { 00587 alpha[i] = C_i; 00588 alpha[j] = sum - C_i; 00589 } 00590 } 00591 else 00592 { 00593 if(alpha[j] < 0) 00594 { 00595 alpha[j] = 0; 00596 alpha[i] = sum; 00597 } 00598 } 00599 if(sum > C_j) 00600 { 00601 if(alpha[j] > C_j) 00602 { 00603 alpha[j] = C_j; 00604 alpha[i] = sum - C_j; 00605 } 00606 } 00607 else 00608 { 00609 if(alpha[i] < 0) 00610 { 00611 alpha[i] = 0; 00612 alpha[j] = sum; 00613 } 00614 } 00615 } 00616 00617 // update G 00618 00619 double delta_alpha_i = alpha[i] - old_alpha_i; 00620 double delta_alpha_j = alpha[j] - old_alpha_j; 00621 00622 for(int k=0;k<active_size;k++) 00623 { 00624 G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j; 00625 } 00626 00627 // update alpha_status and G_bar 00628 00629 { 00630 bool ui = is_upper_bound(i); 00631 bool uj = is_upper_bound(j); 00632 update_alpha_status(i); 00633 update_alpha_status(j); 00634 int k; 00635 if(ui != is_upper_bound(i)) 00636 { 00637 Q_i = Q.get_Q(i,l); 00638 if(ui) 00639 for(k=0;k<l;k++) 00640 G_bar[k] -= C_i * Q_i[k]; 00641 else 00642 for(k=0;k<l;k++) 00643 G_bar[k] += C_i * Q_i[k]; 00644 } 00645 00646 if(uj != is_upper_bound(j)) 00647 { 00648 Q_j = Q.get_Q(j,l); 00649 if(uj) 00650 for(k=0;k<l;k++) 00651 G_bar[k] -= C_j * Q_j[k]; 00652 else 00653 for(k=0;k<l;k++) 00654 G_bar[k] += C_j * Q_j[k]; 00655 } 00656 } 00657 } 00658 00659 // calculate rho 00660 00661 si->rho = calculate_rho(); 00662 00663 // calculate objective value 00664 { 00665 double v = 0; 00666 int i; 00667 for(i=0;i<l;i++) 00668 v += alpha[i] * (G[i] + b[i]); 00669 00670 si->obj = v/2; 00671 } 00672 00673 // put back the solution 00674 { 00675 for(int i=0;i<l;i++) 00676 alpha_[active_set[i]] = alpha[i]; 00677 } 00678 00679 // juggle everything back 00680 /*{ 00681 for(int i=0;i<l;i++) 00682 while(active_set[i] != i) 00683 swap_index(i,active_set[i]); 00684 // or Q.swap_index(i,active_set[i]); 00685 }*/ 00686 00687 si->upper_bound_p = Cp; 00688 si->upper_bound_n = Cn; 00689 00690 info("\noptimization finished, #iter = %d\n",iter); 00691 00692 delete[] b; 00693 delete[] y; 00694 delete[] alpha; 00695 delete[] alpha_status; 00696 delete[] active_set; 00697 delete[] G; 00698 delete[] G_bar; 00699 }
Here is the call graph for this function:
