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