BCLS: Bound Constrained Least Squares

Version 0.1

bcsolver.c

Go to the documentation of this file.
00001 /* bcsolver.c
00002    $Revision: 283 $ $Date: 2006-12-17 20:27:56 -0800 (Sun, 17 Dec 2006) $
00003 
00004    ----------------------------------------------------------------------
00005    This file is part of BCLS (Bound-Constrained Least Squares).
00006 
00007    Copyright (C) 2006 Michael P. Friedlander, Department of Computer
00008    Science, University of British Columbia, Canada. All rights
00009    reserved. E-mail: <mpf@cs.ubc.ca>.
00010    
00011    BCLS is free software; you can redistribute it and/or modify it
00012    under the terms of the GNU Lesser General Public License as
00013    published by the Free Software Foundation; either version 2.1 of the
00014    License, or (at your option) any later version.
00015    
00016    BCLS is distributed in the hope that it will be useful, but WITHOUT
00017    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
00018    or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
00019    Public License for more details.
00020    
00021    You should have received a copy of the GNU Lesser General Public
00022    License along with BCLS; if not, write to the Free Software
00023    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00024    USA
00025    ----------------------------------------------------------------------
00026 */
00032 #ifdef __APPLE__
00033   #include <vecLib/vecLib.h>
00034 #else
00035   #include <cblas.h>
00036 #endif
00037 #include <math.h>
00038 #include <assert.h>
00039 #include <string.h>
00040 
00041 #include "bcls.h"
00042 #include "bclib.h"
00043 #include "bccgls.h"
00044 #include "bclsqr.h"
00045 #include "bcsolver.h"
00046 
00047 // =====================================================================
00048 // Private (static) function declarations.
00049 // =====================================================================
00050 
00051 // Compute a Newton step on the free variables, dx(free).
00052 static int
00053 bcls_newton_step( BCLS *ls, int m, int nFree, int ix[], double damp,
00054                   int itnLim, double tol, double dxFree[], double x[],
00055                   double c[], double r[], int *itns, double *rTol );
00056 
00057 // Find the minimizer along the projected search direction.
00058 static int
00059 bcls_proj_search( BCLS *ls, int m, int n, double x[], double dx[],
00060                   double f, double g[], double aBreak[], int iBreak[],
00061                   int ix[], double bl[], double bu[],
00062                   double ex[], double Ad[], double Ae[], int *totHits );
00063 
00064 // Simple backtracking to satisfy an Armijo condition.
00065 static int
00066 bcls_proj_backtrack( BCLS *ls, int m, int n, double x[], double dx[],
00067                      double f, double g[], double bl[], double bu[],
00068                      double s[], int ix[], double As[], int *nSteps );
00069 
00070 
00094 void
00095 bcls_solver( BCLS *ls, int m, int n, double *bNorm,
00096              double x[], double b[], double c[], double bl[], double bu[],
00097              double r[], double g[], double dx[], double dxFree[],
00098              int ix[], double aBreak[], int iBreak[], double anorm[] )
00099 {
00100     int
00101         err = 0,       // Error indicator.
00102         j,             // Misc. counter.
00103         jInf,          // Variable with largest dual infeasibility.
00104         itnLimk,       // No. of allowed itns for kth subproblem.
00105         itnsk = 0,     // No. of iterations performed on kth subproblem.
00106         nFree,         // No. of variables floating between their bounds.
00107         nSteps = 0;    // No. of projected search steps.
00108 
00109     double
00110         optk = 0.0,    // kth     subproblem achieved optimality.
00111         optkm1 = 0.0,  // (k-1)th ...
00112         subTol,        // Required subproblem optimality.
00113         rTol = 0.1,    // Inexact Newton optimality reduction.
00114         dInf,          // Dual infeasibility.
00115         rNorm,         // Norm of the residual.
00116         xNorm,         // Norm of the current solution estimate.
00117         f;             // Objective value = 1/2 * rNorm^2.
00118 
00119     const int scaled_steepest = anorm != NULL;
00120     const int linear     = c != NULL;
00121     const int damped     = ls->damp > 0.0;
00122     const double damp    = ls->damp;
00123     const double damp2   = damp*damp;
00124     const double rTolMin = ls->optTol;
00125 
00126     // Initialize various variables.
00127     ls->itnMaj = 0;    // Reset the major iteration counter.
00128     ls->itnMin = 0;    // Reset the minor iteration counter.
00129 
00130     // Push x into the box.
00131     bcls_mid( n, x, bl, bu );
00132 
00133     // Convergence will be normalized by ||A'b|| + ||c||.
00134     for (j = 0; j < n; j++) ix[j] = j;
00135     bcls_aprod( ls, BCLS_PROD_At, n, ix, g, b );
00136     const double AbNorm = cblas_dnrm2( n, g, 1 );             // ||A'b||
00137     const double cnorm  = (c) ? cblas_dnrm2( n, c, 1 ) : 0.0; // ||  c||
00138     *bNorm = cblas_dnrm2( m, b, 1 );                          // ||  b||
00139 
00140     // -----------------------------------------------------------------
00141     // Start BCLS major iterations.
00142     // -----------------------------------------------------------------
00143     while (1) {
00144 
00145         // Need to include contribution for *all* cols of A.
00146         for (j = 0; j < n; j++) ix[j] = j;
00147 
00148         // Evaluate the (unregularized) residual:
00149         // r(1:m) = b - Ax
00150         bcls_aprod( ls, BCLS_PROD_A, n, ix, x, r );  // r(1:m) =  Ax
00151         cblas_daxpy( m, -1.0, b, 1, r, 1 );          // r(1:m) = -b + Ax
00152         cblas_dscal( m, -1.0, r, 1 );                // r(1:m) =  b - Ax
00153 
00154         // Evaluate the gradient:
00155         // g = A'(Ax - b)  +  damp^2 x = -A'r + damp^2 x.
00156         bcls_aprod( ls, BCLS_PROD_At, n, ix, g, r ); // g =  A'r
00157         cblas_dscal( n, -1.0, g, 1 );                // g = -A'r
00158         if (damped)
00159             cblas_daxpy( n, damp2, x, 1, g, 1 );     // g = g + damp^2 x
00160         if (linear)
00161             cblas_daxpy( n, 1.0, c, 1, g, 1 );       // g = g + c
00162         
00163         // Norm of the current solution and residual.
00164         xNorm = cblas_dnrm2( n, x, 1 );              //  || x ||
00165         rNorm = cblas_dnrm2( m, r, 1 );              //  || r ||
00166             
00167         // Objective value:  f = 1/2 ||r||^2 + 1/2 damp^2 ||x||^2.
00168         f = 0.5 * rNorm * rNorm;
00169         if (damped)
00170             f = f + 0.5 * damp2 * xNorm * xNorm;
00171         if (linear)
00172             f = f + cblas_ddot( n, c, 1, x, 1 );
00173 
00174         // -------------------------------------------------------------
00175         // Evaluate optimality, active set, and exit conditions.
00176         // -------------------------------------------------------------
00177 
00178         // Store free variable indices in ix[:nFree].
00179         nFree = bcls_free_vars( ls, n, ix, x, g, bl, bu );
00180 
00181         // Evaluate optimality of the current iterate x.
00182         bcls_dual_inf( n, x, g, bl, bu, &dInf, &jInf );
00183 
00184         // Test for exit conditions.  Acted on below.
00185         if ( dInf <= ls->optTol * (1.0 + AbNorm + cnorm) ) {
00186             ls->exit      = BCLS_EXIT_CNVGD;
00187             ls->soln_stat = BCLS_SOLN_OPTIM;    // Optimal soln found.
00188         }
00189         else if ( ls->itnMaj >= ls->itnMajLim ) // Too many majors.
00190             ls->exit      = BCLS_EXIT_MAJOR;
00191         else if ( ls->itnMin >= ls->itnMinLim ) // Too many minors.
00192             ls->exit      = BCLS_EXIT_MINOR;
00193         else if ( bcls_callback( ls ) )         // Exit requested by user.
00194             ls->exit      = BCLS_EXIT_CALLBK;
00195 
00196         // -------------------------------------------------------------
00197         // Print out the "kth" iteration log.
00198         // -------------------------------------------------------------
00199         PRINT1(" %5d  %5d  %11.4e  %9.2e  %11.4e  %9.2e %7d %7d",
00200                ls->itnMaj, itnsk, rNorm, xNorm, dInf, optk, nSteps, nFree );
00201         
00202         // Check if rTol was "tight enough" for the last major.
00203         if ( ls->itnMaj <  1 ) {
00204             PRINT1(" %7.0e", rTol);
00205         }
00206         else if ( optk  <= .2 * optkm1 )
00207             ; // Relax. Optimality has made a reasonable reduction.
00208         else if ( rTol  <= rTolMin ) {
00209               // Inexact Newton factor is alrady small.
00210             PRINT1(" %7.0e", rTol);
00211         }
00212         else if ( itnsk <= 1 ) {
00213             rTol /= 10.0;
00214             PRINT1(" %7.0e", rTol);
00215         }
00216 
00217         // Finish the log (eg, print a new-line).
00218         PRINT1("\n");
00219 
00220         // Execute the exit condition (possibly set above).
00221         if ( ls->exit != BCLS_EXIT_UNDEF ) break;
00222 
00223         // =============================================================
00224         // ----  New major iteration starts here ----
00225         // =============================================================
00226         ls->itnMaj++;
00227 
00228         // Zero out the step direction: dx = 0.
00229         bcls_dload( n, 0.0, dx, 1 );
00230 
00231         // -------------------------------------------------------------
00232         // Compute a Newton step on the free variables: dxFree.
00233         // -------------------------------------------------------------
00234         if (nFree) {
00235 
00236             // Find out the number of remaining allowed minors.
00237             itnLimk = ls->itnMinLim - ls->itnMin;
00238 
00239             // Save the previous subproblem optimality value and
00240             // compute the required optimality for this next
00241             // subproblem.
00242             optkm1 = optk;
00243             if (ls->unconstrained)
00244                 subTol = rTol * ls->optTol;
00245             else
00246                 subTol = rTol;// * bcls_vec_dnorm2( nFree, ix, g );
00247             
00248             // Note that ls->dx is used as workspace.
00249             err = bcls_newton_step( ls, m, nFree, ix, damp, itnLimk,
00250                                     subTol, dxFree, x, c, r,
00251                                     &itnsk, &optk );
00252             
00253             // Accumulate the minor iteration count.
00254             ls->itnMin += itnsk;
00255         }
00256         
00257         // -------------------------------------------------------------
00258         // Load steepest descent into the fixed variables: dx(fixed).
00259         // -------------------------------------------------------------
00260         // Scale steepest descent if user provided column scales.
00261         // Note that this needs to be done *after* bcls_newton_step.
00262         cblas_dcopy( n, g, 1, dx, 1 );      //  dx <-  g
00263         cblas_dscal( n, -1.0, dx, 1 );      //  dx <- -g
00264 
00265         if ( scaled_steepest )
00266             for (j = 0; j < n; j++)
00267                 dx[j] /= anorm[j] * anorm[j];
00268 
00269         // Load Newton step on free variables into dx.
00270         bcls_scatter( nFree, ix, dxFree, dx );
00271 
00272         // Projected search along dx.
00273         // Note that ix will be reset.
00274         if (ls->proj_search == BCLS_PROJ_SEARCH_EXACT)
00275             err = bcls_proj_search( ls, m, n, x, dx, f, g, aBreak, iBreak,
00276                                     ix, bl, bu,
00277                                     ls->wrk_u, ls->wrk_v, ls->wrk_w,
00278                                     &nSteps );
00279         else
00280             err = bcls_proj_backtrack( ls, m, n, x, dx, f, g, bl, bu,
00281                                        ls->wrk_v, ix, ls->wrk_u,
00282                                        &nSteps );
00283 
00284         // -------------------------------------------------------------
00285         // Check for linesearch failures.
00286         // -------------------------------------------------------------
00287 
00288         // Exit if linesearch found a direction of unbounded descent.
00289         if (err == BCLS_EXIT_UNBND) {
00290             ls->exit = err;
00291             break;
00292         }
00293 
00294         // Generic linesearch failure.  Further tests needed.
00295         else if (err == BCLS_EXIT_LFAIL) {
00296 
00297             // If using approx linesearch, revert to exact linesearch.
00298             if (ls->proj_search == BCLS_PROJ_SEARCH_APPROX) {
00299                 ls->proj_search =  BCLS_PROJ_SEARCH_EXACT;
00300                 PRINT1(" W: Too many approx backtrack iterations."
00301                        " Switching to exact linesearch.\n");
00302             }
00303             // Already using exact linesearch.  Exit.
00304             else {
00305                 ls->exit = err;
00306                 break;
00307             }
00308         }
00309 
00310     } // end of while (1) -- the main iteration loop.
00311 
00312     // Collect various solution statistics.
00313     ls->soln_obj   = f;
00314     ls->soln_rNorm = rNorm;
00315     ls->soln_jInf  = jInf;
00316     ls->soln_dInf  = dInf;
00317 
00318     // Exit.
00319     return;
00320 }
00321 
00322 
00406 static int
00407 bcls_newton_step( BCLS *ls, int m, int nFree, int ix[], double damp,
00408                   int itnLim, double tol, double dxFree[], double x[],
00409                   double c[], double r[], int *itns, double *opt )
00410 {
00411     int err;
00412     const int preconditioning = ls->Usolve != NULL;
00413     double *dx = ls->dx;  // Used as workspace.
00414 
00415     // -----------------------------------------------------------------
00416     // Initialize the Aprod routine, and the Usolve routine if given.
00417     // -----------------------------------------------------------------
00418     bcls_aprod ( ls, BCLS_PROD_INIT, nFree, ix, NULL, NULL );
00419     if (preconditioning)
00420         bcls_usolve( ls, BCLS_PRECON_INIT, nFree, ix, NULL, NULL );
00421 
00422     // -----------------------------------------------------------------
00423     // Print dividing line in the subproblem log.
00424     // -----------------------------------------------------------------
00425     if (ls->minor_file) {
00426         fprintf(ls->minor_file,
00427                 "\n"
00428                 "------------------------- "
00429                 "BCLS Major Iteration: %6d "
00430                 "------------------------- "
00431                 "\n\n", ls->itnMaj);
00432         fflush(ls->minor_file);
00433     }
00434    
00435     // -----------------------------------------------------------------
00436     // Call one of the subproblem solvers.
00437     // -----------------------------------------------------------------
00438     if (ls->newton_step == BCLS_NEWTON_STEP_LSQR)
00439 
00440         err = bcls_newton_step_lsqr( ls, m, nFree, ix, damp,
00441                                      itnLim, tol, dxFree, x,
00442                                      c, r, itns, opt );
00443     
00444     else
00445 
00446         err = bcls_newton_step_cgls( ls, m, nFree, ix, damp,
00447                                      itnLim, tol, dxFree, x,
00448                                      c, r, itns, opt );
00449 
00450     // -----------------------------------------------------------------
00451     // If the user is preconditioning, then dxFree currently lives in
00452     // the preconditioned space.  We need one more call to Usolve to
00453     // recover dxFree in the original space.  Take this opportunity to
00454     // call Usolve in its termination mode.
00455     // -----------------------------------------------------------------
00456     if (preconditioning) {
00457 
00458         bcls_usolve( ls, BCLS_PRECON_U, nFree, ix, dx, dxFree );
00459         cblas_dcopy( nFree, dx, 1, dxFree, 1 );  //  dxFree <- dx
00460 
00461         bcls_usolve( ls, BCLS_PRECON_TERM, nFree, ix, NULL, NULL );
00462     }
00463 
00464     // Call Aprod in its termination mode.
00465     bcls_aprod ( ls, BCLS_PROD_TERM, nFree, ix, NULL, NULL );
00466 
00467     // Exit.
00468     return err;
00469 }
00470 
00490 static int
00491 bcls_proj_search( BCLS *ls, int m, int n, double x[], double dx[],
00492                   double f, double g[], double aBreak[], int iBreak[],
00493                   int ix[], double bl[], double bu[],
00494                   double ex[], double Ad[], double Ae[], int *totHits )
00495 {
00496     int j, k;
00497     int err     = 0;  // Error flag.
00498     int nBrkPts = 0;  // No. of breakpoints encountered.
00499     int nFree   = 0;  // No. of variables that are free to move.
00500     int xlower;
00501     int xupper;
00502     int xfree;
00503     int stepstat = 0; // Flag if step is
00504     const int STEP_UND = -1; //   undefined
00505     const int STEP_LHS =  0; //   on LHS
00506     const int STEP_MID =  1; //   to interior of interval
00507     const int STEP_RHS =  2; //   to RHS
00508     const int STEP_INF =  3; //   unbounded descent.
00509     int gdpos;        // Flag if gd is "positive".
00510     int gdzero;       // Flag if gd is "zero".
00511     int dHdNonNeg;    // Flag if dHd is "nonnegative".
00512     int dHdNonPos;    // Flag if dHd is "nonpositive".
00513     int jHits;        // No. of vars that hit ther bound at break point j.
00514 
00515     double lhs;       // The start (LHS) of the current interval.
00516     double rhs;       // The end (RHS) of the current interval.
00517     double step;      // Fraction of the current search interval.
00518     double maxStep;   // The width of the current search interval.
00519     double q;         // The piecewise-quadratic objective function.
00520     double gd;        // Projected gradient of q along current arc.
00521     double ge;        // Projected gradient of q along last arc.
00522     double dHd;       // Projected Hessian  of q along current arc.
00523     double dHe;       // Inner product of Hessian along d and update e.
00524 
00525     const int damped = ls->damp > 0.0;
00526     const double damp2 = ls->damp * ls->damp;
00527     const double eps = ls->eps;
00528     const double epsx = ls->epsx;
00529     const double epsfixed = ls->epsfixed;
00530     const double BigNum = ls->BigNum;
00531 
00532     PRINT6("\n%3s %10s %10s %10s %10s %10s\n",
00533            "Var","Lower","x","Upper","dx","Break");
00534 
00535     // -----------------------------------------------------------------
00536     // Find out which variables are free to move.
00537     // -----------------------------------------------------------------
00538     for (j = 0; j < n; j++) {
00539 
00540         if ( bu[j] - bl[j] < epsfixed ) {
00541             // This variable is fixed.  Zero out its search direction.
00542             dx[j] = 0.0;
00543             continue;
00544         }
00545         
00546         xupper = bu[j] - x[j] <= epsx; // Variable on upper bound.
00547         xlower = x[j] - bl[j] <= epsx; // Variable on lower bound.
00548         xfree  = !xlower && !xupper;   // Variable is floating btwn bnds.
00549         
00550         assert( xupper || xlower || xfree );
00551 
00552         if ( xfree ) {
00553             if ( fabs(dx[j]) > eps ) {
00554                 iBreak[nFree] = j;
00555                 nFree++;
00556             }
00557             else
00558                 dx[j] = 0.0;
00559         }
00560 
00561         else if ( xupper ) {
00562             if ( dx[j] < - eps ) {
00563                 iBreak[nFree] = j;
00564                 nFree++;
00565             }
00566             else
00567                 dx[j] = 0.0;
00568         }
00569 
00570         else { // if ( xlower )
00571             if ( dx[j] > eps ) {
00572                 iBreak[nFree] = j;
00573                 nFree++;
00574             }
00575             else
00576                 dx[j] = 0.0;
00577         }
00578     }
00579     
00580     // -----------------------------------------------------------------
00581     // Compute the breakpoints of the free variables.
00582     //
00583     // Each breakpoint measures the maximum allowable step that will
00584     // take the corresponding free variable to its boundary.
00585     // -----------------------------------------------------------------
00586     for (j = 0; j < nFree; j++) {
00587         k = iBreak[j];
00588         if ( dx[k] > 0.0 )
00589             aBreak[j] = ( bu[k] - x[k] ) / dx[k];
00590         else
00591             aBreak[j] = ( bl[k] - x[k] ) / dx[k];
00592     }
00593 
00594     // -----------------------------------------------------------------
00595     // Diagnostic printing.
00596     // -----------------------------------------------------------------
00597     if (ls->print_level >= 6) {
00598         k = 0;
00599         for (j = 0; j < n; j++) {
00600             PRINT6("%3d %10.3e %10.3e %10.3e %10.3e ",
00601                    j, bl[j], x[j], bu[j], dx[j]);
00602             if ( j == iBreak[k] ) {
00603                 PRINT6("%10.3e\n",aBreak[k]);
00604                 k++;
00605             }
00606             else
00607                 PRINT6("%10s\n","--");
00608         }
00609         PRINT6("\n");
00610     }
00611 
00612     // Exit if no variables are free to move.
00613     if ( !nFree ) return 0;
00614 
00615     // No. of breakpoints remaining is initially the number of free vars.
00616     nBrkPts = nFree;
00617 
00618     // Build the breakpoints into a heap.
00619     bcls_heap_build( nBrkPts, aBreak, iBreak );
00620 
00621     // -----------------------------------------------------------------
00622     // Projected gradient of q along dx.
00623     // -----------------------------------------------------------------
00624     gd = 0.0;
00625     for (j = 0; j < nFree; j++) {
00626         k = iBreak[j];
00627         gd += g[k] * dx[k];
00628     }
00629 
00630     // Quick exit if the directional derivative is positive.
00631     if ( gd > epsx ) return 0;
00632 
00633     // -----------------------------------------------------------------
00634     // Projected Hessian of q along dx:
00635     // dHd = dx' A' A dx = (A dx)' (A dx) = Ad' Ad.
00636     // -----------------------------------------------------------------
00637     bcls_aprod( ls, BCLS_PROD_A, nFree, iBreak, dx, Ad );
00638     dHd = cblas_ddot( m, Ad, 1, Ad, 1 );
00639     if (damped)
00640         dHd += damp2 * cblas_ddot( n, dx, 1, dx, 1 );
00641 
00642     // Initialization.
00643     q        = 0.0; // Value of piecewise quadratic.
00644     rhs      = 0.0; // This initializes the very first breakpoint.
00645     *totHits = 0;   // Initially no hits.
00646 
00647     // -----------------------------------------------------------------
00648     // Explore each piecewise quadratic intervals until non remain.
00649     // -----------------------------------------------------------------
00650     while (nBrkPts > 0) {
00651 
00652         // Set the status of the current step to "undefined".
00653         stepstat = STEP_UND;
00654 
00655         // Grab the boundaries of the current search interval.
00656         lhs     = rhs;
00657         rhs     = aBreak[0];  // The next breakpoint (top of the heap).
00658         maxStep = rhs - lhs;  // The maximum step.
00659         k       = iBreak[0];  // Var corresponding to the current break.
00660 
00661         // The smallest breakpoint has been sampled.  Move it to the
00662         // top of the heap and restore the heap property.
00663         nBrkPts = bcls_heap_del_min( nBrkPts, aBreak, iBreak );
00664 
00665         // -------------------------------------------------------------
00666         // Test the various exit or continuation conditions.
00667         // -------------------------------------------------------------
00668         gdpos     = gd       >    epsx;
00669         gdzero    = fabs(gd) <=   epsx;
00670         dHdNonNeg = dHd      >= - epsx;
00671         dHdNonPos = dHd      <    epsx;
00672 
00673         if ( gdpos || ( gdzero && dHdNonNeg) ) {
00674             // ---------------------------------------------------------
00675             // Minimimum is on LHS. (Captures the  dx = 0  case.)  Done.
00676             // ---------------------------------------------------------
00677             step = 0.0;
00678             stepstat = STEP_LHS;
00679             break;
00680         }       
00681 
00682         // If we got this far, gd is zero or negative.
00683 
00684         else if ( dHdNonPos ) {
00685             // ---------------------------------------------------------
00686             // Found dir of neg curvature.  Minimizer is either on RHS,
00687             // or at -infinity.  Check.
00688             // ---------------------------------------------------------
00689             if ( maxStep >= BigNum ) {
00690                 stepstat = STEP_INF;
00691                 break;
00692             }
00693             else {
00694                 step = maxStep;
00695                 stepstat = STEP_RHS;
00696             }
00697         }
00698         else {
00699             // ---------------------------------------------------------
00700             // Minimum is on RHS or inside the interval.  Check.
00701             // ---------------------------------------------------------
00702             step = - gd / dHd;
00703             if ( step >= BigNum  &&  step < maxStep ) {
00704                 // Step is unbounded.  Exit with unbounded flag.
00705                 stepstat = STEP_INF;
00706                 break;
00707             }
00708             else if ( step <  maxStep ) {
00709                 // Minimum found in the interval's interior.  Done.
00710                 stepstat = STEP_MID;
00711                 break;
00712             }
00713             else if ( step >= maxStep ) {
00714                 // Minimum is on RHS.  Truncate the step and continue.
00715                 step = maxStep;
00716                 stepstat = STEP_RHS;
00717             }
00718         }
00719         // -------------------------------------------------------------
00720         // Take the step.
00721         // -------------------------------------------------------------
00722         cblas_daxpy( n, step, dx, 1, x, 1 );
00723             
00724         // -------------------------------------------------------------
00725         // Other variables may share (essentially) the same breakpoint.
00726         //
00727         // Loop over these (and the current) variables, pushing them
00728         // onto their bounds, zeroing out the corresponding search
00729         // direction.
00730         // -------------------------------------------------------------
00731         jHits = 0; // No. of vars that hit their bound at this breakpoint.
00732         while(1) {
00733             
00734             // There may have been roundoff error.  Push the var that
00735             // hit its breakpoint exactly onto its bound.
00736             if (dx[k] < 0.0)
00737                 x[k] = bl[k];
00738             else
00739                 x[k] = bu[k];
00740 
00741             // Record the index and step of the var that just hit its bound.
00742             ex[k] = dx[k];
00743             ix[jHits] = k;
00744             jHits++;
00745 
00746             // That var no longer contributes.  Zero out its step.
00747             dx[k] = 0.0;
00748             
00749             // Print something out.
00750             PRINT6("%3d %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
00751                    k,bl[k],x[k],bu[k],dx[k],gd,dHd);
00752 
00753             // No more breakpoints left to explore.
00754             if (nBrkPts == 0)
00755                 break;
00756 
00757             // Next breakpoint is significantly different.
00758             else if (aBreak[0] > rhs + 10*eps)
00759                 break;
00760 
00761             // Next breakpoint is essentially the same.
00762             else {
00763                 k = iBreak[0];
00764                 nBrkPts = bcls_heap_del_min( nBrkPts, aBreak, iBreak );
00765             }
00766         }
00767         
00768         // Update the aggregate no. of breakpoints encountered so far.
00769         *totHits += jHits; 
00770 
00771         // -------------------------------------------------------------
00772         // Update the function value.  (Used for printing only.)
00773         // -------------------------------------------------------------
00774         q  = q + step * (gd + 0.5 * step * dHd);
00775 
00776         // -------------------------------------------------------------
00777         // Update the directional derivative and Hessian for the next
00778         // search direction.
00779         // -------------------------------------------------------------
00780         
00781         // Inner products: g'e.
00782         ge = 0.0;
00783         for (j = 0; j < jHits; j++) {
00784             k = ix[j];
00785             ge += g[k]*ex[k];
00786         }
00787 
00788         // Inner product: d' H e.
00789         bcls_aprod( ls, BCLS_PROD_A, jHits, ix, ex, Ae ); //  Ae = A * ex
00790         dHe = cblas_ddot( m, Ad, 1, Ae, 1 );              // dHe = d' A' A e
00791         if (damped) {
00792             for (j = 0; j < jHits; j++) {
00793                 k = ix[j];
00794                 dHe += damp2 * dx[k] * ex[k];
00795             }
00796         }
00797 
00798         // Update projected gradient.
00799         gd  = gd - ge + step * ( dHd - dHe );
00800         
00801         // Update projected Hessian.
00802         cblas_daxpy( m, -1.0, Ae, 1, Ad, 1 );     //  Ad <- Ad - Ae
00803         dHd = cblas_ddot( m, Ad, 1, Ad, 1 );      // dHd <- d' A' A d
00804         if (damped)
00805             dHd = dHd + damp2 * cblas_ddot( n, dx, 1, dx, 1 );
00806 
00807     } // end of while (nBrkPs > 0)
00808 
00809     if ( stepstat == STEP_MID )
00810         // Step into the middle of an interval and exit.
00811         cblas_daxpy( n, step, dx, 1, x, 1 );
00812 
00813     else if ( stepstat == STEP_INF )
00814         err = BCLS_EXIT_UNBND;
00815 
00816     // Return the status of the last step.
00817     return err;
00818 }
00819 
00845 static int 
00846 bcls_proj_backtrack( BCLS *ls, int m, int n, double x[], double dx[],
00847                      double f, double g[], double bl[], double bu[],
00848                      double s[], int ix[], double As[], int *nSteps )
00849 {
00850     int j;
00851     int armijo = 0;   // Flag: Armijo condition satisfied.
00852     int unbounded;    // Flag: Step may be unbounded.
00853     int descent;      // Flag: Descent direction.
00854     int negCurv;      // Flag: Nonpositive curvature.
00855     int xupper;       // Flag: Variable on upper bound.
00856     int xlower;       // Flag: Variable on lower bound.
00857     int xdecr;        // Flag: Variable is decreasing.
00858     int xincr;        // Flag: Variable is increasing.
00859     int nBreakPts;    // No. of break points.
00860     int jBreakMin;    // Index of variable with the nearest breakpoint.
00861     double Break;     // Step length to the break point for variable j.
00862     double BreakMin;  // Step length to the nearest breakpoint.
00863     double q;         // The quadratic objective function.
00864     double gs;        // Projected gradient along current search direction.
00865     double sHs;       // Projected Hessian  along current search direction.
00866     double step = 1;  // Current step length.
00867 
00868     const int    backlimit = ls->backtrack_limit;
00869     const int    damped    = ls->damp > 0.0;
00870     const double damp2     = ls->damp * ls->damp;
00871     const double mu        = ls->mu;
00872     const double backtrack = ls->backtrack;
00873     const double BigNum    = ls->BigNum;
00874     const double eps       = ls->eps;
00875     const double epsx      = ls->epsx;
00876 
00877     *nSteps = 0;      // No. of backtracking steps.
00878 
00879     // -----------------------------------------------------------------
00880     // Need to include contribution for *all* columns of A.
00881     // -----------------------------------------------------------------
00882     for (j = 0; j < n; j++) ix[j] = j;
00883 
00884     // -----------------------------------------------------------------
00885     // Find first breakpoint, and store as BreakMin.
00886     // -----------------------------------------------------------------
00887     nBreakPts = 0;
00888     BreakMin  = BigNum;
00889     for (j = 0; j < n; j++) {
00890         xupper = bu[j] - x[j] <= epsx; // Variable on upper bound.
00891         xlower = x[j] - bl[j] <= epsx; // Variable on lower bound.
00892         xdecr  = dx[j] < -eps;         // Variable is decreasing.
00893         xincr  = dx[j] >  eps;         // Variable is increasing.
00894 
00895         if ( !xlower  &&  xdecr ) {
00896             nBreakPts++;
00897             Break = ( bl[j] - x[j] ) / dx[j];
00898             if ( Break < BreakMin ) {
00899                 BreakMin  = Break;
00900                 jBreakMin = j;
00901             }
00902         }
00903         else if ( !xupper  &&  xincr ) {
00904             nBreakPts++;
00905             Break = ( bu[j] - x[j] ) / dx[j];
00906             if ( Break < BreakMin ) {
00907                 BreakMin  = Break;
00908                 jBreakMin = j;
00909             }
00910         }
00911     }
00912     
00913     if ( nBreakPts == 0 )
00914         BreakMin = 0.0;
00915 
00916     PRINT3("I: Proj Backtrack: Steplength to the nearest bound: %9.2e\n",
00917            BreakMin);
00918     
00919     // -----------------------------------------------------------------
00920     // Backtracking loop.
00921     // -----------------------------------------------------------------
00922     while ( !armijo  &&  step > BreakMin ) {
00923 
00924         if (*nSteps >= backlimit)
00925             return BCLS_EXIT_LFAIL;
00926         
00927         // -------------------------------------------------------------
00928         // Store the projected step into s:  s = P[x + step*dx] - x.
00929         // -------------------------------------------------------------
00930         bcls_project_step( n, s, step, x, dx, bl, bu );
00931 
00932         // -------------------------------------------------------------
00933         // Compute q(s) = g's + 0.5 s' H s.
00934         // -------------------------------------------------------------
00935         bcls_aprod( ls, BCLS_PROD_A, n, ix, s, As ); //  As = A * s
00936         sHs = cblas_ddot( m, As, 1, As, 1 );         // sHs = (A s)'(A s)
00937         if (damped)
00938             sHs = sHs + damp2 * cblas_ddot( n, s, 1, s, 1 );
00939         gs  = cblas_ddot( n, g, 1, s, 1 );    //  gs = g's
00940         q   = gs + 0.5 * sHs;                 //   q = g's + 0.5 s' H s
00941         
00942         // -------------------------------------------------------------
00943         // If sHs <= 0 and gs < 0, there may be a direction of
00944         // unbounded descent.  Check.
00945         // -------------------------------------------------------------
00946         descent  =  gs  <= -epsx;
00947         negCurv  =  sHs <=  epsx;
00948 
00949         if ( descent  &&  negCurv ) {
00950 
00951             unbounded = 1;
00952             for (j = 0; j < n; j++) {
00953                 if (s[j] > eps)
00954                     unbounded = bu[j] >=   BigNum;
00955                 
00956                 else if (s[j] < eps)
00957                     unbounded = bl[j] <= - BigNum;
00958                 
00959                 if (!unbounded)
00960                     break;
00961             }
00962             
00963             if (unbounded)
00964                 return BCLS_EXIT_UNBND;
00965         }
00966 
00967         // -------------------------------------------------------------
00968         // Check if the sufficient decrease criteria was satisfied.
00969         // Otherwise, reduce the step and continue.
00970         // -------------------------------------------------------------
00971         PRINT3("I: Proj Backtrack: step = %9.2e  q = %10.2e  gs = %10.2e\n",
00972                step, q, gs);
00973 
00974         armijo = ( q <= mu * gs );
00975 
00976         if (!armijo) {
00977             step = backtrack * step;
00978             *nSteps += 1;
00979         }
00980     }
00981 
00982     // -----------------------------------------------------------------
00983     // Force a variable onto its nearest bound: If the armijo
00984     // condition is satisfied by a truncated step that falls short of
00985     // the nearest breakpoint, then we need to be a bit more
00986     // aggressive and continue the step until the nearest boundary.
00987     //
00988     // 17 Dec 06: Removed this step.  Why was I doing it?!
00989     // -----------------------------------------------------------------
00990 #ifdef UNDEF
00991     if ( step < 1.0  &&  step < BreakMin ) {
00992         step = BreakMin;
00993         PRINT3("I: Proj Backtrack: Force index %4i onto bound\n",
00994                jBreakMin );
00995     }
00996 #endif
00997 
00998     // Compute the final projected step.
00999     bcls_project_step( n, s, step, x, dx, bl, bu );
01000     
01001     // Take the step: x = x + s.
01002     cblas_daxpy( n, 1.0, s, 1, x, 1 );
01003     
01004     // Push x onto bounds that are very close.
01005     for (j = 0; j < n; j++) {
01006         if      (x[j] < bl[j] + 10*eps)
01007                  x[j] = bl[j];
01008         else if (x[j] > bu[j] - 10*eps)
01009                  x[j] = bu[j];
01010     }
01011 
01012     return 0;
01013 }

Generated on Sun Mar 4 22:50:03 2007 by Doxygen 1.5.1