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