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