Actual source code: brgn.c

  1: #include <../src/tao/leastsquares/impls/brgn/brgn.h>

  3: #define BRGN_REGULARIZATION_USER   0
  4: #define BRGN_REGULARIZATION_L2PROX 1
  5: #define BRGN_REGULARIZATION_L2PURE 2
  6: #define BRGN_REGULARIZATION_L1DICT 3
  7: #define BRGN_REGULARIZATION_LM     4
  8: #define BRGN_REGULARIZATION_TYPES  5

 10: static const char *BRGN_REGULARIZATION_TABLE[64] = {"user", "l2prox", "l2pure", "l1dict", "lm"};

 12: static PetscErrorCode GNHessianProd(Mat H, Vec in, Vec out)
 13: {
 14:   TAO_BRGN *gn;

 16:   MatShellGetContext(H, &gn);
 17:   MatMult(gn->subsolver->ls_jac, in, gn->r_work);
 18:   MatMultTranspose(gn->subsolver->ls_jac, gn->r_work, out);
 19:   switch (gn->reg_type) {
 20:   case BRGN_REGULARIZATION_USER:
 21:     MatMult(gn->Hreg, in, gn->x_work);
 22:     VecAXPY(out, gn->lambda, gn->x_work);
 23:     break;
 24:   case BRGN_REGULARIZATION_L2PURE:
 25:     VecAXPY(out, gn->lambda, in);
 26:     break;
 27:   case BRGN_REGULARIZATION_L2PROX:
 28:     VecAXPY(out, gn->lambda, in);
 29:     break;
 30:   case BRGN_REGULARIZATION_L1DICT:
 31:     /* out = out + lambda*D'*(diag.*(D*in)) */
 32:     if (gn->D) {
 33:       MatMult(gn->D, in, gn->y); /* y = D*in */
 34:     } else {
 35:       VecCopy(in, gn->y);
 36:     }
 37:     VecPointwiseMult(gn->y_work, gn->diag, gn->y); /* y_work = diag.*(D*in), where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */
 38:     if (gn->D) {
 39:       MatMultTranspose(gn->D, gn->y_work, gn->x_work); /* x_work = D'*(diag.*(D*in)) */
 40:     } else {
 41:       VecCopy(gn->y_work, gn->x_work);
 42:     }
 43:     VecAXPY(out, gn->lambda, gn->x_work);
 44:     break;
 45:   case BRGN_REGULARIZATION_LM:
 46:     VecPointwiseMult(gn->x_work, gn->damping, in);
 47:     VecAXPY(out, 1, gn->x_work);
 48:     break;
 49:   }
 50:   return 0;
 51: }
 52: static PetscErrorCode ComputeDamping(TAO_BRGN *gn)
 53: {
 54:   const PetscScalar *diag_ary;
 55:   PetscScalar       *damping_ary;
 56:   PetscInt           i, n;

 58:   /* update damping */
 59:   VecGetArray(gn->damping, &damping_ary);
 60:   VecGetArrayRead(gn->diag, &diag_ary);
 61:   VecGetLocalSize(gn->damping, &n);
 62:   for (i = 0; i < n; i++) damping_ary[i] = PetscClipInterval(diag_ary[i], PETSC_SQRT_MACHINE_EPSILON, PetscSqrtReal(PETSC_MAX_REAL));
 63:   VecScale(gn->damping, gn->lambda);
 64:   VecRestoreArray(gn->damping, &damping_ary);
 65:   VecRestoreArrayRead(gn->diag, &diag_ary);
 66:   return 0;
 67: }

 69: PetscErrorCode TaoBRGNGetDampingVector(Tao tao, Vec *d)
 70: {
 71:   TAO_BRGN *gn = (TAO_BRGN *)tao->data;

 74:   *d = gn->damping;
 75:   return 0;
 76: }

 78: static PetscErrorCode GNObjectiveGradientEval(Tao tao, Vec X, PetscReal *fcn, Vec G, void *ptr)
 79: {
 80:   TAO_BRGN   *gn = (TAO_BRGN *)ptr;
 81:   PetscInt    K; /* dimension of D*X */
 82:   PetscScalar yESum;
 83:   PetscReal   f_reg;

 85:   /* compute objective *fcn*/
 86:   /* compute first term 0.5*||ls_res||_2^2 */
 87:   TaoComputeResidual(tao, X, tao->ls_res);
 88:   VecDot(tao->ls_res, tao->ls_res, fcn);
 89:   *fcn *= 0.5;
 90:   /* compute gradient G */
 91:   TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre);
 92:   MatMultTranspose(tao->ls_jac, tao->ls_res, G);
 93:   /* add the regularization contribution */
 94:   switch (gn->reg_type) {
 95:   case BRGN_REGULARIZATION_USER:
 96:     (*gn->regularizerobjandgrad)(tao, X, &f_reg, gn->x_work, gn->reg_obj_ctx);
 97:     *fcn += gn->lambda * f_reg;
 98:     VecAXPY(G, gn->lambda, gn->x_work);
 99:     break;
100:   case BRGN_REGULARIZATION_L2PURE:
101:     /* compute f = f + lambda*0.5*xk'*xk */
102:     VecDot(X, X, &f_reg);
103:     *fcn += gn->lambda * 0.5 * f_reg;
104:     /* compute G = G + lambda*xk */
105:     VecAXPY(G, gn->lambda, X);
106:     break;
107:   case BRGN_REGULARIZATION_L2PROX:
108:     /* compute f = f + lambda*0.5*(xk - xkm1)'*(xk - xkm1) */
109:     VecAXPBYPCZ(gn->x_work, 1.0, -1.0, 0.0, X, gn->x_old);
110:     VecDot(gn->x_work, gn->x_work, &f_reg);
111:     *fcn += gn->lambda * 0.5 * f_reg;
112:     /* compute G = G + lambda*(xk - xkm1) */
113:     VecAXPBYPCZ(G, gn->lambda, -gn->lambda, 1.0, X, gn->x_old);
114:     break;
115:   case BRGN_REGULARIZATION_L1DICT:
116:     /* compute f = f + lambda*sum(sqrt(y.^2+epsilon^2) - epsilon), where y = D*x*/
117:     if (gn->D) {
118:       MatMult(gn->D, X, gn->y); /* y = D*x */
119:     } else {
120:       VecCopy(X, gn->y);
121:     }
122:     VecPointwiseMult(gn->y_work, gn->y, gn->y);
123:     VecShift(gn->y_work, gn->epsilon * gn->epsilon);
124:     VecSqrtAbs(gn->y_work); /* gn->y_work = sqrt(y.^2+epsilon^2) */
125:     VecSum(gn->y_work, &yESum);
126:     VecGetSize(gn->y, &K);
127:     *fcn += gn->lambda * (yESum - K * gn->epsilon);
128:     /* compute G = G + lambda*D'*(y./sqrt(y.^2+epsilon^2)),where y = D*x */
129:     VecPointwiseDivide(gn->y_work, gn->y, gn->y_work); /* reuse y_work = y./sqrt(y.^2+epsilon^2) */
130:     if (gn->D) {
131:       MatMultTranspose(gn->D, gn->y_work, gn->x_work);
132:     } else {
133:       VecCopy(gn->y_work, gn->x_work);
134:     }
135:     VecAXPY(G, gn->lambda, gn->x_work);
136:     break;
137:   }
138:   return 0;
139: }

