1 #include <../src/tao/leastsquares/impls/brgn/brgn.h> 2 3 #define BRGN_user 0 4 #define BRGN_l2prox 1 5 #define BRGN_l1dict 2 6 #define BRGNRegTypes 3 7 8 static const char *BRGN_Table[64] = {"user","l2prox","l1dict"}; 9 10 static PetscErrorCode GNHessianProd(Mat H,Vec in,Vec out) 11 { 12 TAO_BRGN *gn; 13 PetscErrorCode ierr; 14 15 PetscFunctionBegin; 16 ierr = MatShellGetContext(H,&gn);CHKERRQ(ierr); 17 ierr = MatMult(gn->subsolver->ls_jac,in,gn->r_work);CHKERRQ(ierr); 18 ierr = MatMultTranspose(gn->subsolver->ls_jac,gn->r_work,out);CHKERRQ(ierr); 19 switch (gn->reg_type) { 20 case BRGN_user: 21 ierr = MatMult(gn->Hreg,in,gn->x_work);CHKERRQ(ierr); 22 ierr = VecAXPY(out,gn->lambda,gn->x_work);CHKERRQ(ierr); 23 break; 24 case BRGN_l2prox: 25 ierr = VecAXPY(out,gn->lambda,in);CHKERRQ(ierr); 26 break; 27 case BRGN_l1dict: 28 /* out = out + lambda*D'*(diag.*(D*in)) */ 29 if (gn->D) { 30 ierr = MatMult(gn->D,in,gn->y);CHKERRQ(ierr);/* y = D*in */ 31 } else { 32 ierr = VecCopy(in,gn->y);CHKERRQ(ierr); 33 } 34 ierr = VecPointwiseMult(gn->y_work,gn->diag,gn->y);CHKERRQ(ierr); /* y_work = diag.*(D*in), where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */ 35 if (gn->D) { 36 ierr = MatMultTranspose(gn->D,gn->y_work,gn->x_work);CHKERRQ(ierr); /* x_work = D'*(diag.*(D*in)) */ 37 } else { 38 ierr = VecCopy(gn->y_work,gn->x_work);CHKERRQ(ierr); 39 } 40 ierr = VecAXPY(out,gn->lambda,gn->x_work);CHKERRQ(ierr); 41 break; 42 } 43 44 PetscFunctionReturn(0); 45 } 46 47 static PetscErrorCode GNObjectiveGradientEval(Tao tao,Vec X,PetscReal *fcn,Vec G,void *ptr) 48 { 49 TAO_BRGN *gn = (TAO_BRGN *)ptr; 50 PetscInt K; /* dimension of D*X */ 51 PetscScalar yESum; 52 PetscErrorCode ierr; 53 PetscReal f_reg; 54 55 PetscFunctionBegin; 56 /* compute objective *fcn*/ 57 /* compute first term 0.5*||ls_res||_2^2 */ 58 ierr = TaoComputeResidual(tao,X,tao->ls_res);CHKERRQ(ierr); 59 ierr = VecDot(tao->ls_res,tao->ls_res,fcn);CHKERRQ(ierr); 60 *fcn *= 0.5; 61 /* compute gradient G */ 62 ierr = TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);CHKERRQ(ierr); 63 ierr = MatMultTranspose(tao->ls_jac,tao->ls_res,G);CHKERRQ(ierr); 64 /* add the regularization contribution */ 65 switch (gn->reg_type) { 66 case BRGN_user: 67 ierr = (*gn->regularizerobjandgrad)(tao,X,&f_reg,gn->x_work,gn->reg_obj_ctx);CHKERRQ(ierr); 68 *fcn += gn->lambda*f_reg; 69 ierr = VecAXPY(G,gn->lambda,gn->x_work);CHKERRQ(ierr); 70 break; 71 case BRGN_l2prox: 72 /* compute f = f + lambda*0.5*(xk - xkm1)^T(xk - xkm1) */ 73 ierr = VecAXPBYPCZ(gn->x_work,1.0,-1.0,0.0,X,gn->x_old);CHKERRQ(ierr); /*TODO: no need to use VecAXPBYPCZ for x - xkm1 */ 74 ierr = VecDot(gn->x_work,gn->x_work,&f_reg);CHKERRQ(ierr); 75 *fcn += gn->lambda*0.5*f_reg; 76 /* compute G = G + lambda*(xk - xkm1) */ 77 ierr = VecAXPBYPCZ(G,gn->lambda,-gn->lambda,1.0,X,gn->x_old);CHKERRQ(ierr); 78 break; 79 case BRGN_l1dict: 80 /* compute f = f + lambda*sum(sqrt(y.^2+epsilon^2) - epsilon), where y = D*x*/ 81 if (gn->D) { 82 ierr = MatMult(gn->D,X,gn->y);CHKERRQ(ierr);/* y = D*x */ 83 } else { 84 ierr = VecCopy(X,gn->y);CHKERRQ(ierr); 85 } 86 ierr = VecPointwiseMult(gn->y_work,gn->y,gn->y);CHKERRQ(ierr); 87 ierr = VecShift(gn->y_work,gn->epsilon*gn->epsilon);CHKERRQ(ierr); 88 ierr = VecSqrtAbs(gn->y_work);CHKERRQ(ierr); /* gn->y_work = sqrt(y.^2+epsilon^2) */ 89 ierr = VecSum(gn->y_work,&yESum);CHKERRQ(ierr);CHKERRQ(ierr); 90 ierr = VecGetSize(gn->y,&K);CHKERRQ(ierr); 91 *fcn += gn->lambda*(yESum - K*gn->epsilon); 92 /* compute G = G + lambda*D'*(y./sqrt(y.^2+epsilon^2)),where y = D*x */ 93 ierr = VecPointwiseDivide(gn->y_work,gn->y,gn->y_work);CHKERRQ(ierr); /* reuse y_work = y./sqrt(y.^2+epsilon^2) */ 94 if (gn->D) { 95 ierr = MatMultTranspose(gn->D,gn->y_work,gn->x_work);CHKERRQ(ierr); 96 } else { 97 ierr = VecCopy(gn->y_work,gn->x_work);CHKERRQ(ierr); 98 } 99 ierr = VecAXPY(G,gn->lambda,gn->x_work);CHKERRQ(ierr); 100 break; 101 } 102 PetscFunctionReturn(0); 103 } 104 105 static PetscErrorCode GNComputeHessian(Tao tao,Vec X,Mat H,Mat Hpre,void *ptr) 106 { 107 TAO_BRGN *gn = (TAO_BRGN *)ptr; 108 PetscErrorCode ierr; 109 110 PetscFunctionBegin; 111 ierr = TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);CHKERRQ(ierr); 112 113 switch (gn->reg_type) { 114 case BRGN_user: 115 ierr = (*gn->regularizerhessian)(tao,X,gn->Hreg,gn->reg_hess_ctx);CHKERRQ(ierr); 116 break; 117 case BRGN_l2prox: 118 break; 119 case BRGN_l1dict: 120 /* calculate and store diagonal matrix as a vector: diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3* --> diag = epsilon^2 ./ sqrt(y.^2+epsilon^2).^3,where y = D*x */ 121 if (gn->D) { 122 ierr = MatMult(gn->D,X,gn->y);CHKERRQ(ierr);/* y = D*x */ 123 } else { 124 ierr = VecCopy(X,gn->y);CHKERRQ(ierr); 125 } 126 ierr = VecPointwiseMult(gn->y_work,gn->y,gn->y);CHKERRQ(ierr); 127 ierr = VecShift(gn->y_work,gn->epsilon*gn->epsilon);CHKERRQ(ierr); 128 ierr = VecCopy(gn->y_work,gn->diag);CHKERRQ(ierr); /* gn->diag = y.^2+epsilon^2 */ 129 ierr = VecSqrtAbs(gn->y_work);CHKERRQ(ierr); /* gn->y_work = sqrt(y.^2+epsilon^2) */ 130 ierr = VecPointwiseMult(gn->diag,gn->y_work,gn->diag);CHKERRQ(ierr);/* gn->diag = sqrt(y.^2+epsilon^2).^3 */ 131 ierr = VecReciprocal(gn->diag);CHKERRQ(ierr); 132 ierr = VecScale(gn->diag,gn->epsilon*gn->epsilon);CHKERRQ(ierr); 133 break; 134 } 135 136 PetscFunctionReturn(0); 137 } 138 139 static PetscErrorCode GNHookFunction(Tao tao,PetscInt iter) 140 { 141 TAO_BRGN *gn = (TAO_BRGN *)tao->user_update; 142 PetscErrorCode ierr; 143 144 PetscFunctionBegin; 145 /* Update basic tao information from the subsolver */ 146 gn->parent->nfuncs = tao->nfuncs; 147 gn->parent->ngrads = tao->ngrads; 148 gn->parent->nfuncgrads = tao->nfuncgrads; 149 gn->parent->nhess = tao->nhess; 150 gn->parent->niter = tao->niter; 151 gn->parent->ksp_its = tao->ksp_its; 152 gn->parent->ksp_tot_its = tao->ksp_tot_its; 153 ierr = TaoGetConvergedReason(tao,&gn->parent->reason);CHKERRQ(ierr); 154 /* Update the solution vectors */ 155 if (iter == 0) { 156 ierr = VecSet(gn->x_old,0.0);CHKERRQ(ierr); 157 } else { 158 ierr = VecCopy(tao->solution,gn->x_old);CHKERRQ(ierr); 159 ierr = VecCopy(tao->solution,gn->parent->solution);CHKERRQ(ierr); 160 } 161 /* Update the gradient */ 162 ierr = VecCopy(tao->gradient,gn->parent->gradient);CHKERRQ(ierr); 163 /* Call general purpose update function */ 164 if (gn->parent->ops->update) { 165 ierr = (*gn->parent->ops->update)(gn->parent,gn->parent->niter);CHKERRQ(ierr); 166 } 167 PetscFunctionReturn(0); 168 } 169 170 static PetscErrorCode TaoSolve_BRGN(Tao tao) 171 { 172 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 173 PetscErrorCode ierr; 174 175 PetscFunctionBegin; 176 ierr = TaoSolve(gn->subsolver);CHKERRQ(ierr); 177 /* Update basic tao information from the subsolver */ 178 tao->nfuncs = gn->subsolver->nfuncs; 179 tao->ngrads = gn->subsolver->ngrads; 180 tao->nfuncgrads = gn->subsolver->nfuncgrads; 181 tao->nhess = gn->subsolver->nhess; 182 tao->niter = gn->subsolver->niter; 183 tao->ksp_its = gn->subsolver->ksp_its; 184 tao->ksp_tot_its = gn->subsolver->ksp_tot_its; 185 ierr = TaoGetConvergedReason(gn->subsolver,&tao->reason);CHKERRQ(ierr); 186 /* Update vectors */ 187 ierr = VecCopy(gn->subsolver->solution,tao->solution);CHKERRQ(ierr); 188 ierr = VecCopy(gn->subsolver->gradient,tao->gradient);CHKERRQ(ierr); 189 PetscFunctionReturn(0); 190 } 191 192 static PetscErrorCode TaoSetFromOptions_BRGN(PetscOptionItems *PetscOptionsObject,Tao tao) 193 { 194 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 195 PetscErrorCode ierr; 196 197 PetscFunctionBegin; 198 ierr = PetscOptionsHead(PetscOptionsObject,"least-squares problems with regularizer: ||f(x)||^2 + lambda*g(x), g(x) = ||xk-xkm1||^2 or ||Dx||_1 or user defined function.");CHKERRQ(ierr); 199 ierr = PetscOptionsReal("-tao_brgn_lambda","regularizer weight (default 1e-4)","",gn->lambda,&gn->lambda,NULL);CHKERRQ(ierr); 200 ierr = PetscOptionsReal("-tao_brgn_epsilon","L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)","",gn->epsilon,&gn->epsilon,NULL);CHKERRQ(ierr); 201 ierr = PetscOptionsEList("-tao_brgn_reg_type","regularization type", "",BRGN_Table,BRGNRegTypes,BRGN_Table[gn->reg_type],&gn->reg_type,NULL);CHKERRQ(ierr); 202 ierr = PetscOptionsTail();CHKERRQ(ierr); 203 ierr = TaoSetFromOptions(gn->subsolver);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 static PetscErrorCode TaoView_BRGN(Tao tao,PetscViewer viewer) 208 { 209 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 210 PetscErrorCode ierr; 211 212 PetscFunctionBegin; 213 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 214 ierr = TaoView(gn->subsolver,viewer);CHKERRQ(ierr); 215 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 216 PetscFunctionReturn(0); 217 } 218 219 static PetscErrorCode TaoSetUp_BRGN(Tao tao) 220 { 221 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 222 PetscErrorCode ierr; 223 PetscBool is_bnls,is_bntr,is_bntl; 224 PetscInt i,n,N,K; /* dict has size K*N*/ 225 226 PetscFunctionBegin; 227 if (!tao->ls_res) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualRoutine() must be called before setup!"); 228 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNLS,&is_bnls);CHKERRQ(ierr); 229 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTR,&is_bntr);CHKERRQ(ierr); 230 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTL,&is_bntl);CHKERRQ(ierr); 231 if ((is_bnls || is_bntr || is_bntl) && !tao->ls_jac) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualJacobianRoutine() must be called before setup!"); 232 if (!tao->gradient) { 233 ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr); 234 } 235 if (!gn->x_work) { 236 ierr = VecDuplicate(tao->solution,&gn->x_work);CHKERRQ(ierr); 237 } 238 if (!gn->r_work) { 239 ierr = VecDuplicate(tao->ls_res,&gn->r_work);CHKERRQ(ierr); 240 } 241 if (!gn->x_old) { 242 ierr = VecDuplicate(tao->solution,&gn->x_old);CHKERRQ(ierr); 243 ierr = VecSet(gn->x_old,0.0);CHKERRQ(ierr); 244 } 245 246 /*ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr);*/ 247 if (BRGN_l1dict == gn->reg_type) { 248 if (gn->D) { 249 /* XH: debug: check matrix */ 250 #if 0 251 ierr = PetscPrintf(PETSC_COMM_SELF,"-------- Check D matrix: -------- \n"); CHKERRQ(ierr); 252 ierr = MatView(gn->D,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 253 #endif 254 ierr = MatGetSize(gn->D,&K,&N);CHKERRQ(ierr); /* Shell matrices still must have sizes defined. K = N for identity matrix, K=N-1 or N for gradient matrix */ 255 } else { 256 ierr = VecGetSize(tao->solution,&K);CHKERRQ(ierr); /* If user does not setup dict matrix, use identiy matrix, K=N */ 257 } 258 if (!gn->y) { 259 ierr = VecCreate(PETSC_COMM_SELF,&gn->y);CHKERRQ(ierr); 260 ierr = VecSetSizes(gn->y,PETSC_DECIDE,K);CHKERRQ(ierr); 261 ierr = VecSetFromOptions(gn->y);CHKERRQ(ierr); 262 ierr = VecSet(gn->y,0.0);CHKERRQ(ierr); 263 264 } 265 if (!gn->y_work) { 266 ierr = VecDuplicate(gn->y,&gn->y_work);CHKERRQ(ierr); 267 } 268 if (!gn->diag) { 269 ierr = VecDuplicate(gn->y,&gn->diag);CHKERRQ(ierr); 270 ierr = VecSet(gn->diag,0.0);CHKERRQ(ierr); 271 } 272 } 273 274 275 if (!tao->setupcalled) { 276 /* Hessian setup */ 277 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 278 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 279 ierr = MatSetSizes(gn->H,n,n,N,N);CHKERRQ(ierr); 280 ierr = MatSetType(gn->H,MATSHELL);CHKERRQ(ierr); 281 ierr = MatSetUp(gn->H);CHKERRQ(ierr); 282 ierr = MatShellSetOperation(gn->H,MATOP_MULT,(void (*)(void))GNHessianProd);CHKERRQ(ierr); 283 ierr = MatShellSetContext(gn->H,(void*)gn);CHKERRQ(ierr); 284 /* Subsolver setup,include initial vector and dicttionary D */ 285 ierr = TaoSetUpdate(gn->subsolver,GNHookFunction,(void*)gn);CHKERRQ(ierr); 286 ierr = TaoSetInitialVector(gn->subsolver,tao->solution);CHKERRQ(ierr); 287 if (tao->bounded) { 288 ierr = TaoSetVariableBounds(gn->subsolver,tao->XL,tao->XU);CHKERRQ(ierr); 289 } 290 ierr = TaoSetResidualRoutine(gn->subsolver,tao->ls_res,tao->ops->computeresidual,tao->user_lsresP);CHKERRQ(ierr); 291 ierr = TaoSetJacobianResidualRoutine(gn->subsolver,tao->ls_jac,tao->ls_jac,tao->ops->computeresidualjacobian,tao->user_lsjacP);CHKERRQ(ierr); 292 ierr = TaoSetObjectiveAndGradientRoutine(gn->subsolver,GNObjectiveGradientEval,(void*)gn);CHKERRQ(ierr); 293 ierr = TaoSetHessianRoutine(gn->subsolver,gn->H,gn->H,GNComputeHessian,(void*)gn);CHKERRQ(ierr); 294 /* Propagate some options down */ 295 ierr = TaoSetTolerances(gn->subsolver,tao->gatol,tao->grtol,tao->gttol);CHKERRQ(ierr); 296 ierr = TaoSetMaximumIterations(gn->subsolver,tao->max_it);CHKERRQ(ierr); 297 ierr = TaoSetMaximumFunctionEvaluations(gn->subsolver,tao->max_funcs);CHKERRQ(ierr); 298 for (i=0; i<tao->numbermonitors; ++i) { 299 ierr = TaoSetMonitor(gn->subsolver,tao->monitor[i],tao->monitorcontext[i],tao->monitordestroy[i]);CHKERRQ(ierr); 300 ierr = PetscObjectReference((PetscObject)(tao->monitorcontext[i]));CHKERRQ(ierr); 301 } 302 ierr = TaoSetUp(gn->subsolver);CHKERRQ(ierr); 303 } 304 PetscFunctionReturn(0); 305 } 306 307 static PetscErrorCode TaoDestroy_BRGN(Tao tao) 308 { 309 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 310 PetscErrorCode ierr; 311 312 PetscFunctionBegin; 313 if (tao->setupcalled) { 314 ierr = VecDestroy(&tao->gradient);CHKERRQ(ierr); 315 ierr = VecDestroy(&gn->x_work);CHKERRQ(ierr); 316 ierr = VecDestroy(&gn->r_work);CHKERRQ(ierr); 317 ierr = VecDestroy(&gn->x_old);CHKERRQ(ierr); 318 ierr = VecDestroy(&gn->diag);CHKERRQ(ierr); 319 ierr = VecDestroy(&gn->y);CHKERRQ(ierr); 320 ierr = VecDestroy(&gn->y_work);CHKERRQ(ierr); 321 } 322 ierr = MatDestroy(&gn->H);CHKERRQ(ierr); 323 ierr = MatDestroy(&gn->D);CHKERRQ(ierr); 324 ierr = TaoDestroy(&gn->subsolver);CHKERRQ(ierr); 325 gn->parent = NULL; 326 ierr = PetscFree(tao->data);CHKERRQ(ierr); 327 PetscFunctionReturn(0); 328 } 329 330 /*MC 331 TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares 332 problems with bound constraints. This algorithm is a thin wrapper around TAOBNTL 333 that constructs the Guass-Newton problem with the user-provided least-squares 334 residual and Jacobian. The problem is regularized with an L2-norm proximal point 335 term. 336 337 Options Database Keys: 338 + -tao_bqnk_max_cg_its - maximum number of bounded conjugate-gradient iterations taken in each Newton loop 339 . -tao_bqnk_init_type - trust radius initialization method ("constant", "direction", "interpolation") 340 . -tao_bqnk_update_type - trust radius update method ("step", "direction", "interpolation") 341 - -tao_bqnk_as_type - active-set estimation method ("none", "bertsekas") 342 343 Level: beginner 344 M*/ 345 PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao) 346 { 347 TAO_BRGN *gn; 348 PetscErrorCode ierr; 349 350 PetscFunctionBegin; 351 ierr = PetscNewLog(tao,&gn);CHKERRQ(ierr); 352 353 tao->ops->destroy = TaoDestroy_BRGN; 354 tao->ops->setup = TaoSetUp_BRGN; 355 tao->ops->setfromoptions = TaoSetFromOptions_BRGN; 356 tao->ops->view = TaoView_BRGN; 357 tao->ops->solve = TaoSolve_BRGN; 358 359 tao->data = (void*)gn; 360 gn->lambda = 1e-4; 361 gn->epsilon = 1e-6; 362 gn->parent = tao; 363 364 ierr = MatCreate(PetscObjectComm((PetscObject)tao),&gn->H);CHKERRQ(ierr); 365 ierr = MatSetOptionsPrefix(gn->H,"tao_brgn_hessian_");CHKERRQ(ierr); 366 367 ierr = TaoCreate(PetscObjectComm((PetscObject)tao),&gn->subsolver);CHKERRQ(ierr); 368 ierr = TaoSetType(gn->subsolver,TAOBNLS);CHKERRQ(ierr); 369 ierr = TaoSetOptionsPrefix(gn->subsolver,"tao_brgn_subsolver_");CHKERRQ(ierr); 370 PetscFunctionReturn(0); 371 } 372 373 /*@C 374 TaoBRGNGetSubsolver - Get the pointer to the subsolver inside BRGN 375 376 Collective on Tao 377 378 Level: developer 379 380 Input Parameters: 381 + tao - the Tao solver context 382 - subsolver - the Tao sub-solver context 383 @*/ 384 PetscErrorCode TaoBRGNGetSubsolver(Tao tao,Tao *subsolver) 385 { 386 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 387 388 PetscFunctionBegin; 389 *subsolver = gn->subsolver; 390 PetscFunctionReturn(0); 391 } 392 393 /*@C 394 TaoBRGNSetL1RegularizerWeight - Set the L1-norm regularizer weight for the Gauss-Newton least-squares algorithm 395 396 Collective on Tao 397 398 Level: developer 399 400 Input Parameters: 401 + tao - the Tao solver context 402 - lambda - L1-norm regularizer weight 403 @*/ 404 PetscErrorCode TaoBRGNSetRegularizerWeight(Tao tao,PetscReal lambda) 405 { 406 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 407 408 /* Initialize lambda here */ 409 410 PetscFunctionBegin; 411 gn->lambda = lambda; 412 PetscFunctionReturn(0); 413 } 414 415 /*@C 416 TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm 417 418 Collective on Tao 419 420 Level: developer 421 422 Input Parameters: 423 + tao - the Tao solver context 424 - epsilon - L1-norm smooth approximation parameter 425 @*/ 426 PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao,PetscReal epsilon) 427 { 428 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 429 430 /* Initialize epsilon here */ 431 432 PetscFunctionBegin; 433 gn->epsilon = epsilon; 434 PetscFunctionReturn(0); 435 } 436 437 /*@C 438 TaoBRGNSetDictionaryMatrix - bind the dictionary matrix from user application context to gn->D, for compressed sensing (with least-squares problem) 439 440 Input Parameters: 441 + tao - the Tao context 442 . dict - the user specified dictionary matrix 443 444 Level: developer 445 @*/ 446 PetscErrorCode TaoBRGNSetDictionaryMatrix(Tao tao,Mat dict) 447 { 448 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 449 PetscErrorCode ierr; 450 PetscFunctionBegin; 451 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 452 if (dict) { 453 PetscValidHeaderSpecific(dict,MAT_CLASSID,2); 454 PetscCheckSameComm(tao,1,dict,2); 455 ierr = PetscObjectReference((PetscObject)dict);CHKERRQ(ierr); 456 } 457 ierr = MatDestroy(&gn->D);CHKERRQ(ierr); 458 gn->D = dict; /* We allow to set a null dictionary, which means we just use default identity matrix? */ 459 PetscFunctionReturn(0); 460 } 461 462 /*@C 463 @*/ 464 PetscErrorCode TaoBRGNSetRegularizerObjectiveAndGradientRoutine(Tao tao,PetscErrorCode (*func)(Tao,Vec,PetscReal *,Vec,void*),void *ctx) 465 { 466 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 467 468 PetscFunctionBegin; 469 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 470 if (ctx) { 471 gn->reg_obj_ctx = ctx; 472 } 473 if (func) { 474 gn->regularizerobjandgrad = func; 475 } 476 PetscFunctionReturn(0); 477 } 478 479 /*@C 480 @*/ 481 PetscErrorCode TaoBRGNSetRegularizerHessianRoutine(Tao tao,Mat Hreg,PetscErrorCode (*func)(Tao,Vec,Mat,void*),void *ctx) 482 { 483 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 484 PetscErrorCode ierr; 485 486 PetscFunctionBegin; 487 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 488 if (Hreg) { 489 PetscValidHeaderSpecific(Hreg,MAT_CLASSID,2); 490 PetscCheckSameComm(tao,1,Hreg,2); 491 } else { 492 SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONG,"NULL Hessian detected! User must provide valid Hessian for the regularizer."); 493 } 494 if (ctx) { 495 gn->reg_hess_ctx = ctx; 496 } 497 if (func) { 498 gn->regularizerhessian = func; 499 } 500 if (Hreg) { 501 ierr = PetscObjectReference((PetscObject)Hreg);CHKERRQ(ierr); 502 ierr = MatDestroy(&gn->Hreg);CHKERRQ(ierr); 503 gn->Hreg = Hreg; 504 } 505 PetscFunctionReturn(0); 506 }