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 infinity 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 static PetscErrorCode TaoDestroy_OWLQN(Tao tao) 243 { 244 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data; 245 246 PetscFunctionBegin; 247 if (tao->setupcalled) { 248 PetscCall(VecDestroy(&lmP->Xold)); 249 PetscCall(VecDestroy(&lmP->Gold)); 250 PetscCall(VecDestroy(&lmP->D)); 251 PetscCall(MatDestroy(&lmP->M)); 252 PetscCall(VecDestroy(&lmP->GV)); 253 } 254 PetscCall(PetscFree(tao->data)); 255 PetscFunctionReturn(PETSC_SUCCESS); 256 } 257 258 static PetscErrorCode TaoSetFromOptions_OWLQN(Tao tao, PetscOptionItems PetscOptionsObject) 259 { 260 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data; 261 262 PetscFunctionBegin; 263 PetscOptionsHeadBegin(PetscOptionsObject, "Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization"); 264 PetscCall(PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight", "", 100, &lmP->lambda, NULL)); 265 PetscOptionsHeadEnd(); 266 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch)); 267 PetscFunctionReturn(PETSC_SUCCESS); 268 } 269 270 static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer) 271 { 272 TAO_OWLQN *lm = (TAO_OWLQN *)tao->data; 273 PetscBool isascii; 274 275 PetscFunctionBegin; 276 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 277 if (isascii) { 278 PetscCall(PetscViewerASCIIPushTab(viewer)); 279 PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS steps: %" PetscInt_FMT "\n", lm->bfgs)); 280 PetscCall(PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %" PetscInt_FMT "\n", lm->sgrad)); 281 PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lm->grad)); 282 PetscCall(PetscViewerASCIIPopTab(viewer)); 283 } 284 PetscFunctionReturn(PETSC_SUCCESS); 285 } 286 287 /*MC 288 TAOOWLQN - orthant-wise limited memory quasi-newton algorithm 289 290 . - tao_owlqn_lambda - regulariser weight 291 292 Level: beginner 293 M*/ 294 295 PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao) 296 { 297 TAO_OWLQN *lmP; 298 const char *owarmijo_type = TAOLINESEARCHOWARMIJO; 299 300 PetscFunctionBegin; 301 tao->ops->setup = TaoSetUp_OWLQN; 302 tao->ops->solve = TaoSolve_OWLQN; 303 tao->ops->view = TaoView_OWLQN; 304 tao->ops->setfromoptions = TaoSetFromOptions_OWLQN; 305 tao->ops->destroy = TaoDestroy_OWLQN; 306 307 PetscCall(PetscNew(&lmP)); 308 lmP->D = NULL; 309 lmP->M = NULL; 310 lmP->GV = NULL; 311 lmP->Xold = NULL; 312 lmP->Gold = NULL; 313 lmP->lambda = 1.0; 314 315 tao->data = (void *)lmP; 316 /* Override default settings (unless already changed) */ 317 PetscCall(TaoParametersInitialize(tao)); 318 PetscObjectParameterSetDefault(tao, max_it, 2000); 319 PetscObjectParameterSetDefault(tao, max_funcs, 4000); 320 321 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 322 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 323 PetscCall(TaoLineSearchSetType(tao->linesearch, owarmijo_type)); 324 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao)); 325 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix)); 326 PetscFunctionReturn(PETSC_SUCCESS); 327 } 328