141: static PetscErrorCode GNComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
142: {
143:   TAO_BRGN    *gn = (TAO_BRGN *)ptr;
144:   PetscInt     i, n, cstart, cend;
145:   PetscScalar *cnorms, *diag_ary;

147:   TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre);
148:   if (gn->mat_explicit) MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_REUSE_MATRIX, PETSC_DEFAULT, &gn->H);

150:   switch (gn->reg_type) {
151:   case BRGN_REGULARIZATION_USER:
152:     (*gn->regularizerhessian)(tao, X, gn->Hreg, gn->reg_hess_ctx);
153:     if (gn->mat_explicit) MatAXPY(gn->H, 1.0, gn->Hreg, DIFFERENT_NONZERO_PATTERN);
154:     break;
155:   case BRGN_REGULARIZATION_L2PURE:
156:     if (gn->mat_explicit) MatShift(gn->H, gn->lambda);
157:     break;
158:   case BRGN_REGULARIZATION_L2PROX:
159:     if (gn->mat_explicit) MatShift(gn->H, gn->lambda);
160:     break;
161:   case BRGN_REGULARIZATION_L1DICT:
162:     /* 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 */
163:     if (gn->D) {
164:       MatMult(gn->D, X, gn->y); /* y = D*x */
165:     } else {
166:       VecCopy(X, gn->y);
167:     }
168:     VecPointwiseMult(gn->y_work, gn->y, gn->y);
169:     VecShift(gn->y_work, gn->epsilon * gn->epsilon);
170:     VecCopy(gn->y_work, gn->diag);                    /* gn->diag = y.^2+epsilon^2 */
171:     VecSqrtAbs(gn->y_work);                           /* gn->y_work = sqrt(y.^2+epsilon^2) */
172:     VecPointwiseMult(gn->diag, gn->y_work, gn->diag); /* gn->diag = sqrt(y.^2+epsilon^2).^3 */
173:     VecReciprocal(gn->diag);
174:     VecScale(gn->diag, gn->epsilon * gn->epsilon);
175:     if (gn->mat_explicit) MatDiagonalSet(gn->H, gn->diag, ADD_VALUES);
176:     break;
177:   case BRGN_REGULARIZATION_LM:
178:     /* compute diagonal of J^T J */
179:     MatGetSize(gn->parent->ls_jac, NULL, &n);
180:     PetscMalloc1(n, &cnorms);
181:     MatGetColumnNorms(gn->parent->ls_jac, NORM_2, cnorms);
182:     MatGetOwnershipRangeColumn(gn->parent->ls_jac, &cstart, &cend);
183:     VecGetArray(gn->diag, &diag_ary);
184:     for (i = 0; i < cend - cstart; i++) diag_ary[i] = cnorms[cstart + i] * cnorms[cstart + i];
185:     VecRestoreArray(gn->diag, &diag_ary);
186:     PetscFree(cnorms);
187:     ComputeDamping(gn);
188:     if (gn->mat_explicit) MatDiagonalSet(gn->H, gn->damping, ADD_VALUES);
189:     break;
190:   }
191:   return 0;
192: }

