snobfit.minq

38 """ 39 Minimizes an affine quadratic form subject to simple bounds, 40 using coordinate searches and reduced subspace minimizations 41 using LDL^T factorization updates 42 min fval = gamma + c^T x + 0.5 x^T G x } 43 s.t. x in [xu,xo] 44 45 where G is a symmetric (n x n) matrix, not necessarily definite 46 (if G is indefinite, only a local minimum is found) 47 48 If G is sparse, it is assumed that the ordering is such that 49 a sparse modified Cholesky factorization is feasible 50 """ 51 convex = 0 52 53 # The maximal number of iterative refinement steps 54 nitrefmax = 3 55 56 # The perturbation in last two digits 57 neps = 100*eps 58 59 # The number of subspace steps 60 nsub = 0 61 62 # The allow variables to be freed in csearch? 63 unfix = 1 64 65 # No iterative refinement steps so far 66 nitref = 0 67 68 # Improvement expected 69 improvement = 1 70 71 # Best function value 72 fval = inf 73 74 nfree = 0 75 nfree_old = -1 76

77 - def __init__(self, gamma, c, G, xu, xo, prt=0, x0=None):

78 """ 79 Inputs: 80 ------- 81 gamma a constant. 82 c a colomn vector. 83 G a symmetric (n x n) matrix. 84 xu lower bound 85 xo upper bound 86 87 Optional Inputs: 88 ---------------- 89 prt print level. 90 xx initial guess (optional) 91 """ 92 self.gamma = gamma 93 self.c = c 94 self.G = numpy.mat(G) 95 self.xu = xu 96 self.xo = xo 97 self.prt = prt 98 self.xTry = x0 99 self.n = self.G.shape[0] 100 101 # Initialization 102 if prt>0: 103 self.print_init()

104 105

106 - def print_minq(self):

107 print '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 108 print '!!!!! Minq !!!!!' 109 print '!!!!! incomplete minimization !!!!!' 110 print '!!!!! too many iterations !!!!!' 111 print '!!!!! increase maxit !!!!!' 112 print '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'

113 114

115 - def print_init(self):

116 print '====================================================' 117 print '====================== Minq ======================' 118 print '===================================================='

119 120

121 - def setSearchInit(self):

122 123 # Initialize trial point xx, function value fval and gradient g 124 if self.xTry is None: 125 # Cold start with absolutely smallest feasible point 126 self.xTry = numpy.empty( self.n, 'd' ) 127 128 # Force starting point into the box 129 self.xTry = numpy.maximum( self.xu, numpy.minimum(self.xTry,self.xo) ) 130 131 # Maximal number of iterations 132 # This limits the work to about 1+4*maxit/n matrix multiplies 133 # Usually at most 2*n iterations are needed for convergence 134 self.maxit = 3 * self.n 135 136 # Regularization for low rank problems, perturbation in last two digits 137 for i in xrange( self.n ): 138 self.G[i,i] *= (1+self.neps) 139 140 # Initialize LDL^T factorization of G_KK 141 self.K = numpy.zeros( self.n ) # initially no rows in factorization 142 self.L = numpy.eye(self.n) 143 self.d = numpy.ones( self.n ) 144 145 # Dummy initialization of indicator of free variables 146 # will become correct after first coordinate search 147 self.free = [0] * self.n

148 149

150 - def _quadMin(self, low, high, b, a):

151 """ 152 Get minimizer x in [low, high] for a univariate quadratic 153 154 q(x) = 0.5*a*x^2 + b*x 155 156 Output: 157 ------ 158 x: minimizer x in [low, high] 159 lba: lower bound active 160 uba: upper bound active 161 162 ier: 0 (finite minimizer) 163 1 (unbounded minimum) 164 """ 165 lba = 0 # lower bound is active or not? 166 uba = 0 # upper bound is active or not? 167 168 # Determine unboundedness 169 ier=0 170 if (low == -inf) and ( a<0 or (a==0 and b>0) ): 171 ier=1; lba=1 172 173 if (high == inf) and ( a<0 or (a==0 and b<0) ): 174 ier=1; uba=1 175 176 if ier: 177 return (nan, lba, uba, ier ) 178 179 # Determine activity 180 if a==0 and b==0: 181 x=0 182 183 elif a<=0: 184 # Concave case minimal at a bound 185 if low == -inf: lba = 0 186 elif high == inf: lba = 1 187 else: lba = (2*b+(low+high)*a>0) 188 uba = not lba 189 190 else: 191 x = -float(b)/float(a) # Unconstrained optimal step 192 lba = (x <= low ) # Lower bound is active? 193 uba = (x >= high) # Upper bound is active? 194 195 if lba: x = low 196 if uba: x = high 197 198 return (x, lba, uba, ier )

