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