194: static PetscErrorCode GNHookFunction(Tao tao, PetscInt iter, void *ctx)
195: {
196:   TAO_BRGN *gn = (TAO_BRGN *)ctx;

198:   /* Update basic tao information from the subsolver */
199:   gn->parent->nfuncs      = tao->nfuncs;
200:   gn->parent->ngrads      = tao->ngrads;
201:   gn->parent->nfuncgrads  = tao->nfuncgrads;
202:   gn->parent->nhess       = tao->nhess;
203:   gn->parent->niter       = tao->niter;
204:   gn->parent->ksp_its     = tao->ksp_its;
205:   gn->parent->ksp_tot_its = tao->ksp_tot_its;
206:   gn->parent->fc          = tao->fc;
207:   TaoGetConvergedReason(tao, &gn->parent->reason);
208:   /* Update the solution vectors */
209:   if (iter == 0) {
210:     VecSet(gn->x_old, 0.0);
211:   } else {
212:     VecCopy(tao->solution, gn->x_old);
213:     VecCopy(tao->solution, gn->parent->solution);
214:   }
215:   /* Update the gradient */
216:   VecCopy(tao->gradient, gn->parent->gradient);

218:   /* Update damping parameter for LM */
219:   if (gn->reg_type == BRGN_REGULARIZATION_LM) {
220:     if (iter > 0) {
221:       if (gn->fc_old > tao->fc) {
222:         gn->lambda = gn->lambda * gn->downhill_lambda_change;
223:       } else {
224:         /* uphill step */
225:         gn->lambda = gn->lambda * gn->uphill_lambda_change;
226:       }
227:     }
228:     gn->fc_old = tao->fc;
229:   }

231:   /* Call general purpose update function */
232:   if (gn->parent->ops->update) (*gn->parent->ops->update)(gn->parent, gn->parent->niter, gn->parent->user_update);
233:   return 0;
234: }