199 200

201 - def _Info(self):

202 203 self.info = 0 204 if not self.improvement: # Good termination 205 if self.prt: 206 print 'Terminate: no improvement in coordinate search' 207 208 elif self.nitref > self.nitrefmax: # Good termination 209 if self.prt: print 'Terminate print_minq: nitref>nitrefmax' 210 211 elif self.nitref > 0 and \ 212 self.nfree_old == self.nfree and \ 213 self.fvalnew >= self.fval: # Good termination 214 215 if self.prt: 216 print 'Terminate:nitref>0&nfree_old==nfree&fvalnew>=fval' 217 218 else: 219 self.info = -1 220 221 if self.nitref==0 and self.nsub == self.maxit: 222 if self.prt: self.print_minq() 223 else: print 'iteration limit exceeded' 224 self.info = 99

225 226 227 # -----------------------------------------------------------------------

228 - def ldlrk1(self, L, d, alpha, u):

229 """ 230 Computes LDL^T factorization for LDL^T + alpha * u * u^T 231 if alpha>=0 or if new factorization is definite(both signalled by p=[]) 232 otherwise, the original L,d and 233 a direction p of null or negative curvature are returned 234 235 d contains diag(D) and is assumed positive 236 237 Warning: does not work for dimension 0 238 """ 239 L = numpy.asarray(L) 240 p = numpy.array([]) 241 if alpha==0: 242 return 243 n = u.size 244 245 # Save old factorization 246 L0 = L.copy() 247 d0 = d.copy() 248 249 # Update 250 for k in find(u!=0): 251 252 delta = d[k] + alpha * u[k]**2 253 if alpha<0 and delta <= self.neps: 254 # Update not definite 255 p = numpy.zeros(n, 'd') 256 p[k] = 1 257 p[0:k+1] = toRow ( mldiv( numpy.asmatrix(L[0:k+1,0:k+1].T), 258 toCol( numpy.asarray( p[0:k+1] ) ) 259 ) ) 260 261 # Restore original factorization 262 L = L0 263 d = d0 264 return (L, d, p) 265 266 q = d[k]/delta 267 d[k] = delta 268 269 # ---------------------------------------------- 270 ind = range(k,n) 271 c = L[ind,k]*u[k] 272 L[ind,k] = L[ind,k]*q + (alpha*u[k]/delta)*u[ind] 273 u[ind] -= c 274 275 alpha *= q 276 if alpha==0: 277 break 278 279 return (L, d, p)

280 281

282 - def ldlup(self, L, d, j, g):

