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