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