283 """ 284 Updates LDL^T factorization when a unit j-th row and column 285 are replaced by column g 286 287 If the new matrix is definite (signalled by p=[]); 288 otherwise, the original L,d and 289 a direction p of null or negative curvature are returned 290 291 d contains diag(D) and is assumed positive 292 Note that g must have zeros in other unit rows!!! 293 """ 294 p = numpy.array([]) 295 n = d.size 296 if j==0: 297 delta = g[j] 298 if delta <= self.neps: 299 p = numpy.zeros( n, 'd') 300 p[0] = 1.0 301 return (L,d,p) 302 303 d[j] = delta 304 return (L, d, p) 305 306 # Now j>1, K nonempty 307 LII = L[0:j,0:j] 308 tmat = g[0:j] 309 u = mldiv(LII, tmat ) 310 v = u / d[0:j] 311 delta = g[j] - numpy.dot(u, v) 312 if delta <= self.neps: 313 pLII = mldiv(LII.T,v) 314 p = numpy.zeros( n, 'd') 315 for i in xrange( len(pLII) ): 316 p[i] = pLII[i] 317 p[j] = -1 318 return (L, d, p) 319 320 LKI = L[j+1:n,0:j] 321 LKI_u = LKI*numpy.mat( toCol(u) ) 322 323 w = ( toCol(g[j+1:n]) - LKI_u )/delta 324 325 (LKK,d[j+1:n],q) = self.ldlrk1( L[j+1:n,j+1:n], 326 d[j+1:n], 327 -delta, 328 toRow(w) 329 ) 330 331 if len(q)==0: 332 # Work around expensive sparse L(K,K)=LKK 333 L1 = L[0:j,] 334 L2 = numpy.mat( numpy.zeros(n) ) 335 L3 = numpy.mat(L[j,j+1:n]) 336 L_1 = numpy.mat( 1 ) 337 Lw = numpy.mat( w ) 338 Lv = numpy.mat( v ) 339 L4 = numpy.mat( numpy.zeros( (n-j-1,1) ) ) 340 L = numpy.bmat( 'L1; Lv,L_1,L3; LKI,Lw,LKK' ) 341 d[j]=delta 342 343 else: # Work around expensive sparse L(K,K)=LKK 344 L1 = L[0:j+1,:] 345 L2 = L[j+1:n,j] 346 L = numpy.bmat( 'L1; LKI,L2,LKK' ) 347 pi = crdot( w, q ) 348 LKI_q = LKI.T * numpy.mat( toCol(q) ) 349 pi_v = pi*toCol(v) 350 351 p1 = mldiv( LII.T, pi_v-LKI_q) 352 353 p = numpy.vstack( ( toCol(p1), numpy.array(-pi), toCol(q) ) ) 354 355 return (L, d, p)

356 357 358 359 # ---------------------------------------------------------------------

360 - def ldldown(self, L, d, j):

361 """ 362 Downdates LDL^T factorization 363 when j-th row and column are replaced by j-th unit vector 364 365 d contains diag(D) and is assumed positive 366 """ 367 L = numpy.asmatrix(L) 368 n = d.size 369 370 w = toRow( L[j+1:n,j] ) 371 372 if j<n: 373 [ LKK, d[j+1:n], q ] = self.ldlrk1( L[j+1:n,j+1:n], 374 d[j+1:n], 375 d[j], 376 w 377 ) 378 # Work around expensive sparse L[K,K]=LKK 379 L1 = L[0:j,] 380 L2 = numpy.mat( numpy.zeros(n) ) 381 L3 = L[j+1:n, 0:j] 382 L4 = numpy.mat( numpy.zeros( (n-j-1,1) ) ) 383 L = numpy.bmat( 'L1; L2; L3,L4,LKK' ) 384 L[j,j] = 1 385 386 else: 387 for i in xrange(n-1): 388 L[n,i] = 0.0 389 390 d[j]=1 391 392 return (L,d)

393 394

395 - def coordinateSearch(self):

396 """ 397 Coordinate search 398 """ 399 # The number of consecutive free steps 400 count = 0 401 402 # Current coordinate searched 403 k = -1 404 while 1: 405 406 while count <= self.n: 407 # Find next free index (or next index if unfix) 408 count += 1 409 if k == self.n-1: k=0 410 k += 1 #? 411 412 if self.free[k] or self.unfix: 413 break 414 415 if count > self.n: 416 # Complete sweep performed without fixing a new active bound 417 break 418 419 q = self.G[:,k] 420 alpha_u = self.xu[k] - self.x[k] 421 alpha_o = self.xo[k] - self.x[k] # bounds on step 422 423 # Find step size 424 [alpha, lba, uba, self.info] = \ 425 self._quadMin( alpha_u, alpha_o, self.gradient[k], q[k] ) 426 427 if self.info: 428 x = numpy.zeros( self.n ) 429 if lba: x[k] = -1 430 else: x[k] = 1 431 return 432 433 # New x 434 xnew = self.x[k] + alpha 435 if self.prt and self.nitref>0: 436 print xnew,alpha 437 438 if lba or xnew <= self.xu[k]: # Lower bound active 439 if alpha_u !=0 : 440 self.x[k] = self.xu[k] 441 self.gradient += alpha_u * q 442 count=0 443 self.free[k] = 0 444 445 elif uba or xnew >= self.xo[k]: # Upper bound active 446 if alpha_o != 0: 447 self.x[k] = self.xo[k] 448 self.gradient += alpha_o * q 449 count = 0 450 self.free[k] = 0 451 452 else: # No bound active 453 if alpha != 0: 454 if self.prt>1 and not self.free[k]: 455 unfixstep=[ x[k], alpha ] 456 457 self.x[k] = xnew 458 self.gradient += alpha * q 459 self.free[k] = 1