236: static PetscErrorCode TaoSolve_BRGN(Tao tao)
237: {
238:   TAO_BRGN *gn = (TAO_BRGN *)tao->data;

240:   TaoSolve(gn->subsolver);
241:   /* Update basic tao information from the subsolver */
242:   tao->nfuncs      = gn->subsolver->nfuncs;
243:   tao->ngrads      = gn->subsolver->ngrads;
244:   tao->nfuncgrads  = gn->subsolver->nfuncgrads;
245:   tao->nhess       = gn->subsolver->nhess;
246:   tao->niter       = gn->subsolver->niter;
247:   tao->ksp_its     = gn->subsolver->ksp_its;
248:   tao->ksp_tot_its = gn->subsolver->ksp_tot_its;
249:   TaoGetConvergedReason(gn->subsolver, &tao->reason);
250:   /* Update vectors */
251:   VecCopy(gn->subsolver->solution, tao->solution);
252:   VecCopy(gn->subsolver->gradient, tao->gradient);
253:   return 0;
254: }

256: static PetscErrorCode TaoSetFromOptions_BRGN(Tao tao, PetscOptionItems *PetscOptionsObject)
257: {
258:   TAO_BRGN     *gn = (TAO_BRGN *)tao->data;
259:   TaoLineSearch ls;

261:   PetscOptionsHeadBegin(PetscOptionsObject, "least-squares problems with regularizer: ||f(x)||^2 + lambda*g(x), g(x) = ||xk-xkm1||^2 or ||Dx||_1 or user defined function.");
262:   PetscOptionsBool("-tao_brgn_mat_explicit", "switches the Hessian construction to be an explicit matrix rather than MATSHELL", "", gn->mat_explicit, &gn->mat_explicit, NULL);
263:   PetscOptionsReal("-tao_brgn_regularizer_weight", "regularizer weight (default 1e-4)", "", gn->lambda, &gn->lambda, NULL);
264:   PetscOptionsReal("-tao_brgn_l1_smooth_epsilon", "L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)", "", gn->epsilon, &gn->epsilon, NULL);
265:   PetscOptionsReal("-tao_brgn_lm_downhill_lambda_change", "Factor to decrease trust region by on downhill steps", "", gn->downhill_lambda_change, &gn->downhill_lambda_change, NULL);
266:   PetscOptionsReal("-tao_brgn_lm_uphill_lambda_change", "Factor to increase trust region by on uphill steps", "", gn->uphill_lambda_change, &gn->uphill_lambda_change, NULL);
267:   PetscOptionsEList("-tao_brgn_regularization_type", "regularization type", "", BRGN_REGULARIZATION_TABLE, BRGN_REGULARIZATION_TYPES, BRGN_REGULARIZATION_TABLE[gn->reg_type], &gn->reg_type, NULL);
268:   PetscOptionsHeadEnd();
269:   /* set unit line search direction as the default when using the lm regularizer */
270:   if (gn->reg_type == BRGN_REGULARIZATION_LM) {
271:     TaoGetLineSearch(gn->subsolver, &ls);
272:     TaoLineSearchSetType(ls, TAOLINESEARCHUNIT);
273:   }
274:   TaoSetFromOptions(gn->subsolver);
275:   return 0;
276: }

278: static PetscErrorCode TaoView_BRGN(Tao tao, PetscViewer viewer)
279: {
280:   TAO_BRGN *gn = (TAO_BRGN *)tao->data;

282:   PetscViewerASCIIPushTab(viewer);
283:   TaoView(gn->subsolver, viewer);
284:   PetscViewerASCIIPopTab(viewer);
285:   return 0;
286: }

