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