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