288: static PetscErrorCode TaoSetUp_BRGN(Tao tao)
289: {
290:   TAO_BRGN *gn = (TAO_BRGN *)tao->data;
291:   PetscBool is_bnls, is_bntr, is_bntl;
292:   PetscInt  i, n, N, K; /* dict has size K*N*/

295:   PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNLS, &is_bnls);
296:   PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTR, &is_bntr);
297:   PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTL, &is_bntl);
299:   if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
300:   if (!gn->x_work) VecDuplicate(tao->solution, &gn->x_work);
301:   if (!gn->r_work) VecDuplicate(tao->ls_res, &gn->r_work);
302:   if (!gn->x_old) {
303:     VecDuplicate(tao->solution, &gn->x_old);
304:     VecSet(gn->x_old, 0.0);
305:   }

307:   if (BRGN_REGULARIZATION_L1DICT == gn->reg_type) {
308:     if (!gn->y) {
309:       if (gn->D) {
310:         MatGetSize(gn->D, &K, &N); /* Shell matrices still must have sizes defined. K = N for identity matrix, K=N-1 or N for gradient matrix */
311:         MatCreateVecs(gn->D, NULL, &gn->y);
312:       } else {
313:         VecDuplicate(tao->solution, &gn->y); /* If user does not setup dict matrix, use identiy matrix, K=N */
314:       }
315:       VecSet(gn->y, 0.0);
316:     }
317:     if (!gn->y_work) VecDuplicate(gn->y, &gn->y_work);
318:     if (!gn->diag) {
319:       VecDuplicate(gn->y, &gn->diag);
320:       VecSet(gn->diag, 0.0);
321:     }
322:   }
323:   if (BRGN_REGULARIZATION_LM == gn->reg_type) {
324:     if (!gn->diag) MatCreateVecs(tao->ls_jac, &gn->diag, NULL);
325:     if (!gn->damping) MatCreateVecs(tao->ls_jac, &gn->damping, NULL);
326:   }

328:   if (!tao->setupcalled) {
329:     /* Hessian setup */
330:     if (gn->mat_explicit) {
331:       TaoComputeResidualJacobian(tao, tao->solution, tao->ls_jac, tao->ls_jac_pre);
332:       MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &gn->H);
333:     } else {
334:       VecGetLocalSize(tao->solution, &n);
335:       VecGetSize(tao->solution, &N);
336:       MatCreate(PetscObjectComm((PetscObject)tao), &gn->H);
337:       MatSetSizes(gn->H, n, n, N, N);
338:       MatSetType(gn->H, MATSHELL);
339:       MatSetOption(gn->H, MAT_SYMMETRIC, PETSC_TRUE);
340:       MatShellSetOperation(gn->H, MATOP_MULT, (void (*)(void))GNHessianProd);
341:       MatShellSetContext(gn->H, gn);
342:     }
343:     MatSetUp(gn->H);
344:     /* Subsolver setup,include initial vector and dictionary D */
345:     TaoSetUpdate(gn->subsolver, GNHookFunction, gn);
346:     TaoSetSolution(gn->subsolver, tao->solution);
347:     if (tao->bounded) TaoSetVariableBounds(gn->subsolver, tao->XL, tao->XU);
348:     TaoSetResidualRoutine(gn->subsolver, tao->ls_res, tao->ops->computeresidual, tao->user_lsresP);
349:     TaoSetJacobianResidualRoutine(gn->subsolver, tao->ls_jac, tao->ls_jac, tao->ops->computeresidualjacobian, tao->user_lsjacP);
350:     TaoSetObjectiveAndGradient(gn->subsolver, NULL, GNObjectiveGradientEval, gn);
351:     TaoSetHessian(gn->subsolver, gn->H, gn->H, GNComputeHessian, gn);
352:     /* Propagate some options down */
353:     TaoSetTolerances(gn->subsolver, tao->gatol, tao->grtol, tao->gttol);
354:     TaoSetMaximumIterations(gn->subsolver, tao->max_it);
355:     TaoSetMaximumFunctionEvaluations(gn->subsolver, tao->max_funcs);
356:     for (i = 0; i < tao->numbermonitors; ++i) {
357:       TaoSetMonitor(gn->subsolver, tao->monitor[i], tao->monitorcontext[i], tao->monitordestroy[i]);
358:       PetscObjectReference((PetscObject)(tao->monitorcontext[i]));
359:     }
360:     TaoSetUp(gn->subsolver);
361:   }
362:   return 0;
363: }