460 # ============== End of coordinate search ================ 461 462

463 - def subspace(self):

464 """ 465 Take a subspace step 466 """ 467 self.nsub += 1 468 469 # Downdate factorization 470 for j in find( self.free < self.K ): # list of newly active indices 471 [self.L, self.d] = self.ldldown(self.L, self.d, j) 472 self.K[j] = 0 473 474 # Update factorization or find indefinite search direchtion 475 definite = 1 476 for j in find( self.free > self.K ): # List of newly freed indices 477 p = numpy.zeros( self.n ) 478 ind = find( self.K > 0 ) 479 p[ind] = self.G[ind,j] 480 p[j] = self.G[j,j] 481 [self.L, self.d, pp] = self.ldlup(self.L, self.d, j, p) 482 483 definite = len(p) 484 if not definite: 485 if self.prt: print 'indefinite or illconditioned step' 486 break 487 488 self.K[j]=1 489 490 if definite: # Find reduced Newton direction 491 p = numpy.zeros( self.n ) 492 ind = find( self.K > 0 ) 493 p[ind] = self.gradient[ind] 494 p = -mldiv( self.L.T, mldiv(self.L,p)/self.d ) 495 496 # Set tiny entries to zero 497 p = (self.x + p) - self.x 498 ind = find( p != 0 ) 499 if len(ind)==0: # Zero direction 500 if self.prt: print 'zero direction' 501 self.unfix=1 502 503 # Find range of step sizes 504 pp = p[ind] 505 oo = (self.xo[ind] - self.x[ind]) / pp 506 uu = (self.xu[ind] - self.x[ind]) / pp 507 508 if len(oo[pp<0])==0: max_oo = -inf 509 else: max_oo = max(oo[pp<0]) 510 if len(uu[pp>0])==0: max_uu = -inf 511 else: max_uu = max(uu[pp>0]) 512 alpha_u = numpy.maximum( max_oo, max_uu ) 513 514 if len(oo[pp>0])==0: min_oo = inf 515 else: min_oo = min(oo[pp>0]) 516 if len(uu[pp<0])==0: min_uu = inf 517 else: min_uu = min(uu[pp<0]) 518 alpha_o = numpy.minimum( min_oo, min_uu ) 519 520 if alpha_o<=0 or alpha_u>=0: 521 print 'programming error: no alpha' 522 523 # Find step size 524 gTp = float( numpy.dot(self.gradient.T, p) ) 525 agTp = float( numpy.dot(abs(self.gradient.T), abs(p)) ) 526 527 # Linear term consists of roundoff only 528 if abs(gTp) < self.neps * agTp: 529 gTp=0 530 531 pTGp = numpy.dot(p, self.G*toCol(p) ) 532 533 if self.convex: pTGp = numpy.max(0,pTGp) 534 if (not definite) and (pTGp>0): pTGp = 0 535 536 [alpha, lba, uba, info] = self._quadMin(alpha_u, alpha_o, gTp, pTGp) 537 if info: 538 self.x = numpy.zeros(n,1) 539 if lba: self.x=-p 540 else: self.x=p 541 542 # Allow variables to be freed in csearch? 543 self.unfix = not (lba or uba) 544 545 # Update of xx 546 for k in xrange(len(ind)): 547 # Avoid roundoff for active bounds 548 ik = ind[k] 549 if alpha == uu[k]: 550 self.xTry[ik] = self.xu[ik] 551 self.free[ik] = 0 552 553 elif alpha == oo[k]: 554 self.xTry[ik] = self.xo[ik] 555 self.free[ik] = 0 556 557 else: 558 self.xTry[ik] += alpha * p[ik] 559 560 if abs( self.xTry[ik] ) == inf: 561 print 'infinite xx in Minq' 562 563 self.nfree = numpy.sum( self.free ) 564 self.subdone = 1

