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