365: static PetscErrorCode TaoDestroy_BRGN(Tao tao)
366: {
367:   TAO_BRGN *gn = (TAO_BRGN *)tao->data;

369:   if (tao->setupcalled) {
370:     VecDestroy(&tao->gradient);
371:     VecDestroy(&gn->x_work);
372:     VecDestroy(&gn->r_work);
373:     VecDestroy(&gn->x_old);
374:     VecDestroy(&gn->diag);
375:     VecDestroy(&gn->y);
376:     VecDestroy(&gn->y_work);
377:   }
378:   VecDestroy(&gn->damping);
379:   VecDestroy(&gn->diag);
380:   MatDestroy(&gn->H);
381:   MatDestroy(&gn->D);
382:   MatDestroy(&gn->Hreg);
383:   TaoDestroy(&gn->subsolver);
384:   gn->parent = NULL;
385:   PetscFree(tao->data);
386:   return 0;
387: }

389: /*MC
390:   TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares
391:             problems with bound constraints. This algorithm is a thin wrapper around TAOBNTL
392:             that constructs the Gauss-Newton problem with the user-provided least-squares
393:             residual and Jacobian. The algorithm offers an L2-norm ("l2pure"), L2-norm proximal point ("l2prox")
394:             regularizer, and L1-norm dictionary regularizer ("l1dict"), where we approximate the
395:             L1-norm ||x||_1 by sum_i(sqrt(x_i^2+epsilon^2)-epsilon) with a small positive number epsilon.
396:             Also offered is the "lm" regularizer which uses a scaled diagonal of J^T J.
397:             With the "lm" regularizer, BRGN is a Levenberg-Marquardt optimizer.
398:             The user can also provide own regularization function.

400:   Options Database Keys:
401: + -tao_brgn_regularization_type - regularization type ("user", "l2prox", "l2pure", "l1dict", "lm") (default "l2prox")
402: . -tao_brgn_regularizer_weight  - regularizer weight (default 1e-4)
403: - -tao_brgn_l1_smooth_epsilon   - L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)

405:   Level: beginner
406: M*/
407: PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao)
408: {
409:   TAO_BRGN *gn;

411:   PetscNew(&gn);

413:   tao->ops->destroy        = TaoDestroy_BRGN;
414:   tao->ops->setup          = TaoSetUp_BRGN;
415:   tao->ops->setfromoptions = TaoSetFromOptions_BRGN;
416:   tao->ops->view           = TaoView_BRGN;
417:   tao->ops->solve          = TaoSolve_BRGN;

419:   tao->data                  = gn;
420:   gn->reg_type               = BRGN_REGULARIZATION_L2PROX;
421:   gn->lambda                 = 1e-4;
422:   gn->epsilon                = 1e-6;
423:   gn->downhill_lambda_change = 1. / 5.;
424:   gn->uphill_lambda_change   = 1.5;
425:   gn->parent                 = tao;

427:   TaoCreate(PetscObjectComm((PetscObject)tao), &gn->subsolver);
428:   TaoSetType(gn->subsolver, TAOBNLS);
429:   TaoSetOptionsPrefix(gn->subsolver, "tao_brgn_subsolver_");
430:   return 0;
431: }

433: /*@
434:   TaoBRGNGetSubsolver - Get the pointer to the subsolver inside BRGN

436:   Collective

438:   Level: advanced

440:   Input Parameters:
441: +  tao - the Tao solver context
442: -  subsolver - the Tao sub-solver context
443: @*/
444: PetscErrorCode TaoBRGNGetSubsolver(Tao tao, Tao *subsolver)
445: {
446:   TAO_BRGN *gn = (TAO_BRGN *)tao->data;

448:   *subsolver = gn->subsolver;
449:   return 0;
450: }