565 566 567

568 - def subspaceSearch(self):

569 """ 570 Subspace search 571 """ 572 if not self.improvement or self.nitref > self.nitrefmax: 573 # Optimal point found - nothing done 574 pass 575 576 elif self.nitref > self.nitrefmax: 577 # Enough refinement steps - nothing done 578 pass 579 580 elif self.nfree == 0: 581 # No free variables - no subspace step taken 582 if self.prt>0: 583 print 'no free variables-no subspace step taken' 584 self.unfix=1 585 586 else: 587 self.subspace()

588 589

590 - def calcGain(self):

591 """ Calculate gain """ 592 return self.fval - self.gamma - \ 593 0.5* numpy.dot( self.x, self.c + toRow(self.gradient) )

594 595

596 - def calc_g(self):

597 """ Calculate gradient """ 598 return self.G * toCol(self.xTry) + toCol(self.c)

599 600

601 - def calc_newfval(self):

602 return self.gamma + \ 603 0.5 * numpy.dot( self.c + toRow(self.gradient), self.xTry)

604 605

606 - def pr01(self, name,x):

607 """ 608 Prints a (0,1) profile of x and returns the number of nonzeros 609 610 x: a numpy array 611 """ 612 # Check the size of the vector 613 n = x.size 614 text = name + ': ' 615 summe = 0 616 for k in xrange(n): 617 if x[k]: 618 text += '1' 619 summe += 1 620 else: 621 text += '0' 622 623 print text + ' ', summe, ' nonzeros' 624 return summe

625 626

627 - def search(self):

628 """ 629 Main loop: alternating coordinate and subspace searches 630 631 Outputs: 632 -------- 633 x minimizer (but unbounded direction if info=1) 634 fval optimal function value 635 nsub the number of subspace steps 636 info 0 (local minimizer found) 637 1 (unbounded below) 638 99 (maxiter exceeded) 639 """ 640 self.setSearchInit() 641 while 1: 642 643 self.gradient = self.calc_g() 644 self.fvalnew = self.calc_newfval() 645 646 if self.nitref==0: 647 self.fval = numpy.minimum(self.fval, self.fvalnew) 648 else: # More accurate g and hence f if nitref>0 649 self.fval = self.fvalnew 650 651 self.x = self.xTry 652 self._Info() 653 if self.info == 0 or self.info == 99: break 654 655 # ------------------------------------------------------------ 656 # Coordinate search 657 # ------------------------------------------------------------ 658 self.coordinateSearch() 659 if self.info: 660 return (self.x, self.fval, self.info, self.nsub) 661 662 # ============================================================ 663 self.nfree = numpy.sum( self.free ) 664 if self.unfix and self.nfree_old == self.nfree: 665 # In exact arithmetic, we are already optimal 666 # recompute gradient for iterative refinement 667 self.gradient = self.calc_g() 668 self.nitref += 1 669 if self.prt>0: print 'Optimum found.iterative refinement tried' 670 else: 671 self.nitref = 0 672 673 self.nfree_old = self.nfree 674 675 # Calculate gain 676 self.gain_cs = self.calcGain() 677 678 # Is improvement? 679 self.improvement = (self.gain_cs>0) or (not self.unfix) 680 681 if self.prt: 682 # print (0,1) profile of free and return the number of nonzeros 683 self.nfree = self.pr01('csrch ',self.free) 684 print self.gain_cs 685 686 # ------------------------------------------------------------- 687 # Subspace search 688 # ------------------------------------------------------------- 689 self.subspaceSearch() 690 if self.info: 691 return (self.x, self.fval, self.info, self.nsub) 692 693 # Show it 694 if self.prt: 695 # Print (0,1) profile of free and return the number of nonzeros 696 self.nfree = self.pr01('ssrch ', self.free) 697 if self.unfix and sum(self.nfree) < self.n: 698 print 'Bounds may be freed in next csearch' 699 700 return (self.x, self.fval, self.info, self.nsub)

Source Code for Module snobfit.minq