1 #include <petsctaolinesearch.h> 2 #include <../src/tao/matrix/lmvmmat.h> 3 #include <../src/tao/unconstrained/impls/lmvm/lmvm.h> 4 5 #define LMVM_BFGS 0 6 #define LMVM_SCALED_GRADIENT 1 7 #define LMVM_GRADIENT 2 8 9 #undef __FUNCT__ 10 #define __FUNCT__ "TaoSolve_LMVM" 11 static PetscErrorCode TaoSolve_LMVM(Tao tao) 12 { 13 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 14 PetscReal f, fold, gdx, gnorm; 15 PetscReal step = 1.0; 16 PetscReal delta; 17 PetscErrorCode ierr; 18 PetscInt stepType; 19 PetscInt iter = 0; 20 TaoConvergedReason reason = TAO_CONTINUE_ITERATING; 21 TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 22 23 PetscFunctionBegin; 24 25 if (tao->XL || tao->XU || tao->ops->computebounds) { 26 ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n");CHKERRQ(ierr); 27 } 28 29 /* Check convergence criteria */ 30 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr); 31 ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr); 32 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 33 34 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason);CHKERRQ(ierr); 35 if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 36 37 /* Set initial scaling for the function */ 38 if (f != 0.0) { 39 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 40 } else { 41 delta = 2.0 / (gnorm*gnorm); 42 } 43 ierr = MatLMVMSetDelta(lmP->M,delta);CHKERRQ(ierr); 44 45 /* Set counter for gradient/reset steps */ 46 lmP->bfgs = 0; 47 lmP->sgrad = 0; 48 lmP->grad = 0; 49 50 /* Have not converged; continue with Newton method */ 51 while (reason == TAO_CONTINUE_ITERATING) { 52 /* Compute direction */ 53 ierr = MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);CHKERRQ(ierr); 54 ierr = MatLMVMSolve(lmP->M, tao->gradient, lmP->D);CHKERRQ(ierr); 55 ++lmP->bfgs; 56 57 /* Check for success (descent direction) */ 58 ierr = VecDot(lmP->D, tao->gradient, &gdx);CHKERRQ(ierr); 59 if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) { 60 /* Step is not descent or direction produced not a number 61 We can assert bfgsUpdates > 1 in this case because 62 the first solve produces the scaled gradient direction, 63 which is guaranteed to be descent 64 65 Use steepest descent direction (scaled) 66 */ 67 68 ++lmP->grad; 69 70 if (f != 0.0) { 71 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 72 } else { 73 delta = 2.0 / (gnorm*gnorm); 74 } 75 ierr = MatLMVMSetDelta(lmP->M, delta);CHKERRQ(ierr); 76 ierr = MatLMVMReset(lmP->M);CHKERRQ(ierr); 77 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 78 ierr = MatLMVMSolve(lmP->M,tao->gradient, lmP->D);CHKERRQ(ierr); 79 80 /* On a reset, the direction cannot be not a number; it is a 81 scaled gradient step. No need to check for this condition. */ 82 83 lmP->bfgs = 1; 84 ++lmP->sgrad; 85 stepType = LMVM_SCALED_GRADIENT; 86 } else { 87 if (1 == lmP->bfgs) { 88 /* The first BFGS direction is always the scaled gradient */ 89 ++lmP->sgrad; 90 stepType = LMVM_SCALED_GRADIENT; 91 } else { 92 ++lmP->bfgs; 93 stepType = LMVM_BFGS; 94 } 95 } 96 ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr); 97 98 /* Perform the linesearch */ 99 fold = f; 100 ierr = VecCopy(tao->solution, lmP->Xold);CHKERRQ(ierr); 101 ierr = VecCopy(tao->gradient, lmP->Gold);CHKERRQ(ierr); 102 103 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step,&ls_status);CHKERRQ(ierr); 104 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 105 106 while (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_GRADIENT)) { 107 /* Linesearch failed */ 108 /* Reset factors and use scaled gradient step */ 109 f = fold; 110 ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr); 111 ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr); 112 113 switch(stepType) { 114 case LMVM_BFGS: 115 /* Failed to obtain acceptable iterate with BFGS step */ 116 /* Attempt to use the scaled gradient direction */ 117 118 if (f != 0.0) { 119 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 120 } else { 121 delta = 2.0 / (gnorm*gnorm); 122 } 123 ierr = MatLMVMSetDelta(lmP->M, delta);CHKERRQ(ierr); 124 ierr = MatLMVMReset(lmP->M);CHKERRQ(ierr); 125 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 126 ierr = MatLMVMSolve(lmP->M, tao->gradient, lmP->D);CHKERRQ(ierr); 127 128 /* On a reset, the direction cannot be not a number; it is a 129 scaled gradient step. No need to check for this condition. */ 130 131 lmP->bfgs = 1; 132 ++lmP->sgrad; 133 stepType = LMVM_SCALED_GRADIENT; 134 break; 135 136 case LMVM_SCALED_GRADIENT: 137 /* The scaled gradient step did not produce a new iterate; 138 attempt to use the gradient direction. 139 Need to make sure we are not using a different diagonal scaling */ 140 ierr = MatLMVMSetDelta(lmP->M, 1.0);CHKERRQ(ierr); 141 ierr = MatLMVMReset(lmP->M);CHKERRQ(ierr); 142 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 143 ierr = MatLMVMSolve(lmP->M, tao->gradient, lmP->D);CHKERRQ(ierr); 144 145 lmP->bfgs = 1; 146 ++lmP->grad; 147 stepType = LMVM_GRADIENT; 148 break; 149 } 150 ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr); 151 152 /* Perform the linesearch */ 153 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);CHKERRQ(ierr); 154 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 155 } 156 157 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 158 /* Failed to find an improving point */ 159 f = fold; 160 ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr); 161 ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr); 162 step = 0.0; 163 reason = TAO_DIVERGED_LS_FAILURE; 164 tao->reason = TAO_DIVERGED_LS_FAILURE; 165 } 166 /* Check for termination */ 167 ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr); 168 iter++; 169 ierr = TaoMonitor(tao,iter,f,gnorm,0.0,step,&reason);CHKERRQ(ierr); 170 } 171 PetscFunctionReturn(0); 172 } 173 174 #undef __FUNCT__ 175 #define __FUNCT__ "TaoSetUp_LMVM" 176 static PetscErrorCode TaoSetUp_LMVM(Tao tao) 177 { 178 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 179 PetscInt n,N; 180 PetscErrorCode ierr; 181 182 PetscFunctionBegin; 183 /* Existence of tao->solution checked in TaoSetUp() */ 184 if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr); } 185 if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr); } 186 if (!lmP->D) {ierr = VecDuplicate(tao->solution,&lmP->D);CHKERRQ(ierr); } 187 if (!lmP->Xold) {ierr = VecDuplicate(tao->solution,&lmP->Xold);CHKERRQ(ierr); } 188 if (!lmP->Gold) {ierr = VecDuplicate(tao->solution,&lmP->Gold);CHKERRQ(ierr); } 189 190 /* Create matrix for the limited memory approximation */ 191 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 192 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 193 ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&lmP->M);CHKERRQ(ierr); 194 ierr = MatLMVMAllocateVectors(lmP->M,tao->solution);CHKERRQ(ierr); 195 PetscFunctionReturn(0); 196 } 197 198 /* ---------------------------------------------------------- */ 199 #undef __FUNCT__ 200 #define __FUNCT__ "TaoDestroy_LMVM" 201 static PetscErrorCode TaoDestroy_LMVM(Tao tao) 202 { 203 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 204 PetscErrorCode ierr; 205 206 PetscFunctionBegin; 207 if (tao->setupcalled) { 208 ierr = VecDestroy(&lmP->Xold);CHKERRQ(ierr); 209 ierr = VecDestroy(&lmP->Gold);CHKERRQ(ierr); 210 ierr = VecDestroy(&lmP->D);CHKERRQ(ierr); 211 ierr = MatDestroy(&lmP->M);CHKERRQ(ierr); 212 } 213 ierr = PetscFree(tao->data);CHKERRQ(ierr); 214 PetscFunctionReturn(0); 215 } 216 217 /*------------------------------------------------------------*/ 218 #undef __FUNCT__ 219 #define __FUNCT__ "TaoSetFromOptions_LMVM" 220 static PetscErrorCode TaoSetFromOptions_LMVM(Tao tao) 221 { 222 PetscErrorCode ierr; 223 224 PetscFunctionBegin; 225 ierr = PetscOptionsHead("Limited-memory variable-metric method for unconstrained optimization");CHKERRQ(ierr); 226 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 227 ierr = PetscOptionsTail();CHKERRQ(ierr); 228 PetscFunctionReturn(0); 229 PetscFunctionReturn(0); 230 } 231 232 /*------------------------------------------------------------*/ 233 #undef __FUNCT__ 234 #define __FUNCT__ "TaoView_LMVM" 235 static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer) 236 { 237 TAO_LMVM *lm = (TAO_LMVM *)tao->data; 238 PetscBool isascii; 239 PetscErrorCode ierr; 240 241 PetscFunctionBegin; 242 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 243 if (isascii) { 244 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 245 ierr = PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);CHKERRQ(ierr); 246 ierr = PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);CHKERRQ(ierr); 247 ierr = PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);CHKERRQ(ierr); 248 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 249 } 250 PetscFunctionReturn(0); 251 } 252 253 /* ---------------------------------------------------------- */ 254 255 EXTERN_C_BEGIN 256 #undef __FUNCT__ 257 #define __FUNCT__ "TaoCreate_LMVM" 258 PetscErrorCode TaoCreate_LMVM(Tao tao) 259 { 260 TAO_LMVM *lmP; 261 const char *morethuente_type = TAOLINESEARCHMT; 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 tao->ops->setup = TaoSetUp_LMVM; 266 tao->ops->solve = TaoSolve_LMVM; 267 tao->ops->view = TaoView_LMVM; 268 tao->ops->setfromoptions = TaoSetFromOptions_LMVM; 269 tao->ops->destroy = TaoDestroy_LMVM; 270 271 ierr = PetscNewLog(tao,&lmP);CHKERRQ(ierr); 272 lmP->D = 0; 273 lmP->M = 0; 274 lmP->Xold = 0; 275 lmP->Gold = 0; 276 277 tao->data = (void*)lmP; 278 tao->max_it = 2000; 279 tao->max_funcs = 4000; 280 tao->fatol = 1e-4; 281 tao->frtol = 1e-4; 282 283 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);CHKERRQ(ierr); 284 ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr); 285 ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr); 286 PetscFunctionReturn(0); 287 } 288 EXTERN_C_END 289