1 #include <../src/tao/leastsquares/impls/brgn/brgn.h> 2 3 /* Old code 4 static PetscErrorCode GNHessianProd(Mat H, Vec in, Vec out) 5 { 6 TAO_BRGN *gn; 7 PetscErrorCode ierr; 8 9 PetscFunctionBegin; 10 ierr = MatShellGetContext(H, &gn);CHKERRQ(ierr); 11 ierr = MatMult(gn->subsolver->ls_jac, in, gn->r_work);CHKERRQ(ierr); 12 ierr = MatMultTranspose(gn->subsolver->ls_jac, gn->r_work, out);CHKERRQ(ierr); 13 ierr = VecAXPY(out, gn->lambda, in);CHKERRQ(ierr); 14 PetscFunctionReturn(0); 15 } 16 17 static PetscErrorCode GNObjectiveGradientEval(Tao tao, Vec X, PetscReal *fcn, Vec G, void *ptr) 18 { 19 TAO_BRGN *gn = (TAO_BRGN *)ptr; 20 PetscScalar xnorm2; 21 PetscErrorCode ierr; 22 23 PetscFunctionBegin; 24 ierr = TaoComputeResidual(tao, X, tao->ls_res);CHKERRQ(ierr); 25 ierr = VecDotBegin(tao->ls_res, tao->ls_res, fcn);CHKERRQ(ierr); 26 ierr = VecAXPBYPCZ(gn->x_work, 1.0, -1.0, 0.0, X, gn->x_old);CHKERRQ(ierr); 27 ierr = VecDotBegin(gn->x_work, gn->x_work, &xnorm2);CHKERRQ(ierr); 28 ierr = VecDotEnd(tao->ls_res, tao->ls_res, fcn);CHKERRQ(ierr); 29 ierr = VecDotEnd(gn->x_work, gn->x_work, &xnorm2);CHKERRQ(ierr); 30 *fcn = 0.5*(*fcn) + 0.5*gn->lambda*xnorm2; 31 32 ierr = TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre);CHKERRQ(ierr); 33 ierr = MatMultTranspose(tao->ls_jac, tao->ls_res, G);CHKERRQ(ierr); 34 ierr = VecAXPBYPCZ(G, gn->lambda, -gn->lambda, 1.0, X, gn->x_old);CHKERRQ(ierr); 35 PetscFunctionReturn(0); 36 } 37 */ 38 39 static PetscErrorCode GNHessianProd(Mat H, Vec in, Vec out) 40 { 41 TAO_BRGN *gn; 42 PetscErrorCode ierr; 43 44 PetscFunctionBegin; 45 ierr = MatShellGetContext(H, &gn);CHKERRQ(ierr); 46 ierr = MatMult(gn->subsolver->ls_jac, in, gn->r_work);CHKERRQ(ierr); 47 ierr = MatMultTranspose(gn->subsolver->ls_jac, gn->r_work, out);CHKERRQ(ierr); 48 /* out = out + lambda*in.*diag*/ 49 ierr = VecPointwiseMult(gn->x_work, in, gn->diag);CHKERRQ(ierr); /* gn->x_work = in.*diag, where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */ 50 ierr = VecAXPY(out, gn->lambda, gn->x_work);CHKERRQ(ierr); 51 52 PetscFunctionReturn(0); 53 } 54 55 static PetscErrorCode GNObjectiveGradientEval(Tao tao, Vec X, PetscReal *fcn, Vec G, void *ptr) 56 { 57 TAO_BRGN *gn = (TAO_BRGN *)ptr; 58 PetscInt N; /* dimension of X */ 59 PetscScalar xESum; 60 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 /* compute objective */ 64 /* compute first term ||ls_res||^2 */ 65 ierr = TaoComputeResidual(tao, X, tao->ls_res);CHKERRQ(ierr); 66 ierr = VecDotBegin(tao->ls_res, tao->ls_res, fcn);CHKERRQ(ierr); 67 ierr = VecDotEnd(tao->ls_res, tao->ls_res, fcn);CHKERRQ(ierr); 68 /* add the second term lambda*sum(sqrt(x.^2+epsilon^2) - epsilon)*/ 69 ierr = VecPointwiseMult(gn->x_work, X, X);CHKERRQ(ierr); 70 ierr = VecShift(gn->x_work, gn->epsilon*gn->epsilon);CHKERRQ(ierr); 71 ierr = VecSqrtAbs(gn->x_work);CHKERRQ(ierr); /* gn->x_work = sqrt(x.^2+epsilon^2) */ 72 ierr = VecSum(gn->x_work, &xESum);CHKERRQ(ierr);CHKERRQ(ierr); 73 ierr = VecGetSize(X, &N);CHKERRQ(ierr); 74 *fcn = 0.5*(*fcn) + gn->lambda*(xESum - N*gn->epsilon); 75 76 /* compute gradient G */ 77 ierr = TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre);CHKERRQ(ierr); 78 ierr = MatMultTranspose(tao->ls_jac, tao->ls_res, G);CHKERRQ(ierr); 79 /* compute G = G + lambda*(x./sqrt(x.^2+epsilon^2)) */ 80 ierr = VecPointwiseDivide(gn->x_work, X, gn->x_work);CHKERRQ(ierr); /* reuse x_work = x./sqrt(x.^2+epsilon^2) */ 81 ierr = VecAXPY(G, gn->lambda, gn->x_work);CHKERRQ(ierr); 82 83 PetscFunctionReturn(0); 84 } 85 86 87 static PetscErrorCode GNComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr) 88 { 89 TAO_BRGN *gn = (TAO_BRGN *)ptr; 90 PetscErrorCode ierr; 91 92 PetscFunctionBegin; 93 ierr = TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre);CHKERRQ(ierr); 94 95 /* calculate and store diagonal matrix as a vector: diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3*/ 96 ierr = VecPointwiseMult(gn->x_work, X, X);CHKERRQ(ierr); 97 ierr = VecShift(gn->x_work, gn->epsilon*gn->epsilon);CHKERRQ(ierr); 98 ierr = VecCopy(gn->x_work, gn->diag);CHKERRQ(ierr); /* gn->diag = x.^2+epsilon^2 */ 99 ierr = VecSqrtAbs(gn->x_work);CHKERRQ(ierr); /* gn->x_work = sqrt(x.^2+epsilon^2) */ 100 ierr = VecPointwiseMult(gn->diag, gn->x_work, gn->diag);CHKERRQ(ierr); /* gn->diag = sqrt(x.^2+epsilon^2).^3 */ 101 ierr = VecReciprocal(gn->diag);CHKERRQ(ierr); 102 ierr = VecScale(gn->diag, gn->epsilon*gn->epsilon);CHKERRQ(ierr); 103 104 PetscFunctionReturn(0); 105 } 106 107 static PetscErrorCode GNHookFunction(Tao tao, PetscInt iter) 108 { 109 TAO_BRGN *gn = (TAO_BRGN *)tao->user_update; 110 PetscErrorCode ierr; 111 112 PetscFunctionBegin; 113 /* Update basic tao information from the subsolver */ 114 gn->parent->nfuncs = tao->nfuncs; 115 gn->parent->ngrads = tao->ngrads; 116 gn->parent->nfuncgrads = tao->nfuncgrads; 117 gn->parent->nhess = tao->nhess; 118 gn->parent->niter = tao->niter; 119 gn->parent->ksp_its = tao->ksp_its; 120 gn->parent->ksp_tot_its = tao->ksp_tot_its; 121 ierr = TaoGetConvergedReason(tao, &gn->parent->reason);CHKERRQ(ierr); 122 /* Update the solution vectors */ 123 if (iter == 0) { 124 ierr = VecSet(gn->x_old, 0.0);CHKERRQ(ierr); 125 } else { 126 ierr = VecCopy(tao->solution, gn->x_old);CHKERRQ(ierr); 127 ierr = VecCopy(tao->solution, gn->parent->solution);CHKERRQ(ierr); 128 } 129 /* Update the gradient */ 130 ierr = VecCopy(tao->gradient, gn->parent->gradient);CHKERRQ(ierr); 131 /* Call general purpose update function */ 132 if (gn->parent->ops->update) { 133 ierr = (*gn->parent->ops->update)(gn->parent, gn->parent->niter);CHKERRQ(ierr); 134 } 135 PetscFunctionReturn(0); 136 } 137 138 static PetscErrorCode TaoSolve_BRGN(Tao tao) 139 { 140 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 141 PetscErrorCode ierr; 142 143 PetscFunctionBegin; 144 ierr = TaoSolve(gn->subsolver);CHKERRQ(ierr); 145 /* Update basic tao information from the subsolver */ 146 tao->nfuncs = gn->subsolver->nfuncs; 147 tao->ngrads = gn->subsolver->ngrads; 148 tao->nfuncgrads = gn->subsolver->nfuncgrads; 149 tao->nhess = gn->subsolver->nhess; 150 tao->niter = gn->subsolver->niter; 151 tao->ksp_its = gn->subsolver->ksp_its; 152 tao->ksp_tot_its = gn->subsolver->ksp_tot_its; 153 ierr = TaoGetConvergedReason(gn->subsolver, &tao->reason);CHKERRQ(ierr); 154 /* Update vectors */ 155 ierr = VecCopy(gn->subsolver->solution, tao->solution);CHKERRQ(ierr); 156 ierr = VecCopy(gn->subsolver->gradient, tao->gradient);CHKERRQ(ierr); 157 PetscFunctionReturn(0); 158 } 159 160 static PetscErrorCode TaoSetFromOptions_BRGN(PetscOptionItems *PetscOptionsObject,Tao tao) 161 { 162 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 163 PetscErrorCode ierr; 164 165 PetscFunctionBegin; 166 /* old Tikhonov regularization code 167 ierr = PetscOptionsHead(PetscOptionsObject,"Gauss-Newton method for least-squares problems using Tikhonov regularization");CHKERRQ(ierr); 168 ierr = PetscOptionsReal("-tao_brgn_lambda", "Tikhonov regularization factor", "", gn->lambda, &gn->lambda, NULL);CHKERRQ(ierr); 169 */ 170 ierr = PetscOptionsHead(PetscOptionsObject,"least-squares problems with L1 regularizer: ||f(x)||^2 + lambda*||x||_1. Currently L1-norm is approximated with smooth form");CHKERRQ(ierr); 171 ierr = PetscOptionsReal("-tao_brgn_lambda", "L1-norm regularizer weight", "", gn->lambda, &gn->lambda, NULL);CHKERRQ(ierr); 172 ierr = PetscOptionsReal("-tao_brgn_epsilon", "L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon)", "", gn->epsilon, &gn->epsilon, NULL);CHKERRQ(ierr); 173 ierr = PetscOptionsTail();CHKERRQ(ierr); 174 ierr = TaoSetFromOptions(gn->subsolver);CHKERRQ(ierr); 175 PetscFunctionReturn(0); 176 } 177 178 static PetscErrorCode TaoView_BRGN(Tao tao, PetscViewer viewer) 179 { 180 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 185 ierr = TaoView(gn->subsolver, viewer);CHKERRQ(ierr); 186 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 static PetscErrorCode TaoSetUp_BRGN(Tao tao) 191 { 192 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 193 PetscErrorCode ierr; 194 PetscBool is_bnls, is_bntr, is_bntl; 195 PetscInt i, nx, Nx; 196 197 PetscFunctionBegin; 198 if (!tao->ls_res) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetResidualRoutine() must be called before setup!"); 199 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNLS, &is_bnls);CHKERRQ(ierr); 200 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTR, &is_bntr);CHKERRQ(ierr); 201 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTL, &is_bntl);CHKERRQ(ierr); 202 if ((is_bnls || is_bntr || is_bntl) && !tao->ls_jac) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetResidualJacobianRoutine() must be called before setup!"); 203 if (!tao->gradient){ 204 ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr); 205 } 206 if (!gn->x_work){ 207 ierr = VecDuplicate(tao->solution, &gn->x_work);CHKERRQ(ierr); 208 } 209 if (!gn->r_work){ 210 ierr = VecDuplicate(tao->ls_res, &gn->r_work);CHKERRQ(ierr); 211 } 212 if (!gn->x_old) { 213 ierr = VecDuplicate(tao->solution, &gn->x_old);CHKERRQ(ierr); 214 ierr = VecSet(gn->x_old, 0.0);CHKERRQ(ierr); 215 } 216 if (!gn->diag){ 217 ierr = VecDuplicate(tao->solution, &gn->diag);CHKERRQ(ierr); 218 ierr = VecSet(gn->diag, 0.0);CHKERRQ(ierr); 219 } 220 221 if (!tao->setupcalled) { 222 /* Hessian setup */ 223 ierr = VecGetLocalSize(tao->solution, &nx);CHKERRQ(ierr); 224 ierr = VecGetSize(tao->solution, &Nx);CHKERRQ(ierr); 225 ierr = MatSetSizes(gn->H, nx, nx, Nx, Nx);CHKERRQ(ierr); 226 ierr = MatSetType(gn->H, MATSHELL);CHKERRQ(ierr); 227 ierr = MatSetUp(gn->H);CHKERRQ(ierr); 228 ierr = MatShellSetOperation(gn->H, MATOP_MULT, (void (*)(void))GNHessianProd);CHKERRQ(ierr); 229 ierr = MatShellSetContext(gn->H, (void*)gn);CHKERRQ(ierr); 230 /* Subsolver setup */ 231 ierr = TaoSetUpdate(gn->subsolver, GNHookFunction, (void*)gn);CHKERRQ(ierr); 232 ierr = TaoSetInitialVector(gn->subsolver, tao->solution);CHKERRQ(ierr); 233 if (tao->bounded) { 234 ierr = TaoSetVariableBounds(gn->subsolver, tao->XL, tao->XU);CHKERRQ(ierr); 235 } 236 ierr = TaoSetResidualRoutine(gn->subsolver, tao->ls_res, tao->ops->computeresidual, tao->user_lsresP);CHKERRQ(ierr); 237 ierr = TaoSetJacobianResidualRoutine(gn->subsolver, tao->ls_jac, tao->ls_jac, tao->ops->computeresidualjacobian, tao->user_lsjacP);CHKERRQ(ierr); 238 ierr = TaoSetObjectiveAndGradientRoutine(gn->subsolver, GNObjectiveGradientEval, (void*)gn);CHKERRQ(ierr); 239 ierr = TaoSetHessianRoutine(gn->subsolver, gn->H, gn->H, GNComputeHessian, (void*)gn);CHKERRQ(ierr); 240 /* Propagate some options down */ 241 ierr = TaoSetTolerances(gn->subsolver, tao->gatol, tao->grtol, tao->gttol);CHKERRQ(ierr); 242 ierr = TaoSetMaximumIterations(gn->subsolver, tao->max_it);CHKERRQ(ierr); 243 ierr = TaoSetMaximumFunctionEvaluations(gn->subsolver, tao->max_funcs);CHKERRQ(ierr); 244 for (i=0; i<tao->numbermonitors; ++i) { 245 ierr = TaoSetMonitor(gn->subsolver, tao->monitor[i], tao->monitorcontext[i], tao->monitordestroy[i]);CHKERRQ(ierr); 246 ierr = PetscObjectReference((PetscObject)(tao->monitorcontext[i]));CHKERRQ(ierr); 247 } 248 ierr = TaoSetUp(gn->subsolver);CHKERRQ(ierr); 249 } 250 PetscFunctionReturn(0); 251 } 252 253 static PetscErrorCode TaoDestroy_BRGN(Tao tao) 254 { 255 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 256 PetscErrorCode ierr; 257 258 PetscFunctionBegin; 259 if (tao->setupcalled) { 260 ierr = VecDestroy(&tao->gradient);CHKERRQ(ierr); 261 ierr = VecDestroy(&gn->x_work);CHKERRQ(ierr); 262 ierr = VecDestroy(&gn->r_work);CHKERRQ(ierr); 263 ierr = VecDestroy(&gn->x_old);CHKERRQ(ierr); 264 ierr = VecDestroy(&gn->diag);CHKERRQ(ierr); 265 } 266 ierr = MatDestroy(&gn->H);CHKERRQ(ierr); 267 ierr = TaoDestroy(&gn->subsolver);CHKERRQ(ierr); 268 gn->parent = NULL; 269 ierr = PetscFree(tao->data);CHKERRQ(ierr); 270 PetscFunctionReturn(0); 271 } 272 273 /*MC 274 TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares 275 problems with bound constraints. This algorithm is a thin wrapper around TAOBNTL 276 that constructs the Guass-Newton problem with the user-provided least-squares 277 residual and Jacobian. The problem is regularized with an L2-norm proximal point 278 term. 279 280 Options Database Keys: 281 + -tao_bqnk_max_cg_its - maximum number of bounded conjugate-gradient iterations taken in each Newton loop 282 . -tao_bqnk_init_type - trust radius initialization method ("constant", "direction", "interpolation") 283 . -tao_bqnk_update_type - trust radius update method ("step", "direction", "interpolation") 284 - -tao_bqnk_as_type - active-set estimation method ("none", "bertsekas") 285 286 Level: beginner 287 M*/ 288 PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao) 289 { 290 TAO_BRGN *gn; 291 PetscErrorCode ierr; 292 293 PetscFunctionBegin; 294 ierr = PetscNewLog(tao,&gn);CHKERRQ(ierr); 295 296 tao->ops->destroy = TaoDestroy_BRGN; 297 tao->ops->setup = TaoSetUp_BRGN; 298 tao->ops->setfromoptions = TaoSetFromOptions_BRGN; 299 tao->ops->view = TaoView_BRGN; 300 tao->ops->solve = TaoSolve_BRGN; 301 302 tao->data = (void*)gn; 303 gn->lambda = 1e-4; 304 gn->epsilon = 1e-6; 305 gn->parent = tao; 306 307 ierr = MatCreate(PetscObjectComm((PetscObject)tao), &gn->H);CHKERRQ(ierr); 308 ierr = MatSetOptionsPrefix(gn->H, "tao_brgn_hessian_");CHKERRQ(ierr); 309 310 ierr = TaoCreate(PetscObjectComm((PetscObject)tao), &gn->subsolver);CHKERRQ(ierr); 311 ierr = TaoSetType(gn->subsolver, TAOBNLS);CHKERRQ(ierr); 312 ierr = TaoSetOptionsPrefix(gn->subsolver, "tao_brgn_subsolver_");CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 /*@C 317 TaoBRGNGetSubsolver - Get the pointer to the subsolver inside BRGN 318 319 Collective on Tao 320 321 Level: developer 322 323 Input Parameters: 324 + tao - the Tao solver context 325 - subsolver - the Tao sub-solver context 326 @*/ 327 PetscErrorCode TaoBRGNGetSubsolver(Tao tao, Tao *subsolver) 328 { 329 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 330 331 PetscFunctionBegin; 332 *subsolver = gn->subsolver; 333 PetscFunctionReturn(0); 334 } 335 336 /*@C 337 TaoBRGNSetTikhonovLambda - Set the Tikhonov regularization factor for the Gauss-Newton least-squares algorithm 338 339 Collective on Tao 340 341 Level: developer 342 343 Input Parameters: 344 + tao - the Tao solver context 345 - lambda - Tikhonov regularization factor 346 @*/ 347 PetscErrorCode TaoBRGNSetTikhonovLambda(Tao tao, PetscReal lambda) 348 { 349 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 350 351 /* Initialize lambda here */ 352 353 PetscFunctionBegin; 354 gn->lambda = lambda; 355 PetscFunctionReturn(0); 356 } 357 358 /*@C 359 TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm 360 361 Collective on Tao 362 363 Level: developer 364 365 Input Parameters: 366 + tao - the Tao solver context 367 - epsilon - L1-norm smooth approximation parameter 368 @*/ 369 PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao, PetscReal epsilon) 370 { 371 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 372 373 /* Initialize epsilon here */ 374 375 PetscFunctionBegin; 376 gn->epsilon = epsilon; 377 PetscFunctionReturn(0); 378 } 379 /* XH: done TaoBRGNSetL1SmoothEpsilon by copy TaoBRGNSetTikhonovLambda in peststao.h, brgn.c and zbrgnf.c. 380 maybe change the name of Tikhonov in TaoBRGNSetTikhonovLambda() etc, as lambda is no longer the Tikhonov regularizer weight but the L1 regularizer weight */ 381