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) { 95 PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 96 } 97 98 /* Compute direction */ 99 PetscCall(MatLMVMUpdate(lmP->M,tao->solution,tao->gradient)); 100 PetscCall(MatSolve(lmP->M, lmP->GV, lmP->D)); 101 102 PetscCall(ProjDirect_OWLQN(lmP->D,lmP->GV)); 103 104 ++lmP->bfgs; 105 106 /* Check for success (descent direction) */ 107 PetscCall(VecDot(lmP->D, lmP->GV , &gdx)); 108 if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) { 109 110 /* Step is not descent or direction produced not a number 111 We can assert bfgsUpdates > 1 in this case because 112 the first solve produces the scaled gradient direction, 113 which is guaranteed to be descent 114 115 Use steepest descent direction (scaled) */ 116 ++lmP->grad; 117 118 delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm); 119 PetscCall(MatLMVMSetJ0Scale(lmP->M, delta)); 120 PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE)); 121 PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient)); 122 PetscCall(MatSolve(lmP->M,lmP->GV, lmP->D)); 123 124 PetscCall(ProjDirect_OWLQN(lmP->D,lmP->GV)); 125 126 lmP->bfgs = 1; 127 ++lmP->sgrad; 128 stepType = OWLQN_SCALED_GRADIENT; 129 } else { 130 if (1 == lmP->bfgs) { 131 /* The first BFGS direction is always the scaled gradient */ 132 ++lmP->sgrad; 133 stepType = OWLQN_SCALED_GRADIENT; 134 } else { 135 ++lmP->bfgs; 136 stepType = OWLQN_BFGS; 137 } 138 } 139 140 PetscCall(VecScale(lmP->D, -1.0)); 141 142 /* Perform the linesearch */ 143 fold = f; 144 PetscCall(VecCopy(tao->solution, lmP->Xold)); 145 PetscCall(VecCopy(tao->gradient, lmP->Gold)); 146 147 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step,&ls_status)); 148 PetscCall(TaoAddLineSearchCounts(tao)); 149 150 while (((int)ls_status < 0) && (stepType != OWLQN_GRADIENT)) { 151 152 /* Reset factors and use scaled gradient step */ 153 f = fold; 154 PetscCall(VecCopy(lmP->Xold, tao->solution)); 155 PetscCall(VecCopy(lmP->Gold, tao->gradient)); 156 PetscCall(VecCopy(tao->gradient, lmP->GV)); 157 158 PetscCall(ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda)); 159 160 switch(stepType) { 161 case OWLQN_BFGS: 162 /* Failed to obtain acceptable iterate with BFGS step 163 Attempt to use the scaled gradient direction */ 164 165 delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm); 166 PetscCall(MatLMVMSetJ0Scale(lmP->M, delta)); 167 PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE)); 168 PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient)); 169 PetscCall(MatSolve(lmP->M, lmP->GV, lmP->D)); 170 171 PetscCall(ProjDirect_OWLQN(lmP->D,lmP->GV)); 172 173 lmP->bfgs = 1; 174 ++lmP->sgrad; 175 stepType = OWLQN_SCALED_GRADIENT; 176 break; 177 178 case OWLQN_SCALED_GRADIENT: 179 /* The scaled gradient step did not produce a new iterate; 180 attempt to use the gradient direction. 181 Need to make sure we are not using a different diagonal scaling */ 182 PetscCall(MatLMVMSetJ0Scale(lmP->M, 1.0)); 183 PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE)); 184 PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient)); 185 PetscCall(MatSolve(lmP->M, lmP->GV, lmP->D)); 186 187 PetscCall(ProjDirect_OWLQN(lmP->D,lmP->GV)); 188 189 lmP->bfgs = 1; 190 ++lmP->grad; 191 stepType = OWLQN_GRADIENT; 192 break; 193 } 194 PetscCall(VecScale(lmP->D, -1.0)); 195 196 /* Perform the linesearch */ 197 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status)); 198 PetscCall(TaoAddLineSearchCounts(tao)); 199 } 200 201 if ((int)ls_status < 0) { 202 /* Failed to find an improving point*/ 203 f = fold; 204 PetscCall(VecCopy(lmP->Xold, tao->solution)); 205 PetscCall(VecCopy(lmP->Gold, tao->gradient)); 206 PetscCall(VecCopy(tao->gradient, lmP->GV)); 207 step = 0.0; 208 } else { 209 /* a little hack here, because that gv is used to store g */ 210 PetscCall(VecCopy(lmP->GV, tao->gradient)); 211 } 212 213 PetscCall(ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda)); 214 215 /* Check for termination */ 216 217 PetscCall(VecNorm(lmP->GV,NORM_2,&gnorm)); 218 219 iter++; 220 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 221 PetscCall(TaoMonitor(tao,iter,f,gnorm,0.0,step)); 222 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 223 224 if ((int)ls_status < 0) break; 225 } 226 PetscFunctionReturn(0); 227 } 228 229 static PetscErrorCode TaoSetUp_OWLQN(Tao tao) 230 { 231 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data; 232 PetscInt n,N; 233 234 PetscFunctionBegin; 235 /* Existence of tao->solution checked in TaoSetUp() */ 236 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution,&tao->gradient)); 237 if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution,&tao->stepdirection)); 238 if (!lmP->D) PetscCall(VecDuplicate(tao->solution,&lmP->D)); 239 if (!lmP->GV) PetscCall(VecDuplicate(tao->solution,&lmP->GV)); 240 if (!lmP->Xold) PetscCall(VecDuplicate(tao->solution,&lmP->Xold)); 241 if (!lmP->Gold) PetscCall(VecDuplicate(tao->solution,&lmP->Gold)); 242 243 /* Create matrix for the limited memory approximation */ 244 PetscCall(VecGetLocalSize(tao->solution,&n)); 245 PetscCall(VecGetSize(tao->solution,&N)); 246 PetscCall(MatCreateLMVMBFGS(((PetscObject)tao)->comm,n,N,&lmP->M)); 247 PetscCall(MatLMVMAllocate(lmP->M,tao->solution,tao->gradient)); 248 PetscFunctionReturn(0); 249 } 250 251 /* ---------------------------------------------------------- */ 252 static PetscErrorCode TaoDestroy_OWLQN(Tao tao) 253 { 254 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data; 255 256 PetscFunctionBegin; 257 if (tao->setupcalled) { 258 PetscCall(VecDestroy(&lmP->Xold)); 259 PetscCall(VecDestroy(&lmP->Gold)); 260 PetscCall(VecDestroy(&lmP->D)); 261 PetscCall(MatDestroy(&lmP->M)); 262 PetscCall(VecDestroy(&lmP->GV)); 263 } 264 PetscCall(PetscFree(tao->data)); 265 PetscFunctionReturn(0); 266 } 267 268 /*------------------------------------------------------------*/ 269 static PetscErrorCode TaoSetFromOptions_OWLQN(PetscOptionItems *PetscOptionsObject,Tao tao) 270 { 271 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data; 272 273 PetscFunctionBegin; 274 PetscOptionsHeadBegin(PetscOptionsObject,"Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization"); 275 PetscCall(PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight","", 100,&lmP->lambda,NULL)); 276 PetscOptionsHeadEnd(); 277 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch)); 278 PetscFunctionReturn(0); 279 } 280 281 /*------------------------------------------------------------*/ 282 static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer) 283 { 284 TAO_OWLQN *lm = (TAO_OWLQN *)tao->data; 285 PetscBool isascii; 286 287 PetscFunctionBegin; 288 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 289 if (isascii) { 290 PetscCall(PetscViewerASCIIPushTab(viewer)); 291 PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS steps: %" PetscInt_FMT "\n", lm->bfgs)); 292 PetscCall(PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %" PetscInt_FMT "\n", lm->sgrad)); 293 PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lm->grad)); 294 PetscCall(PetscViewerASCIIPopTab(viewer)); 295 } 296 PetscFunctionReturn(0); 297 } 298 299 /* ---------------------------------------------------------- */ 300 /*MC 301 TAOOWLQN - orthant-wise limited memory quasi-newton algorithm 302 303 . - tao_owlqn_lambda - regulariser weight 304 305 Level: beginner 306 M*/ 307 308 PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao) 309 { 310 TAO_OWLQN *lmP; 311 const char *owarmijo_type = TAOLINESEARCHOWARMIJO; 312 313 PetscFunctionBegin; 314 tao->ops->setup = TaoSetUp_OWLQN; 315 tao->ops->solve = TaoSolve_OWLQN; 316 tao->ops->view = TaoView_OWLQN; 317 tao->ops->setfromoptions = TaoSetFromOptions_OWLQN; 318 tao->ops->destroy = TaoDestroy_OWLQN; 319 320 PetscCall(PetscNewLog(tao,&lmP)); 321 lmP->D = NULL; 322 lmP->M = NULL; 323 lmP->GV = NULL; 324 lmP->Xold = NULL; 325 lmP->Gold = NULL; 326 lmP->lambda = 1.0; 327 328 tao->data = (void*)lmP; 329 /* Override default settings (unless already changed) */ 330 if (!tao->max_it_changed) tao->max_it = 2000; 331 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 332 333 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch)); 334 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 335 PetscCall(TaoLineSearchSetType(tao->linesearch,owarmijo_type)); 336 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch,tao)); 337 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix)); 338 PetscFunctionReturn(0); 339 } 340