452: /*@
453:   TaoBRGNSetRegularizerWeight - Set the regularizer weight for the Gauss-Newton least-squares algorithm

455:   Collective

457:   Input Parameters:
458: +  tao - the Tao solver context
459: -  lambda - L1-norm regularizer weight

461:   Level: beginner
462: @*/
463: PetscErrorCode TaoBRGNSetRegularizerWeight(Tao tao, PetscReal lambda)
464: {
465:   TAO_BRGN *gn = (TAO_BRGN *)tao->data;

467:   /* Initialize lambda here */

469:   gn->lambda = lambda;
470:   return 0;
471: }

473: /*@
474:   TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm

476:   Collective

478:   Input Parameters:
479: +  tao - the Tao solver context
480: -  epsilon - L1-norm smooth approximation parameter

482:   Level: advanced
483: @*/
484: PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao, PetscReal epsilon)
485: {
486:   TAO_BRGN *gn = (TAO_BRGN *)tao->data;

488:   /* Initialize epsilon here */

490:   gn->epsilon = epsilon;
491:   return 0;
492: }

494: /*@
495:    TaoBRGNSetDictionaryMatrix - bind the dictionary matrix from user application context to gn->D, for compressed sensing (with least-squares problem)

497:    Input Parameters:
498: +  tao  - the Tao context
499: -  dict - the user specified dictionary matrix.  We allow to set a null dictionary, which means identity matrix by default

501:     Level: advanced
502: @*/
503: PetscErrorCode TaoBRGNSetDictionaryMatrix(Tao tao, Mat dict)
504: {
505:   TAO_BRGN *gn = (TAO_BRGN *)tao->data;
507:   if (dict) {
510:     PetscObjectReference((PetscObject)dict);
511:   }
512:   MatDestroy(&gn->D);
513:   gn->D = dict;
514:   return 0;
515: }

517: /*@C
518:    TaoBRGNSetRegularizerObjectiveAndGradientRoutine - Sets the user-defined regularizer call-back
519:    function into the algorithm.

521:    Input Parameters:
522: + tao - the Tao context
523: . func - function pointer for the regularizer value and gradient evaluation
524: - ctx - user context for the regularizer

526:    Level: advanced
527: @*/
528: PetscErrorCode TaoBRGNSetRegularizerObjectiveAndGradientRoutine(Tao tao, PetscErrorCode (*func)(Tao, Vec, PetscReal *, Vec, void *), void *ctx)
529: {
530:   TAO_BRGN *gn = (TAO_BRGN *)tao->data;

533:   if (ctx) gn->reg_obj_ctx = ctx;
534:   if (func) gn->regularizerobjandgrad = func;
535:   return 0;
536: }

538: /*@C
539:    TaoBRGNSetRegularizerHessianRoutine - Sets the user-defined regularizer call-back
540:    function into the algorithm.

542:    Input Parameters:
543: + tao - the Tao context
544: . Hreg - user-created matrix for the Hessian of the regularization term
545: . func - function pointer for the regularizer Hessian evaluation
546: - ctx - user context for the regularizer Hessian

548:    Level: advanced
549: @*/
550: PetscErrorCode TaoBRGNSetRegularizerHessianRoutine(Tao tao, Mat Hreg, PetscErrorCode (*func)(Tao, Vec, Mat, void *), void *ctx)
551: {
552:   TAO_BRGN *gn = (TAO_BRGN *)tao->data;

555:   if (Hreg) {
558:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "NULL Hessian detected! User must provide valid Hessian for the regularizer.");
559:   if (ctx) gn->reg_hess_ctx = ctx;
560:   if (func) gn->regularizerhessian = func;
561:   if (Hreg) {
562:     PetscObjectReference((PetscObject)Hreg);
563:     MatDestroy(&gn->Hreg);
564:     gn->Hreg = Hreg;
565:   }
566:   return 0;
567: }