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