1 #include "taolinesearch.h" 2 #include "../src/tao/matrix/lmvmmat.h" 3 #include "owlqn.h" 4 5 #define OWLQN_BFGS 0 6 #define OWLQN_SCALED_GRADIENT 1 7 #define OWLQN_GRADIENT 2 8 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "ProjDirect_OWLQN" 12 static PetscErrorCode ProjDirect_OWLQN(Vec d, Vec g) 13 { 14 PetscErrorCode ierr; 15 PetscReal *gptr,*dptr; 16 PetscInt low,high,low1,high1,i; 17 18 PetscFunctionBegin; 19 20 ierr=VecGetOwnershipRange(d,&low,&high);CHKERRQ(ierr); 21 ierr=VecGetOwnershipRange(g,&low1,&high1);CHKERRQ(ierr); 22 23 ierr = VecGetArray(g,&gptr);CHKERRQ(ierr); 24 ierr = VecGetArray(d,&dptr);CHKERRQ(ierr); 25 for (i = 0; i < high-low; i++) { 26 if (dptr[i] * gptr[i] <= 0.0 ) 27 { 28 dptr[i] = 0.0; 29 } 30 } 31 ierr = VecRestoreArray(d,&dptr);CHKERRQ(ierr); 32 ierr = VecRestoreArray(g,&gptr);CHKERRQ(ierr); 33 34 35 PetscFunctionReturn(0); 36 } 37 38 #undef __FUNCT__ 39 #define __FUNCT__ "ComputePseudoGrad_OWLQN" 40 static PetscErrorCode ComputePseudoGrad_OWLQN(Vec x, Vec gv, PetscReal lambda) 41 { 42 43 PetscErrorCode ierr; 44 PetscReal *xptr,*gptr; 45 PetscInt low,high,low1,high1,i; 46 47 48 PetscFunctionBegin; 49 50 ierr=VecGetOwnershipRange(x,&low,&high);CHKERRQ(ierr); 51 ierr=VecGetOwnershipRange(gv,&low1,&high1);CHKERRQ(ierr); 52 53 ierr = VecGetArray(x,&xptr);CHKERRQ(ierr); 54 ierr = VecGetArray(gv,&gptr);CHKERRQ(ierr); 55 for (i = 0; i < high-low; i++) { 56 if (xptr[i] < 0.0) 57 gptr[i] = gptr[i] - lambda; 58 else if (xptr[i] > 0.0) 59 gptr[i] = gptr[i] + lambda; 60 else if (gptr[i] + lambda < 0.0) 61 gptr[i] = gptr[i] + lambda; 62 else if (gptr[i] - lambda > 0.0) 63 gptr[i] = gptr[i] - lambda; 64 else 65 gptr[i] = 0.0; 66 } 67 ierr = VecRestoreArray(gv,&gptr);CHKERRQ(ierr); 68 ierr = VecRestoreArray(x,&xptr);CHKERRQ(ierr); 69 70 PetscFunctionReturn(0); 71 } 72 73 74 75 76 #undef __FUNCT__ 77 #define __FUNCT__ "TaoSolve_OWLQN" 78 static PetscErrorCode TaoSolve_OWLQN(TaoSolver tao) 79 { 80 81 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data; 82 83 PetscReal f, fold, gdx, gnorm; 84 PetscReal step = 1.0; 85 86 PetscReal delta; 87 88 PetscErrorCode ierr; 89 PetscInt stepType; 90 PetscInt iter = 0; 91 TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING; 92 TaoLineSearchTerminationReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 93 94 PetscFunctionBegin; 95 96 if (tao->XL || tao->XU || tao->ops->computebounds) { 97 ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by owlqn algorithm\n"); CHKERRQ(ierr); 98 } 99 100 /* Check convergence criteria */ 101 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient); CHKERRQ(ierr); 102 103 ierr = VecCopy(tao->gradient, lmP->GV); CHKERRQ(ierr); 104 105 ierr = ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);CHKERRQ(ierr); 106 107 ierr = VecNorm(lmP->GV,NORM_2,&gnorm); CHKERRQ(ierr); 108 109 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) { 110 SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 111 } 112 113 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason); CHKERRQ(ierr); 114 if (reason != TAO_CONTINUE_ITERATING) { 115 PetscFunctionReturn(0); 116 } 117 118 /* Set initial scaling for the function */ 119 if (f != 0.0) { 120 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 121 } 122 else { 123 delta = 2.0 / (gnorm*gnorm); 124 } 125 ierr = MatLMVMSetDelta(lmP->M,delta); CHKERRQ(ierr); 126 127 /* Set counter for gradient/reset steps */ 128 lmP->bfgs = 0; 129 lmP->sgrad = 0; 130 lmP->grad = 0; 131 132 /* Have not converged; continue with Newton method */ 133 while (reason == TAO_CONTINUE_ITERATING) { 134 135 /* Compute direction */ 136 ierr = MatLMVMUpdate(lmP->M,tao->solution,tao->gradient); CHKERRQ(ierr); 137 ierr = MatLMVMSolve(lmP->M, lmP->GV, lmP->D); CHKERRQ(ierr); 138 139 ierr = ProjDirect_OWLQN(lmP->D,lmP->GV);CHKERRQ(ierr); 140 141 ++lmP->bfgs; 142 143 /* Check for success (descent direction) */ 144 ierr = VecDot(lmP->D, lmP->GV , &gdx); CHKERRQ(ierr); 145 if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) { 146 147 /* Step is not descent or direction produced not a number 148 We can assert bfgsUpdates > 1 in this case because 149 the first solve produces the scaled gradient direction, 150 which is guaranteed to be descent 151 152 Use steepest descent direction (scaled) */ 153 ++lmP->grad; 154 155 if (f != 0.0) { 156 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 157 } 158 else { 159 delta = 2.0 / (gnorm*gnorm); 160 } 161 ierr = MatLMVMSetDelta(lmP->M, delta); CHKERRQ(ierr); 162 ierr = MatLMVMReset(lmP->M); CHKERRQ(ierr); 163 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient); CHKERRQ(ierr); 164 ierr = MatLMVMSolve(lmP->M,lmP->GV, lmP->D); CHKERRQ(ierr); 165 166 ierr = ProjDirect_OWLQN(lmP->D,lmP->GV);CHKERRQ(ierr); 167 /* On a reset, the direction cannot be not a number; it is a 168 scaled gradient step. No need to check for this condition. 169 info = D->Norm2(&dnorm); CHKERRQ(info); 170 if (PetscIsInfOrNanReal(dnorm)) { 171 SETERRQ(PETSC_COMM_SELF,1, "Direction generated Not-a-Number"); 172 } */ 173 174 lmP->bfgs = 1; 175 ++lmP->sgrad; 176 stepType = OWLQN_SCALED_GRADIENT; 177 } 178 else { 179 if (1 == lmP->bfgs) { 180 /* The first BFGS direction is always the scaled gradient */ 181 ++lmP->sgrad; 182 stepType = OWLQN_SCALED_GRADIENT; 183 } 184 else { 185 ++lmP->bfgs; 186 stepType = OWLQN_BFGS; 187 } 188 } 189 190 ierr = VecScale(lmP->D, -1.0); CHKERRQ(ierr); 191 192 /* Perform the linesearch */ 193 fold = f; 194 ierr = VecCopy(tao->solution, lmP->Xold); CHKERRQ(ierr); 195 ierr = VecCopy(tao->gradient, lmP->Gold); CHKERRQ(ierr); 196 197 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step,&ls_status); CHKERRQ(ierr); 198 ierr = TaoAddLineSearchCounts(tao); CHKERRQ(ierr); 199 200 while (((int)ls_status < 0) && (stepType != OWLQN_GRADIENT)) { 201 202 /* Reset factors and use scaled gradient step */ 203 f = fold; 204 ierr = VecCopy(lmP->Xold, tao->solution); CHKERRQ(ierr); 205 ierr = VecCopy(lmP->Gold, tao->gradient); CHKERRQ(ierr); 206 ierr = VecCopy(tao->gradient, lmP->GV); CHKERRQ(ierr); 207 208 ierr = ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);CHKERRQ(ierr); 209 210 switch(stepType) { 211 case OWLQN_BFGS: 212 /* Failed to obtain acceptable iterate with BFGS step 213 Attempt to use the scaled gradient direction */ 214 215 if (f != 0.0) { 216 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 217 } 218 else { 219 delta = 2.0 / (gnorm*gnorm); 220 } 221 ierr = MatLMVMSetDelta(lmP->M, delta); CHKERRQ(ierr); 222 ierr = MatLMVMReset(lmP->M); CHKERRQ(ierr); 223 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient); CHKERRQ(ierr); 224 ierr = MatLMVMSolve(lmP->M, lmP->GV, lmP->D); CHKERRQ(ierr); 225 226 ierr = ProjDirect_OWLQN(lmP->D,lmP->GV);CHKERRQ(ierr); 227 /* On a reset, the direction cannot be not a number; it is a 228 scaled gradient step. No need to check for this condition. 229 info = D->Norm2(&dnorm); CHKERRQ(info); 230 if (PetscIsInfOrNanReal(dnorm)) { 231 SETERRQ(PETSC_COMM_SELF,1, "Direction generated Not-a-Number"); 232 }*/ 233 234 lmP->bfgs = 1; 235 ++lmP->sgrad; 236 stepType = OWLQN_SCALED_GRADIENT; 237 break; 238 239 case OWLQN_SCALED_GRADIENT: 240 /* The scaled gradient step did not produce a new iterate; 241 attempt to use the gradient direction. 242 Need to make sure we are not using a different diagonal scaling */ 243 ierr = MatLMVMSetDelta(lmP->M, 1.0); CHKERRQ(ierr); 244 ierr = MatLMVMReset(lmP->M); CHKERRQ(ierr); 245 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient); CHKERRQ(ierr); 246 ierr = MatLMVMSolve(lmP->M, lmP->GV, lmP->D); CHKERRQ(ierr); 247 248 ierr = ProjDirect_OWLQN(lmP->D,lmP->GV);CHKERRQ(ierr); 249 250 lmP->bfgs = 1; 251 ++lmP->grad; 252 stepType = OWLQN_GRADIENT; 253 break; 254 } 255 ierr = VecScale(lmP->D, -1.0); CHKERRQ(ierr); 256 257 258 /* Perform the linesearch */ 259 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status); CHKERRQ(ierr); 260 ierr = TaoAddLineSearchCounts(tao); CHKERRQ(ierr); 261 } 262 263 if ((int)ls_status < 0) { 264 /* Failed to find an improving point*/ 265 f = fold; 266 ierr = VecCopy(lmP->Xold, tao->solution); CHKERRQ(ierr); 267 ierr = VecCopy(lmP->Gold, tao->gradient); CHKERRQ(ierr); 268 ierr = VecCopy(tao->gradient, lmP->GV); CHKERRQ(ierr); 269 step = 0.0; 270 } 271 else { 272 /* a little hack here, because that gv is used to store g */ 273 ierr = VecCopy(lmP->GV, tao->gradient); CHKERRQ(ierr); 274 } 275 276 ierr = ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);CHKERRQ(ierr); 277 278 /* Check for termination */ 279 280 ierr = VecNorm(lmP->GV,NORM_2,&gnorm); CHKERRQ(ierr); 281 282 iter++; 283 ierr = TaoMonitor(tao,iter,f,gnorm,0.0,step,&reason); CHKERRQ(ierr); 284 285 if ((int)ls_status < 0) 286 break; 287 } 288 PetscFunctionReturn(0); 289 } 290 291 292 #undef __FUNCT__ 293 #define __FUNCT__ "TaoSetUp_OWLQN" 294 static PetscErrorCode TaoSetUp_OWLQN(TaoSolver tao) 295 { 296 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data; 297 PetscInt n,N; 298 PetscErrorCode ierr; 299 300 PetscFunctionBegin; 301 /* Existence of tao->solution checked in TaoSolverSetUp() */ 302 if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient); CHKERRQ(ierr); } 303 if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection); CHKERRQ(ierr); } 304 if (!lmP->D) {ierr = VecDuplicate(tao->solution,&lmP->D); CHKERRQ(ierr); } 305 if (!lmP->GV) {ierr = VecDuplicate(tao->solution,&lmP->GV); CHKERRQ(ierr); } 306 if (!lmP->Xold) {ierr = VecDuplicate(tao->solution,&lmP->Xold); CHKERRQ(ierr); } 307 if (!lmP->Gold) {ierr = VecDuplicate(tao->solution,&lmP->Gold); CHKERRQ(ierr); } 308 309 /* Create matrix for the limited memory approximation */ 310 ierr = VecGetLocalSize(tao->solution,&n); CHKERRQ(ierr); 311 ierr = VecGetSize(tao->solution,&N); CHKERRQ(ierr); 312 ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&lmP->M); CHKERRQ(ierr); 313 ierr = MatLMVMAllocateVectors(lmP->M,tao->solution); CHKERRQ(ierr); 314 315 316 PetscFunctionReturn(0); 317 } 318 319 /* ---------------------------------------------------------- */ 320 #undef __FUNCT__ 321 #define __FUNCT__ "TaoDestroy_OWLQN" 322 static PetscErrorCode TaoDestroy_OWLQN(TaoSolver tao) 323 { 324 325 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data; 326 PetscErrorCode ierr; 327 328 PetscFunctionBegin; 329 if (tao->setupcalled) { 330 ierr = VecDestroy(&lmP->Xold); CHKERRQ(ierr); 331 ierr = VecDestroy(&lmP->Gold); CHKERRQ(ierr); 332 ierr = VecDestroy(&lmP->D); CHKERRQ(ierr); 333 ierr = MatDestroy(&lmP->M); CHKERRQ(ierr); 334 ierr = VecDestroy(&lmP->GV);CHKERRQ(ierr); 335 } 336 ierr = PetscFree(tao->data); CHKERRQ(ierr); 337 tao->data = PETSC_NULL; 338 339 PetscFunctionReturn(0); 340 } 341 342 /*------------------------------------------------------------*/ 343 #undef __FUNCT__ 344 #define __FUNCT__ "TaoSetFromOptions_OWLQN" 345 static PetscErrorCode TaoSetFromOptions_OWLQN(TaoSolver tao) 346 { 347 348 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data; 349 PetscBool flg; 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 ierr = PetscOptionsHead("Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization"); CHKERRQ(ierr); 354 ierr = PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight","", 100,&lmP->lambda,&flg); CHKERRQ(ierr); 355 ierr = PetscOptionsTail(); CHKERRQ(ierr); 356 ierr = TaoLineSearchSetFromOptions(tao->linesearch); CHKERRQ(ierr); 357 PetscFunctionReturn(0); 358 } 359 360 /*------------------------------------------------------------*/ 361 #undef __FUNCT__ 362 #define __FUNCT__ "TaoView_OWLQN" 363 static PetscErrorCode TaoView_OWLQN(TaoSolver tao, PetscViewer viewer) 364 { 365 366 TAO_OWLQN *lm = (TAO_OWLQN *)tao->data; 367 PetscBool isascii; 368 PetscErrorCode ierr; 369 370 371 PetscFunctionBegin; 372 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii); CHKERRQ(ierr); 373 if (isascii) { 374 375 ierr = PetscViewerASCIIPushTab(viewer); CHKERRQ(ierr); 376 ierr = PetscViewerASCIIPrintf(viewer, "BFGS steps: %d\n", lm->bfgs); CHKERRQ(ierr); 377 ierr = PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %d\n", lm->sgrad); CHKERRQ(ierr); 378 ierr = PetscViewerASCIIPrintf(viewer, "Gradient steps: %d\n", lm->grad); CHKERRQ(ierr); 379 ierr = PetscViewerASCIIPopTab(viewer); CHKERRQ(ierr); 380 } else { 381 SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO OWLQN",((PetscObject)viewer)->type_name); 382 } 383 PetscFunctionReturn(0); 384 } 385 386 /* ---------------------------------------------------------- */ 387 388 EXTERN_C_BEGIN 389 #undef __FUNCT__ 390 #define __FUNCT__ "TaoCreate_OWLQN" 391 PetscErrorCode TaoCreate_OWLQN(TaoSolver tao) 392 { 393 394 TAO_OWLQN *lmP; 395 const char *owarmijo_type = TAOLINESEARCH_OWARMIJO; 396 PetscErrorCode ierr; 397 398 PetscFunctionBegin; 399 tao->ops->setup = TaoSetUp_OWLQN; 400 tao->ops->solve = TaoSolve_OWLQN; 401 tao->ops->view = TaoView_OWLQN; 402 tao->ops->setfromoptions = TaoSetFromOptions_OWLQN; 403 tao->ops->destroy = TaoDestroy_OWLQN; 404 405 ierr = PetscNewLog(tao,TAO_OWLQN, &lmP); CHKERRQ(ierr); 406 lmP->D = 0; 407 lmP->M = 0; 408 lmP->GV = 0; 409 lmP->Xold = 0; 410 lmP->Gold = 0; 411 lmP->lambda = 1.0; 412 413 tao->data = (void*)lmP; 414 tao->max_it = 2000; 415 tao->max_funcs = 4000; 416 tao->fatol = 1e-4; 417 tao->frtol = 1e-4; 418 419 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch); CHKERRQ(ierr); 420 ierr = TaoLineSearchSetType(tao->linesearch,owarmijo_type); CHKERRQ(ierr); 421 ierr = TaoLineSearchUseTaoSolverRoutines(tao->linesearch,tao); CHKERRQ(ierr); 422 423 424 PetscFunctionReturn(0); 425 } 426 EXTERN_C_END 427 428