Actual source code: admm.c
1: #include <../src/tao/constrained/impls/admm/admm.h>
2: #include <petsctao.h>
3: #include <petsc/private/petscimpl.h>
5: /* Updates terminating criteria
6: *
7: * 1 ||r_k|| = ||Ax+Bz-c|| =< catol_admm* max{||Ax||,||Bz||,||c||}
8: *
9: * 2. Updates dual residual, d_k
10: *
11: * 3. ||d_k|| = ||mu*A^T*B(z_k-z_{k-1})|| =< gatol_admm * ||A^Ty|| */
13: static PetscBool cited = PETSC_FALSE;
14: static const char citation[] =
15: "@misc{xu2017adaptive,\n"
16: " title={Adaptive Relaxed ADMM: Convergence Theory and Practical Implementation},\n"
17: " author={Zheng Xu and Mario A. T. Figueiredo and Xiaoming Yuan and Christoph Studer and Tom Goldstein},\n"
18: " year={2017},\n"
19: " eprint={1704.02712},\n"
20: " archivePrefix={arXiv},\n"
21: " primaryClass={cs.CV}\n"
22: "} \n";
24: static PetscErrorCode TaoADMMToleranceUpdate(Tao tao)
25: {
26: TAO_ADMM *am = (TAO_ADMM*)tao->data;
28: PetscReal Axnorm,Bznorm,ATynorm,temp;
29: Vec tempJR,tempL;
30: Tao mis;
33: mis = am->subsolverX;
34: tempJR = am->workJacobianRight;
35: tempL = am->workLeft;
36: /* ATy */
37: TaoComputeJacobianEquality(mis, am->y, mis->jacobian_equality, mis->jacobian_equality_pre);
38: MatMultTranspose(mis->jacobian_equality,am->y,tempJR);
39: VecNorm(tempJR,NORM_2,&ATynorm);
40: /* dualres = mu * ||AT(Bz-Bzold)||_2 */
41: VecWAXPY(tempJR,-1.,am->Bzold,am->Bz);
42: MatMultTranspose(mis->jacobian_equality,tempJR,tempL);
43: VecNorm(tempL,NORM_2,&am->dualres);
44: am->dualres *= am->mu;
46: /* ||Ax||_2, ||Bz||_2 */
47: VecNorm(am->Ax,NORM_2,&Axnorm);
48: VecNorm(am->Bz,NORM_2,&Bznorm);
50: /* Set catol to be catol_admm * max{||Ax||,||Bz||,||c||} *
51: * Set gatol to be gatol_admm * ||A^Ty|| *
52: * while cnorm is ||r_k||_2, and gnorm is ||d_k||_2 */
53: temp = am->catol_admm * PetscMax(Axnorm, (!am->const_norm) ? Bznorm : PetscMax(Bznorm,am->const_norm));
54: TaoSetConstraintTolerances(tao,temp,PETSC_DEFAULT);
55: TaoSetTolerances(tao, am->gatol_admm*ATynorm, PETSC_DEFAULT,PETSC_DEFAULT);
56: return(0);
57: }
59: /* Penaly Update for Adaptive ADMM. */
60: static PetscErrorCode AdaptiveADMMPenaltyUpdate(Tao tao)
61: {
62: TAO_ADMM *am = (TAO_ADMM*)tao->data;
64: PetscReal ydiff_norm, yhatdiff_norm, Axdiff_norm, Bzdiff_norm, Axyhat, Bzy, a_sd, a_mg, a_k, b_sd, b_mg, b_k;
65: PetscBool hflag, gflag;
66: Vec tempJR,tempJR2;
69: tempJR = am->workJacobianRight;
70: tempJR2 = am->workJacobianRight2;
71: hflag = PETSC_FALSE;
72: gflag = PETSC_FALSE;
73: a_k = -1;
74: b_k = -1;
76: VecWAXPY(tempJR,-1.,am->Axold,am->Ax);
77: VecWAXPY(tempJR2,-1.,am->yhatold,am->yhat);
78: VecNorm(tempJR,NORM_2,&Axdiff_norm);
79: VecNorm(tempJR2,NORM_2,&yhatdiff_norm);
80: VecDot(tempJR,tempJR2,&Axyhat);
82: VecWAXPY(tempJR,-1.,am->Bz0,am->Bz);
83: VecWAXPY(tempJR2,-1.,am->y,am->y0);
84: VecNorm(tempJR,NORM_2,&Bzdiff_norm);
85: VecNorm(tempJR2,NORM_2,&ydiff_norm);
86: VecDot(tempJR,tempJR2,&Bzy);
88: if (Axyhat > am->orthval*Axdiff_norm*yhatdiff_norm + am->mueps) {
89: hflag = PETSC_TRUE;
90: a_sd = PetscSqr(yhatdiff_norm)/Axyhat; /* alphaSD */
91: a_mg = Axyhat/PetscSqr(Axdiff_norm); /* alphaMG */
92: a_k = (a_mg/a_sd) > 0.5 ? a_mg : a_sd - 0.5*a_mg;
93: }
94: if (Bzy > am->orthval*Bzdiff_norm*ydiff_norm + am->mueps) {
95: gflag = PETSC_TRUE;
96: b_sd = PetscSqr(ydiff_norm)/Bzy; /* betaSD */
97: b_mg = Bzy/PetscSqr(Bzdiff_norm); /* betaMG */
98: b_k = (b_mg/b_sd) > 0.5 ? b_mg : b_sd - 0.5*b_mg;
99: }
100: am->muold = am->mu;
101: if (gflag && hflag) {
102: am->mu = PetscSqrtReal(a_k*b_k);
103: } else if (hflag) {
104: am->mu = a_k;
105: } else if (gflag) {
106: am->mu = b_k;
107: }
108: if (am->mu > am->muold) {
109: am->mu = am->muold;
110: }
111: if (am->mu < am->mumin) {
112: am->mu = am->mumin;
113: }
114: return(0);
115: }
117: static PetscErrorCode TaoADMMSetRegularizerType_ADMM(Tao tao, TaoADMMRegularizerType type)
118: {
119: TAO_ADMM *am = (TAO_ADMM*)tao->data;
122: am->regswitch = type;
123: return(0);
124: }
126: static PetscErrorCode TaoADMMGetRegularizerType_ADMM(Tao tao, TaoADMMRegularizerType *type)
127: {
128: TAO_ADMM *am = (TAO_ADMM*)tao->data;
131: *type = am->regswitch;
132: return(0);
133: }
135: static PetscErrorCode TaoADMMSetUpdateType_ADMM(Tao tao, TaoADMMUpdateType type)
136: {
137: TAO_ADMM *am = (TAO_ADMM*)tao->data;
140: am->update = type;
141: return(0);
142: }
144: static PetscErrorCode TaoADMMGetUpdateType_ADMM(Tao tao, TaoADMMUpdateType *type)
145: {
146: TAO_ADMM *am = (TAO_ADMM*)tao->data;
149: *type = am->update;
150: return(0);
151: }
153: /* This routine updates Jacobians with new x,z vectors,
154: * and then updates Ax and Bz vectors, then computes updated residual vector*/
155: static PetscErrorCode ADMMUpdateConstraintResidualVector(Tao tao, Vec x, Vec z, Vec Ax, Vec Bz, Vec residual)
156: {
157: TAO_ADMM *am = (TAO_ADMM*)tao->data;
158: Tao mis,reg;
162: mis = am->subsolverX;
163: reg = am->subsolverZ;
164: TaoComputeJacobianEquality(mis, x, mis->jacobian_equality, mis->jacobian_equality_pre);
165: MatMult(mis->jacobian_equality,x,Ax);
166: TaoComputeJacobianEquality(reg, z, reg->jacobian_equality, reg->jacobian_equality_pre);
167: MatMult(reg->jacobian_equality,z,Bz);
169: VecWAXPY(residual,1.,Bz,Ax);
170: if (am->constraint != NULL) {
171: VecAXPY(residual,-1.,am->constraint);
172: }
173: return(0);
174: }
176: /* Updates Augmented Lagrangians to given routines *
177: * For subsolverX, routine needs to be ComputeObjectiveAndGraidnet
178: * Separate Objective and Gradient routines are not supported. */
179: static PetscErrorCode SubObjGradUpdate(Tao tao, Vec x, PetscReal *f, Vec g, void *ptr)
180: {
181: Tao parent = (Tao)ptr;
182: TAO_ADMM *am = (TAO_ADMM*)parent->data;
184: PetscReal temp,temp2;
185: Vec tempJR;
188: tempJR = am->workJacobianRight;
189: ADMMUpdateConstraintResidualVector(parent, x, am->subsolverZ->solution, am->Ax, am->Bz, am->residual);
190: (*am->ops->misfitobjgrad)(am->subsolverX,x,f,g,am->misfitobjgradP);
192: am->last_misfit_val = *f;
193: /* Objective Add + yT(Ax+Bz-c) + mu/2*||Ax+Bz-c||_2^2 */
194: VecTDot(am->residual,am->y,&temp);
195: VecTDot(am->residual,am->residual,&temp2);
196: *f += temp + (am->mu/2)*temp2;
198: /* Gradient. Add + mu*AT(Ax+Bz-c) + yTA*/
199: MatMultTranspose(tao->jacobian_equality,am->residual,tempJR);
200: VecAXPY(g,am->mu,tempJR);
201: MatMultTranspose(tao->jacobian_equality,am->y,tempJR);
202: VecAXPY(g,1.,tempJR);
203: return(0);
204: }
206: /* Updates Augmented Lagrangians to given routines
207: * For subsolverZ, routine needs to be ComputeObjectiveAndGraidnet
208: * Separate Objective and Gradient routines are not supported. */
209: static PetscErrorCode RegObjGradUpdate(Tao tao, Vec z, PetscReal *f, Vec g, void *ptr)
210: {
211: Tao parent = (Tao)ptr;
212: TAO_ADMM *am = (TAO_ADMM*)parent->data;
214: PetscReal temp,temp2;
215: Vec tempJR;
218: tempJR = am->workJacobianRight;
219: ADMMUpdateConstraintResidualVector(parent, am->subsolverX->solution, z, am->Ax, am->Bz, am->residual);
220: (*am->ops->regobjgrad)(am->subsolverZ,z,f,g,am->regobjgradP);
221: am->last_reg_val= *f;
222: /* Objective Add + yT(Ax+Bz-c) + mu/2*||Ax+Bz-c||_2^2 */
223: VecTDot(am->residual,am->y,&temp);
224: VecTDot(am->residual,am->residual,&temp2);
225: *f += temp + (am->mu/2)*temp2;
227: /* Gradient. Add + mu*BT(Ax+Bz-c) + yTB*/
228: MatMultTranspose(am->subsolverZ->jacobian_equality,am->residual,tempJR);
229: VecAXPY(g,am->mu,tempJR);
230: MatMultTranspose(am->subsolverZ->jacobian_equality,am->y,tempJR);
231: VecAXPY(g,1.,tempJR);
232: return(0);
233: }
235: /* Computes epsilon padded L1 norm lambda*sum(sqrt(x^2+eps^2)-eps */
236: static PetscErrorCode ADMML1EpsilonNorm(Tao tao, Vec x, PetscReal eps, PetscReal *norm)
237: {
238: TAO_ADMM *am = (TAO_ADMM*)tao->data;
239: PetscInt N;
243: VecGetSize(am->workLeft,&N);
244: VecPointwiseMult(am->workLeft,x,x);
245: VecShift(am->workLeft,am->l1epsilon*am->l1epsilon);
246: VecSqrtAbs(am->workLeft);
247: VecSum(am->workLeft,norm);
248: *norm += N*am->l1epsilon;
249: *norm *= am->lambda;
250: return(0);
251: }
253: static PetscErrorCode ADMMInternalHessianUpdate(Mat H, Mat Constraint, PetscBool Identity, void *ptr)
254: {
255: TAO_ADMM *am = (TAO_ADMM*)ptr;
259: switch (am->update) {
260: case (TAO_ADMM_UPDATE_BASIC):
261: break;
262: case (TAO_ADMM_UPDATE_ADAPTIVE):
263: case (TAO_ADMM_UPDATE_ADAPTIVE_RELAXED):
264: if (H && (am->muold != am->mu)) {
265: if (!Identity) {
266: MatAXPY(H,am->mu-am->muold,Constraint,DIFFERENT_NONZERO_PATTERN);
267: } else {
268: MatShift(H,am->mu-am->muold);
269: }
270: }
271: break;
272: }
273: return(0);
274: }
276: /* Updates Hessian - adds second derivative of augmented Lagrangian
277: * H \gets H + \rho*ATA
278: * Here, \rho does not change in TAO_ADMM_UPDATE_BASIC - thus no-op
279: * For ADAPTAIVE,ADAPTIVE_RELAXED,
280: * H \gets H + (\rho-\rhoold)*ATA
281: * Here, we assume that A is linear constraint i.e., doesnt change.
282: * Thus, for both ADAPTIVE, and RELAXED, ATA matrix is pre-set (except for A=I (null case)) see TaoSetUp_ADMM */
283: static PetscErrorCode SubHessianUpdate(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
284: {
285: Tao parent = (Tao)ptr;
286: TAO_ADMM *am = (TAO_ADMM*)parent->data;
290: if (am->Hxchange) {
291: /* Case where Hessian gets updated with respect to x vector input. */
292: (*am->ops->misfithess)(am->subsolverX,x,H,Hpre,am->misfithessP);
293: ADMMInternalHessianUpdate(am->subsolverX->hessian,am->ATA,am->xJI,am);
294: } else if (am->Hxbool) {
295: /* Hessian doesn't get updated. H(x) = c */
296: /* Update Lagrangian only only per TAO call */
297: ADMMInternalHessianUpdate(am->subsolverX->hessian,am->ATA,am->xJI,am);
298: am->Hxbool = PETSC_FALSE;
299: }
300: return(0);
301: }
303: /* Same as SubHessianUpdate, except for B matrix instead of A matrix */
304: static PetscErrorCode RegHessianUpdate(Tao tao, Vec z, Mat H, Mat Hpre, void *ptr)
305: {
306: Tao parent = (Tao)ptr;
307: TAO_ADMM *am = (TAO_ADMM*)parent->data;
312: if (am->Hzchange) {
313: /* Case where Hessian gets updated with respect to x vector input. */
314: (*am->ops->reghess)(am->subsolverZ,z,H,Hpre,am->reghessP);
315: ADMMInternalHessianUpdate(am->subsolverZ->hessian,am->BTB,am->zJI,am);
316: } else if (am->Hzbool) {
317: /* Hessian doesn't get updated. H(x) = c */
318: /* Update Lagrangian only only per TAO call */
319: ADMMInternalHessianUpdate(am->subsolverZ->hessian,am->BTB,am->zJI,am);
320: am->Hzbool = PETSC_FALSE;
321: }
322: return(0);
323: }
325: /* Shell Matrix routine for A matrix.
326: * This gets used when user puts NULL for
327: * TaoSetJacobianEqualityRoutine(tao, NULL,NULL, ...)
328: * Essentially sets A=I*/
329: static PetscErrorCode JacobianIdentity(Mat mat,Vec in,Vec out)
330: {
334: VecCopy(in,out);
335: return(0);
336: }
338: /* Shell Matrix routine for B matrix.
339: * This gets used when user puts NULL for
340: * TaoADMMSetRegularizerConstraintJacobian(tao, NULL,NULL, ...)
341: * Sets B=-I */
342: static PetscErrorCode JacobianIdentityB(Mat mat,Vec in,Vec out)
343: {
347: VecCopy(in,out);
348: VecScale(out,-1.);
349: return(0);
350: }
352: /* Solve f(x) + g(z) s.t. Ax + Bz = c */
353: static PetscErrorCode TaoSolve_ADMM(Tao tao)
354: {
355: TAO_ADMM *am = (TAO_ADMM*)tao->data;
357: PetscInt N;
358: PetscReal reg_func;
359: PetscBool is_reg_shell;
360: Vec tempL;
363: if (am->regswitch != TAO_ADMM_REGULARIZER_SOFT_THRESH) {
364: if (!am->subsolverX->ops->computejacobianequality) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONGSTATE,"Must call TaoADMMSetMisfitConstraintJacobian() first");
365: if (!am->subsolverZ->ops->computejacobianequality) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONGSTATE,"Must call TaoADMMSetRegularizerConstraintJacobian() first");
366: if (am->constraint != NULL) {
367: VecNorm(am->constraint,NORM_2,&am->const_norm);
368: }
369: }
370: tempL = am->workLeft;
371: VecGetSize(tempL,&N);
373: if (am->Hx && am->ops->misfithess) {
374: TaoSetHessianRoutine(am->subsolverX, am->Hx, am->Hx, SubHessianUpdate, tao);
375: }
377: if (!am->zJI) {
378: /* Currently, B is assumed to be a linear system, i.e., not getting updated*/
379: MatTransposeMatMult(am->JB,am->JB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&(am->BTB));
380: }
381: if (!am->xJI) {
382: /* Currently, A is assumed to be a linear system, i.e., not getting updated*/
383: MatTransposeMatMult(am->subsolverX->jacobian_equality,am->subsolverX->jacobian_equality,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&(am->ATA));
384: }
386: is_reg_shell = PETSC_FALSE;
388: PetscObjectTypeCompare((PetscObject)am->subsolverZ, TAOSHELL, &is_reg_shell);
390: if (!is_reg_shell) {
391: switch (am->regswitch) {
392: case (TAO_ADMM_REGULARIZER_USER):
393: if (!am->ops->regobjgrad) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONGSTATE,"Must call TaoADMMSetRegularizerObjectiveAndGradientRoutine() first if one wishes to use TAO_ADMM_REGULARIZER_USER with non-TAOSHELL type");
394: break;
395: case (TAO_ADMM_REGULARIZER_SOFT_THRESH):
396: /* Soft Threshold. */
397: break;
398: }
399: if (am->ops->regobjgrad) {
400: TaoSetObjectiveAndGradientRoutine(am->subsolverZ, RegObjGradUpdate, tao);
401: }
402: if (am->Hz && am->ops->reghess) {
403: TaoSetHessianRoutine(am->subsolverZ, am->Hz, am->Hzpre, RegHessianUpdate, tao);
404: }
405: }
407: switch (am->update) {
408: case TAO_ADMM_UPDATE_BASIC:
409: if (am->subsolverX->hessian) {
410: /* In basic case, Hessian does not get updated w.r.t. to spectral penalty
411: * Here, when A is set, i.e., am->xJI, add mu*ATA to Hessian*/
412: if (!am->xJI) {
413: MatAXPY(am->subsolverX->hessian,am->mu,am->ATA,DIFFERENT_NONZERO_PATTERN);
414: } else {
415: MatShift(am->subsolverX->hessian,am->mu);
416: }
417: }
418: if (am->subsolverZ->hessian && am->regswitch == TAO_ADMM_REGULARIZER_USER) {
419: if (am->regswitch == TAO_ADMM_REGULARIZER_USER && !am->zJI) {
420: MatAXPY(am->subsolverZ->hessian,am->mu,am->BTB,DIFFERENT_NONZERO_PATTERN);
421: } else {
422: MatShift(am->subsolverZ->hessian,am->mu);
423: }
424: }
425: break;
426: case TAO_ADMM_UPDATE_ADAPTIVE:
427: case TAO_ADMM_UPDATE_ADAPTIVE_RELAXED:
428: break;
429: }
431: PetscCitationsRegister(citation,&cited);
432: tao->reason = TAO_CONTINUE_ITERATING;
434: while (tao->reason == TAO_CONTINUE_ITERATING) {
435: if (tao->ops->update) {
436: (*tao->ops->update)(tao, tao->niter, tao->user_update);
437: }
438: VecCopy(am->Bz, am->Bzold);
440: /* x update */
441: TaoSolve(am->subsolverX);
442: TaoComputeJacobianEquality(am->subsolverX, am->subsolverX->solution, am->subsolverX->jacobian_equality, am->subsolverX->jacobian_equality_pre);
443: MatMult(am->subsolverX->jacobian_equality, am->subsolverX->solution,am->Ax);
445: am->Hxbool = PETSC_TRUE;
447: /* z update */
448: switch (am->regswitch) {
449: case TAO_ADMM_REGULARIZER_USER:
450: TaoSolve(am->subsolverZ);
451: break;
452: case TAO_ADMM_REGULARIZER_SOFT_THRESH:
453: /* L1 assumes A,B jacobians are identity nxn matrix */
454: VecWAXPY(am->workJacobianRight,1/am->mu,am->y,am->Ax);
455: TaoSoftThreshold(am->workJacobianRight,-am->lambda/am->mu,am->lambda/am->mu,am->subsolverZ->solution);
456: break;
457: }
458: am->Hzbool = PETSC_TRUE;
459: /* Returns Ax + Bz - c with updated Ax,Bz vectors */
460: ADMMUpdateConstraintResidualVector(tao, am->subsolverX->solution, am->subsolverZ->solution, am->Ax, am->Bz, am->residual);
461: /* Dual variable, y += y + mu*(Ax+Bz-c) */
462: VecWAXPY(am->y, am->mu, am->residual, am->yold);
464: /* stopping tolerance update */
465: TaoADMMToleranceUpdate(tao);
467: /* Updating Spectral Penalty */
468: switch (am->update) {
469: case TAO_ADMM_UPDATE_BASIC:
470: am->muold = am->mu;
471: break;
472: case TAO_ADMM_UPDATE_ADAPTIVE:
473: case TAO_ADMM_UPDATE_ADAPTIVE_RELAXED:
474: if (tao->niter == 0) {
475: VecCopy(am->y, am->y0);
476: VecWAXPY(am->residual, 1., am->Ax, am->Bzold);
477: if (am->constraint) {
478: VecAXPY(am->residual, -1., am->constraint);
479: }
480: VecWAXPY(am->yhatold, -am->mu, am->residual, am->yold);
481: VecCopy(am->Ax, am->Axold);
482: VecCopy(am->Bz, am->Bz0);
483: am->muold = am->mu;
484: } else if (tao->niter % am->T == 1) {
485: /* we have compute Bzold in a previous iteration, and we computed Ax above */
486: VecWAXPY(am->residual, 1., am->Ax, am->Bzold);
487: if (am->constraint) {
488: VecAXPY(am->residual, -1., am->constraint);
489: }
490: VecWAXPY(am->yhat, -am->mu, am->residual, am->yold);
491: AdaptiveADMMPenaltyUpdate(tao);
492: VecCopy(am->Ax, am->Axold);
493: VecCopy(am->Bz, am->Bz0);
494: VecCopy(am->yhat, am->yhatold);
495: VecCopy(am->y, am->y0);
496: } else {
497: am->muold = am->mu;
498: }
499: break;
500: default:
501: break;
502: }
503: tao->niter++;
505: /* Calculate original function values. misfit part was done in TaoADMMToleranceUpdate*/
506: switch (am->regswitch) {
507: case TAO_ADMM_REGULARIZER_USER:
508: if (is_reg_shell) {
509: ADMML1EpsilonNorm(tao,am->subsolverZ->solution,am->l1epsilon,®_func);
510: } else {
511: (*am->ops->regobjgrad)(am->subsolverZ,am->subsolverX->solution,®_func,tempL,am->regobjgradP);
512: }
513: break;
514: case TAO_ADMM_REGULARIZER_SOFT_THRESH:
515: ADMML1EpsilonNorm(tao,am->subsolverZ->solution,am->l1epsilon,®_func);
516: break;
517: }
518: VecCopy(am->y,am->yold);
519: ADMMUpdateConstraintResidualVector(tao, am->subsolverX->solution, am->subsolverZ->solution, am->Ax, am->Bz, am->residual);
520: VecNorm(am->residual,NORM_2,&am->resnorm);
521: TaoLogConvergenceHistory(tao,am->last_misfit_val + reg_func,am->dualres,am->resnorm,tao->ksp_its);
523: TaoMonitor(tao,tao->niter,am->last_misfit_val + reg_func,am->dualres,am->resnorm,1.0);
524: (*tao->ops->convergencetest)(tao,tao->cnvP);
525: }
526: /* Update vectors */
527: VecCopy(am->subsolverX->solution,tao->solution);
528: VecCopy(am->subsolverX->gradient,tao->gradient);
529: PetscObjectCompose((PetscObject)am->subsolverX,"TaoGetADMMParentTao_ADMM", NULL);
530: PetscObjectCompose((PetscObject)am->subsolverZ,"TaoGetADMMParentTao_ADMM", NULL);
531: PetscObjectComposeFunction((PetscObject)tao,"TaoADMMSetRegularizerType_C",NULL);
532: PetscObjectComposeFunction((PetscObject)tao,"TaoADMMGetRegularizerType_C",NULL);
533: PetscObjectComposeFunction((PetscObject)tao,"TaoADMMSetUpdateType_C",NULL);
534: PetscObjectComposeFunction((PetscObject)tao,"TaoADMMGetUpdateType_C",NULL);
535: return(0);
536: }
538: static PetscErrorCode TaoSetFromOptions_ADMM(PetscOptionItems *PetscOptionsObject,Tao tao)
539: {
540: TAO_ADMM *am = (TAO_ADMM*)tao->data;
544: PetscOptionsHead(PetscOptionsObject,"ADMM problem that solves f(x) in a form of f(x) + g(z) subject to x - z = 0. Norm 1 and 2 are supported. Different subsolver routines can be selected. ");
545: PetscOptionsReal("-tao_admm_regularizer_coefficient","regularizer constant","",am->lambda,&am->lambda,NULL);
546: PetscOptionsReal("-tao_admm_spectral_penalty","Constant for Augmented Lagrangian term.","",am->mu,&am->mu,NULL);
547: PetscOptionsReal("-tao_admm_relaxation_parameter","x relaxation parameter for Z update.","",am->gamma,&am->gamma,NULL);
548: PetscOptionsReal("-tao_admm_tolerance_update_factor","ADMM dynamic tolerance update factor.","",am->tol,&am->tol,NULL);
549: PetscOptionsReal("-tao_admm_spectral_penalty_update_factor","ADMM spectral penalty update curvature safeguard value.","",am->orthval,&am->orthval,NULL);
550: PetscOptionsReal("-tao_admm_minimum_spectral_penalty","Set ADMM minimum spectral penalty.","",am->mumin,&am->mumin,NULL);
551: PetscOptionsEnum("-tao_admm_dual_update","Lagrangian dual update policy","TaoADMMUpdateType",
552: TaoADMMUpdateTypes,(PetscEnum)am->update,(PetscEnum*)&am->update,NULL);
553: PetscOptionsEnum("-tao_admm_regularizer_type","ADMM regularizer update rule","TaoADMMRegularizerType",
554: TaoADMMRegularizerTypes,(PetscEnum)am->regswitch,(PetscEnum*)&am->regswitch,NULL);
555: PetscOptionsTail();
556: TaoSetFromOptions(am->subsolverX);
557: if (am->regswitch != TAO_ADMM_REGULARIZER_SOFT_THRESH) {
558: TaoSetFromOptions(am->subsolverZ);
559: }
560: return(0);
561: }
563: static PetscErrorCode TaoView_ADMM(Tao tao,PetscViewer viewer)
564: {
565: TAO_ADMM *am = (TAO_ADMM*)tao->data;
569: PetscViewerASCIIPushTab(viewer);
570: TaoView(am->subsolverX,viewer);
571: TaoView(am->subsolverZ,viewer);
572: PetscViewerASCIIPopTab(viewer);
573: return(0);
574: }
576: static PetscErrorCode TaoSetUp_ADMM(Tao tao)
577: {
578: TAO_ADMM *am = (TAO_ADMM*)tao->data;
580: PetscInt n,N,M;
583: VecGetLocalSize(tao->solution,&n);
584: VecGetSize(tao->solution,&N);
585: /* If Jacobian is given as NULL, it means Jacobian is identity matrix with size of solution vector */
586: if (!am->JB) {
587: am->zJI = PETSC_TRUE;
588: MatCreateShell(PetscObjectComm((PetscObject)tao),n,n,PETSC_DETERMINE,PETSC_DETERMINE,NULL,&am->JB);
589: MatShellSetOperation(am->JB,MATOP_MULT,(void (*)(void))JacobianIdentityB);
590: MatShellSetOperation(am->JB,MATOP_MULT_TRANSPOSE,(void (*)(void))JacobianIdentityB);
591: MatShellSetOperation(am->JB,MATOP_TRANSPOSE_MAT_MULT,(void (*)(void))JacobianIdentityB);
592: am->JBpre = am->JB;
593: }
594: if (!am->JA) {
595: am->xJI = PETSC_TRUE;
596: MatCreateShell(PetscObjectComm((PetscObject)tao),n,n,PETSC_DETERMINE,PETSC_DETERMINE,NULL,&am->JA);
597: MatShellSetOperation(am->JA,MATOP_MULT,(void (*)(void))JacobianIdentity);
598: MatShellSetOperation(am->JA,MATOP_MULT_TRANSPOSE,(void (*)(void))JacobianIdentity);
599: MatShellSetOperation(am->JA,MATOP_TRANSPOSE_MAT_MULT,(void (*)(void))JacobianIdentity);
600: am->JApre = am->JA;
601: }
602: MatCreateVecs(am->JA,NULL,&am->Ax);
603: if (!tao->gradient) {
604: VecDuplicate(tao->solution,&tao->gradient);
605: }
606: TaoSetInitialVector(am->subsolverX, tao->solution);
607: if (!am->z) {
608: VecDuplicate(tao->solution,&am->z);
609: VecSet(am->z,0.0);
610: }
611: TaoSetInitialVector(am->subsolverZ, am->z);
612: if (!am->workLeft) {
613: VecDuplicate(tao->solution,&am->workLeft);
614: }
615: if (!am->Axold) {
616: VecDuplicate(am->Ax,&am->Axold);
617: }
618: if (!am->workJacobianRight) {
619: VecDuplicate(am->Ax,&am->workJacobianRight);
620: }
621: if (!am->workJacobianRight2) {
622: VecDuplicate(am->Ax,&am->workJacobianRight2);
623: }
624: if (!am->Bz) {
625: VecDuplicate(am->Ax,&am->Bz);
626: }
627: if (!am->Bzold) {
628: VecDuplicate(am->Ax,&am->Bzold);
629: }
630: if (!am->Bz0) {
631: VecDuplicate(am->Ax,&am->Bz0);
632: }
633: if (!am->y) {
634: VecDuplicate(am->Ax,&am->y);
635: VecSet(am->y,0.0);
636: }
637: if (!am->yold) {
638: VecDuplicate(am->Ax,&am->yold);
639: VecSet(am->yold,0.0);
640: }
641: if (!am->y0) {
642: VecDuplicate(am->Ax,&am->y0);
643: VecSet(am->y0,0.0);
644: }
645: if (!am->yhat) {
646: VecDuplicate(am->Ax,&am->yhat);
647: VecSet(am->yhat,0.0);
648: }
649: if (!am->yhatold) {
650: VecDuplicate(am->Ax,&am->yhatold);
651: VecSet(am->yhatold,0.0);
652: }
653: if (!am->residual) {
654: VecDuplicate(am->Ax,&am->residual);
655: VecSet(am->residual,0.0);
656: }
657: if (!am->constraint) {
658: am->constraint = NULL;
659: } else {
660: VecGetSize(am->constraint,&M);
661: if (M != N) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONGSTATE,"Solution vector and constraint vector must be of same size!");
662: }
664: /* Save changed tao tolerance for adaptive tolerance */
665: if (tao->gatol_changed) {
666: am->gatol_admm = tao->gatol;
667: }
668: if (tao->catol_changed) {
669: am->catol_admm = tao->catol;
670: }
672: /*Update spectral and dual elements to X subsolver */
673: TaoSetObjectiveAndGradientRoutine(am->subsolverX, SubObjGradUpdate, tao);
674: TaoSetJacobianEqualityRoutine(am->subsolverX,am->JA,am->JApre, am->ops->misfitjac, am->misfitjacobianP);
675: TaoSetJacobianEqualityRoutine(am->subsolverZ,am->JB,am->JBpre, am->ops->regjac, am->regjacobianP);
676: return(0);
677: }
679: static PetscErrorCode TaoDestroy_ADMM(Tao tao)
680: {
681: TAO_ADMM *am = (TAO_ADMM*)tao->data;
685: VecDestroy(&am->z);
686: VecDestroy(&am->Ax);
687: VecDestroy(&am->Axold);
688: VecDestroy(&am->Bz);
689: VecDestroy(&am->Bzold);
690: VecDestroy(&am->Bz0);
691: VecDestroy(&am->residual);
692: VecDestroy(&am->y);
693: VecDestroy(&am->yold);
694: VecDestroy(&am->y0);
695: VecDestroy(&am->yhat);
696: VecDestroy(&am->yhatold);
697: VecDestroy(&am->workLeft);
698: VecDestroy(&am->workJacobianRight);
699: VecDestroy(&am->workJacobianRight2);
701: MatDestroy(&am->JA);
702: MatDestroy(&am->JB);
703: if (!am->xJI) {
704: MatDestroy(&am->JApre);
705: }
706: if (!am->zJI) {
707: MatDestroy(&am->JBpre);
708: }
709: if (am->Hx) {
710: MatDestroy(&am->Hx);
711: MatDestroy(&am->Hxpre);
712: }
713: if (am->Hz) {
714: MatDestroy(&am->Hz);
715: MatDestroy(&am->Hzpre);
716: }
717: MatDestroy(&am->ATA);
718: MatDestroy(&am->BTB);
719: TaoDestroy(&am->subsolverX);
720: TaoDestroy(&am->subsolverZ);
721: am->parent = NULL;
722: PetscFree(tao->data);
723: return(0);
724: }
726: /*MC
728: TAOADMM - Alternating direction method of multipliers method fo solving linear problems with
729: constraints. in a min_x f(x) + g(z) s.t. Ax+Bz=c.
730: This algorithm employs two sub Tao solvers, of which type can be specificed
731: by the user. User need to provide ObjectiveAndGradient routine, and/or HessianRoutine for both subsolvers.
732: Hessians can be given boolean flag determining whether they change with respect to a input vector. This can be set via
733: TaoADMMSet{Misfit,Regularizer}HessianChangeStatus.
734: Second subsolver does support TAOSHELL. It should be noted that L1-norm is used for objective value for TAOSHELL type.
735: There is option to set regularizer option, and currently soft-threshold is implemented. For spectral penalty update,
736: currently there are baisc option and adaptive option.
737: Constraint is set at Ax+Bz=c, and A and B can be set with TaoADMMSet{Misfit,Regularizer}ConstraintJacobian.
738: c can be set with TaoADMMSetConstraintVectorRHS.
739: The user can also provide regularizer weight for second subsolver.
741: References:
742: . 1. - Xu, Zheng and Figueiredo, Mario A. T. and Yuan, Xiaoming and Studer, Christoph and Goldstein, Tom
743: "Adaptive Relaxed ADMM: Convergence Theory and Practical Implementation"
744: The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), July, 2017.
746: Options Database Keys:
747: + -tao_admm_regularizer_coefficient - regularizer constant (default 1.e-6)
748: . -tao_admm_spectral_penalty - Constant for Augmented Lagrangian term (default 1.)
749: . -tao_admm_relaxation_parameter - relaxation parameter for Z update (default 1.)
750: . -tao_admm_tolerance_update_factor - ADMM dynamic tolerance update factor (default 1.e-12)
751: . -tao_admm_spectral_penalty_update_factor - ADMM spectral penalty update curvature safeguard value (default 0.2)
752: . -tao_admm_minimum_spectral_penalty - Set ADMM minimum spectral penalty (default 0)
753: . -tao_admm_dual_update - Lagrangian dual update policy ("basic","adaptive","adaptive-relaxed") (default "basic")
754: - -tao_admm_regularizer_type - ADMM regularizer update rule ("user","soft-threshold") (default "soft-threshold")
756: Level: beginner
758: .seealso: TaoADMMSetMisfitHessianChangeStatus(), TaoADMMSetRegHessianChangeStatus(), TaoADMMGetSpectralPenalty(),
759: TaoADMMGetMisfitSubsolver(), TaoADMMGetRegularizationSubsolver(), TaoADMMSetConstraintVectorRHS(),
760: TaoADMMSetMinimumSpectralPenalty(), TaoADMMSetRegularizerCoefficient(),
761: TaoADMMSetRegularizerConstraintJacobian(), TaoADMMSetMisfitConstraintJacobian(),
762: TaoADMMSetMisfitObjectiveAndGradientRoutine(), TaoADMMSetMisfitHessianRoutine(),
763: TaoADMMSetRegularizerObjectiveAndGradientRoutine(), TaoADMMSetRegularizerHessianRoutine(),
764: TaoGetADMMParentTao(), TaoADMMGetDualVector(), TaoADMMSetRegularizerType(),
765: TaoADMMGetRegularizerType(), TaoADMMSetUpdateType(), TaoADMMGetUpdateType()
766: M*/
768: PETSC_EXTERN PetscErrorCode TaoCreate_ADMM(Tao tao)
769: {
770: TAO_ADMM *am;
774: PetscNewLog(tao,&am);
776: tao->ops->destroy = TaoDestroy_ADMM;
777: tao->ops->setup = TaoSetUp_ADMM;
778: tao->ops->setfromoptions = TaoSetFromOptions_ADMM;
779: tao->ops->view = TaoView_ADMM;
780: tao->ops->solve = TaoSolve_ADMM;
782: tao->data = (void*)am;
783: am->l1epsilon = 1e-6;
784: am->lambda = 1e-4;
785: am->mu = 1.;
786: am->muold = 0.;
787: am->mueps = PETSC_MACHINE_EPSILON;
788: am->mumin = 0.;
789: am->orthval = 0.2;
790: am->T = 2;
791: am->parent = tao;
792: am->update = TAO_ADMM_UPDATE_BASIC;
793: am->regswitch = TAO_ADMM_REGULARIZER_SOFT_THRESH;
794: am->tol = PETSC_SMALL;
795: am->const_norm = 0;
796: am->resnorm = 0;
797: am->dualres = 0;
798: am->ops->regobjgrad = NULL;
799: am->ops->reghess = NULL;
800: am->gamma = 1;
801: am->regobjgradP = NULL;
802: am->reghessP = NULL;
803: am->gatol_admm = 1e-8;
804: am->catol_admm = 0;
805: am->Hxchange = PETSC_TRUE;
806: am->Hzchange = PETSC_TRUE;
807: am->Hzbool = PETSC_TRUE;
808: am->Hxbool = PETSC_TRUE;
810: TaoCreate(PetscObjectComm((PetscObject)tao),&am->subsolverX);
811: TaoSetOptionsPrefix(am->subsolverX,"misfit_");
812: PetscObjectIncrementTabLevel((PetscObject)am->subsolverX,(PetscObject)tao,1);
813: TaoCreate(PetscObjectComm((PetscObject)tao),&am->subsolverZ);
814: TaoSetOptionsPrefix(am->subsolverZ,"reg_");
815: PetscObjectIncrementTabLevel((PetscObject)am->subsolverZ,(PetscObject)tao,1);
817: TaoSetType(am->subsolverX,TAONLS);
818: TaoSetType(am->subsolverZ,TAONLS);
819: PetscObjectCompose((PetscObject)am->subsolverX,"TaoGetADMMParentTao_ADMM", (PetscObject) tao);
820: PetscObjectCompose((PetscObject)am->subsolverZ,"TaoGetADMMParentTao_ADMM", (PetscObject) tao);
821: PetscObjectComposeFunction((PetscObject)tao,"TaoADMMSetRegularizerType_C",TaoADMMSetRegularizerType_ADMM);
822: PetscObjectComposeFunction((PetscObject)tao,"TaoADMMGetRegularizerType_C",TaoADMMGetRegularizerType_ADMM);
823: PetscObjectComposeFunction((PetscObject)tao,"TaoADMMSetUpdateType_C",TaoADMMSetUpdateType_ADMM);
824: PetscObjectComposeFunction((PetscObject)tao,"TaoADMMGetUpdateType_C",TaoADMMGetUpdateType_ADMM);
825: return(0);
826: }
828: /*@
829: TaoADMMSetMisfitHessianChangeStatus - Set boolean that determines whether Hessian matrix of misfit subsolver changes with respect to input vector.
831: Collective on Tao
833: Input Parameters:
834: + tao - the Tao solver context.
835: - b - the Hessian matrix change status boolean, PETSC_FALSE when the Hessian matrix does not change, TRUE otherwise.
837: Level: advanced
839: .seealso: TAOADMM
841: @*/
842: PetscErrorCode TaoADMMSetMisfitHessianChangeStatus(Tao tao, PetscBool b)
843: {
844: TAO_ADMM *am = (TAO_ADMM*)tao->data;
847: am->Hxchange = b;
848: return(0);
849: }
851: /*@
852: TaoADMMSetRegHessianChangeStatus - Set boolean that determines whether Hessian matrix of regularization subsolver changes with respect to input vector.
855: Collective on Tao
857: Input Parameters:
858: + tao - the Tao solver context
859: - b - the Hessian matrix change status boolean, PETSC_FALSE when the Hessian matrix does not change, TRUE otherwise.
861: Level: advanced
863: .seealso: TAOADMM
865: @*/
866: PetscErrorCode TaoADMMSetRegHessianChangeStatus(Tao tao, PetscBool b)
867: {
868: TAO_ADMM *am = (TAO_ADMM*)tao->data;
871: am->Hzchange = b;
872: return(0);
873: }
875: /*@
876: TaoADMMSetSpectralPenalty - Set the spectral penalty (mu) value
878: Collective on Tao
880: Input Parameters:
881: + tao - the Tao solver context
882: - mu - spectral penalty
884: Level: advanced
886: .seealso: TaoADMMSetMinimumSpectralPenalty(), TAOADMM
887: @*/
888: PetscErrorCode TaoADMMSetSpectralPenalty(Tao tao, PetscReal mu)
889: {
890: TAO_ADMM *am = (TAO_ADMM*)tao->data;
893: am->mu = mu;
894: return(0);
895: }
897: /*@
898: TaoADMMGetSpectralPenalty - Get the spectral penalty (mu) value
900: Collective on Tao
902: Input Parameter:
903: . tao - the Tao solver context
905: Output Parameter:
906: . mu - spectral penalty
908: Level: advanced
910: .seealso: TaoADMMSetMinimumSpectralPenalty(), TaoADMMSetSpectralPenalty(), TAOADMM
911: @*/
912: PetscErrorCode TaoADMMGetSpectralPenalty(Tao tao, PetscReal *mu)
913: {
914: TAO_ADMM *am = (TAO_ADMM*)tao->data;
919: *mu = am->mu;
920: return(0);
921: }
923: /*@
924: TaoADMMGetMisfitSubsolver - Get the pointer to the misfit subsolver inside ADMM
926: Collective on Tao
928: Input Parameter:
929: . tao - the Tao solver context
931: Output Parameter:
932: . misfit - the Tao subsolver context
934: Level: advanced
936: .seealso: TAOADMM
938: @*/
939: PetscErrorCode TaoADMMGetMisfitSubsolver(Tao tao, Tao *misfit)
940: {
941: TAO_ADMM *am = (TAO_ADMM*)tao->data;
944: *misfit = am->subsolverX;
945: return(0);
946: }
948: /*@
949: TaoADMMGetRegularizationSubsolver - Get the pointer to the regularization subsolver inside ADMM
951: Collective on Tao
953: Input Parameter:
954: . tao - the Tao solver context
956: Output Parameter:
957: . reg - the Tao subsolver context
959: Level: advanced
961: .seealso: TAOADMM
963: @*/
964: PetscErrorCode TaoADMMGetRegularizationSubsolver(Tao tao, Tao *reg)
965: {
966: TAO_ADMM *am = (TAO_ADMM*)tao->data;
969: *reg = am->subsolverZ;
970: return(0);
971: }
973: /*@
974: TaoADMMSetConstraintVectorRHS - Set the RHS constraint vector for ADMM
976: Collective on Tao
978: Input Parameters:
979: + tao - the Tao solver context
980: - c - RHS vector
982: Level: advanced
984: .seealso: TAOADMM
986: @*/
987: PetscErrorCode TaoADMMSetConstraintVectorRHS(Tao tao, Vec c)
988: {
989: TAO_ADMM *am = (TAO_ADMM*)tao->data;
992: am->constraint = c;
993: return(0);
994: }
996: /*@
997: TaoADMMSetMinimumSpectralPenalty - Set the minimum value for the spectral penalty
999: Collective on Tao
1001: Input Parameters:
1002: + tao - the Tao solver context
1003: - mu - minimum spectral penalty value
1005: Level: advanced
1007: .seealso: TaoADMMGetSpectralPenalty(), TAOADMM
1008: @*/
1009: PetscErrorCode TaoADMMSetMinimumSpectralPenalty(Tao tao, PetscReal mu)
1010: {
1011: TAO_ADMM *am = (TAO_ADMM*)tao->data;
1014: am->mumin= mu;
1015: return(0);
1016: }
1018: /*@
1019: TaoADMMSetRegularizerCoefficient - Set the regularization coefficient lambda for L1 norm regularization case
1021: Collective on Tao
1023: Input Parameters:
1024: + tao - the Tao solver context
1025: - lambda - L1-norm regularizer coefficient
1027: Level: advanced
1029: .seealso: TaoADMMSetMisfitConstraintJacobian(), TaoADMMSetRegularizerConstraintJacobian(), TAOADMM
1031: @*/
1032: PetscErrorCode TaoADMMSetRegularizerCoefficient(Tao tao, PetscReal lambda)
1033: {
1034: TAO_ADMM *am = (TAO_ADMM*)tao->data;
1037: am->lambda = lambda;
1038: return(0);
1039: }
1041: /*@C
1042: TaoADMMSetMisfitConstraintJacobian - Set the constraint matrix B for the ADMM algorithm. Matrix B constrains the z variable.
1044: Collective on Tao
1046: Input Parameters:
1047: + tao - the Tao solver context
1048: . J - user-created regularizer constraint Jacobian matrix
1049: . Jpre - user-created regularizer Jacobian constraint preconditioner matrix
1050: . func - function pointer for the regularizer constraint Jacobian update function
1051: - ctx - user context for the regularizer Hessian
1053: Level: advanced
1055: .seealso: TaoADMMSetRegularizerCoefficient(), TaoADMMSetRegularizerConstraintJacobian(), TAOADMM
1057: @*/
1058: PetscErrorCode TaoADMMSetMisfitConstraintJacobian(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
1059: {
1060: TAO_ADMM *am = (TAO_ADMM*)tao->data;
1065: if (J) {
1068: }
1069: if (Jpre) {
1072: }
1073: if (ctx) am->misfitjacobianP = ctx;
1074: if (func) am->ops->misfitjac = func;
1076: if (J) {
1077: PetscObjectReference((PetscObject)J);
1078: MatDestroy(&am->JA);
1079: am->JA = J;
1080: }
1081: if (Jpre) {
1082: PetscObjectReference((PetscObject)Jpre);
1083: MatDestroy(&am->JApre);
1084: am->JApre = Jpre;
1085: }
1086: return(0);
1087: }
1089: /*@C
1090: TaoADMMSetRegularizerConstraintJacobian - Set the constraint matrix B for ADMM algorithm. Matrix B constraints z variable.
1092: Collective on Tao
1094: Input Parameters:
1095: + tao - the Tao solver context
1096: . J - user-created regularizer constraint Jacobian matrix
1097: . Jpre - user-created regularizer Jacobian constraint preconditioner matrix
1098: . func - function pointer for the regularizer constraint Jacobian update function
1099: - ctx - user context for the regularizer Hessian
1101: Level: advanced
1103: .seealso: TaoADMMSetRegularizerCoefficient(), TaoADMMSetMisfitConstraintJacobian(), TAOADMM
1105: @*/
1106: PetscErrorCode TaoADMMSetRegularizerConstraintJacobian(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
1107: {
1108: TAO_ADMM *am = (TAO_ADMM*)tao->data;
1113: if (J) {
1116: }
1117: if (Jpre) {
1120: }
1121: if (ctx) am->regjacobianP = ctx;
1122: if (func) am->ops->regjac = func;
1124: if (J) {
1125: PetscObjectReference((PetscObject)J);
1126: MatDestroy(&am->JB);
1127: am->JB = J;
1128: }
1129: if (Jpre) {
1130: PetscObjectReference((PetscObject)Jpre);
1131: MatDestroy(&am->JBpre);
1132: am->JBpre = Jpre;
1133: }
1134: return(0);
1135: }
1137: /*@C
1138: TaoADMMSetMisfitObjectiveAndGradientRoutine - Sets the user-defined misfit call-back function
1140: Collective on tao
1142: Input Parameters:
1143: + tao - the Tao context
1144: . func - function pointer for the misfit value and gradient evaluation
1145: - ctx - user context for the misfit
1147: Level: advanced
1149: .seealso: TAOADMM
1151: @*/
1152: PetscErrorCode TaoADMMSetMisfitObjectiveAndGradientRoutine(Tao tao, PetscErrorCode (*func)(Tao, Vec, PetscReal*, Vec, void*), void *ctx)
1153: {
1154: TAO_ADMM *am = (TAO_ADMM*)tao->data;
1158: am->misfitobjgradP = ctx;
1159: am->ops->misfitobjgrad = func;
1160: return(0);
1161: }
1163: /*@C
1164: TaoADMMSetMisfitHessianRoutine - Sets the user-defined misfit Hessian call-back
1165: function into the algorithm, to be used for subsolverX.
1167: Collective on tao
1169: Input Parameters:
1170: + tao - the Tao context
1171: . H - user-created matrix for the Hessian of the misfit term
1172: . Hpre - user-created matrix for the preconditioner of Hessian of the misfit term
1173: . func - function pointer for the misfit Hessian evaluation
1174: - ctx - user context for the misfit Hessian
1176: Level: advanced
1178: .seealso: TAOADMM
1180: @*/
1181: PetscErrorCode TaoADMMSetMisfitHessianRoutine(Tao tao, Mat H, Mat Hpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
1182: {
1183: TAO_ADMM *am = (TAO_ADMM*)tao->data;
1188: if (H) {
1191: }
1192: if (Hpre) {
1195: }
1196: if (ctx) {
1197: am->misfithessP = ctx;
1198: }
1199: if (func) {
1200: am->ops->misfithess = func;
1201: }
1202: if (H) {
1203: PetscObjectReference((PetscObject)H);
1204: MatDestroy(&am->Hx);
1205: am->Hx = H;
1206: }
1207: if (Hpre) {
1208: PetscObjectReference((PetscObject)Hpre);
1209: MatDestroy(&am->Hxpre);
1210: am->Hxpre = Hpre;
1211: }
1212: return(0);
1213: }
1215: /*@C
1216: TaoADMMSetRegularizerObjectiveAndGradientRoutine - Sets the user-defined regularizer call-back function
1218: Collective on tao
1220: Input Parameters:
1221: + tao - the Tao context
1222: . func - function pointer for the regularizer value and gradient evaluation
1223: - ctx - user context for the regularizer
1225: Level: advanced
1227: .seealso: TAOADMM
1229: @*/
1230: PetscErrorCode TaoADMMSetRegularizerObjectiveAndGradientRoutine(Tao tao, PetscErrorCode (*func)(Tao, Vec, PetscReal*, Vec, void*), void *ctx)
1231: {
1232: TAO_ADMM *am = (TAO_ADMM*)tao->data;
1236: am->regobjgradP = ctx;
1237: am->ops->regobjgrad = func;
1238: return(0);
1239: }
1241: /*@C
1242: TaoADMMSetRegularizerHessianRoutine - Sets the user-defined regularizer Hessian call-back
1243: function, to be used for subsolverZ.
1245: Collective on tao
1247: Input Parameters:
1248: + tao - the Tao context
1249: . H - user-created matrix for the Hessian of the regularization term
1250: . Hpre - user-created matrix for the preconditioner of Hessian of the regularization term
1251: . func - function pointer for the regularizer Hessian evaluation
1252: - ctx - user context for the regularizer Hessian
1254: Level: advanced
1256: .seealso: TAOADMM
1258: @*/
1259: PetscErrorCode TaoADMMSetRegularizerHessianRoutine(Tao tao, Mat H, Mat Hpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
1260: {
1261: TAO_ADMM *am = (TAO_ADMM*)tao->data;
1266: if (H) {
1269: }
1270: if (Hpre) {
1273: }
1274: if (ctx) {
1275: am->reghessP = ctx;
1276: }
1277: if (func) {
1278: am->ops->reghess = func;
1279: }
1280: if (H) {
1281: PetscObjectReference((PetscObject)H);
1282: MatDestroy(&am->Hz);
1283: am->Hz = H;
1284: }
1285: if (Hpre) {
1286: PetscObjectReference((PetscObject)Hpre);
1287: MatDestroy(&am->Hzpre);
1288: am->Hzpre = Hpre;
1289: }
1290: return(0);
1291: }
1293: /*@
1294: TaoGetADMMParentTao - Gets pointer to parent ADMM tao, used by inner subsolver.
1296: Collective on tao
1298: Input Parameter:
1299: . tao - the Tao context
1301: Output Parameter:
1302: . admm_tao - the parent Tao context
1304: Level: advanced
1306: .seealso: TAOADMM
1308: @*/
1309: PetscErrorCode TaoGetADMMParentTao(Tao tao, Tao *admm_tao)
1310: {
1315: PetscObjectQuery((PetscObject)tao,"TaoGetADMMParentTao_ADMM", (PetscObject*) admm_tao);
1316: return(0);
1317: }
1319: /*@
1320: TaoADMMGetDualVector - Returns the dual vector associated with the current TAOADMM state
1322: Not Collective
1324: Input Parameter:
1325: . tao - the Tao context
1327: Output Parameter:
1328: . Y - the current solution
1330: Level: intermediate
1332: .seealso: TAOADMM
1334: @*/
1335: PetscErrorCode TaoADMMGetDualVector(Tao tao, Vec *Y)
1336: {
1337: TAO_ADMM *am = (TAO_ADMM*)tao->data;
1341: *Y = am->y;
1342: return(0);
1343: }
1345: /*@
1346: TaoADMMSetRegularizerType - Set regularizer type for ADMM routine
1348: Not Collective
1350: Input Parameter:
1351: + tao - the Tao context
1352: - type - regularizer type
1354: Options Database:
1355: . -tao_admm_regularizer_type <admm_regularizer_user,admm_regularizer_soft_thresh>
1357: Level: intermediate
1359: .seealso: TaoADMMGetRegularizerType(), TaoADMMRegularizerType, TAOADMM
1360: @*/
1361: PetscErrorCode TaoADMMSetRegularizerType(Tao tao, TaoADMMRegularizerType type)
1362: {
1368: PetscTryMethod(tao,"TaoADMMSetRegularizerType_C",(Tao,TaoADMMRegularizerType),(tao,type));
1369: return(0);
1370: }
1372: /*@
1373: TaoADMMGetRegularizerType - Gets the type of regularizer routine for ADMM
1375: Not Collective
1377: Input Parameter:
1378: . tao - the Tao context
1380: Output Parameter:
1381: . type - the type of regularizer
1383: Options Database:
1384: . -tao_admm_regularizer_type <admm_regularizer_user,admm_regularizer_soft_thresh>
1386: Level: intermediate
1388: .seealso: TaoADMMSetRegularizerType(), TaoADMMRegularizerType, TAOADMM
1389: @*/
1390: PetscErrorCode TaoADMMGetRegularizerType(Tao tao, TaoADMMRegularizerType *type)
1391: {
1396: PetscUseMethod(tao,"TaoADMMGetRegularizerType_C",(Tao,TaoADMMRegularizerType*),(tao,type));
1397: return(0);
1398: }
1400: /*@
1401: TaoADMMSetUpdateType - Set update routine for ADMM routine
1403: Not Collective
1405: Input Parameter:
1406: + tao - the Tao context
1407: - type - spectral parameter update type
1409: Level: intermediate
1411: .seealso: TaoADMMGetUpdateType(), TaoADMMUpdateType, TAOADMM
1412: @*/
1413: PetscErrorCode TaoADMMSetUpdateType(Tao tao, TaoADMMUpdateType type)
1414: {
1420: PetscTryMethod(tao,"TaoADMMSetUpdateType_C",(Tao,TaoADMMUpdateType),(tao,type));
1421: return(0);
1422: }
1424: /*@
1425: TaoADMMGetUpdateType - Gets the type of spectral penalty update routine for ADMM
1427: Not Collective
1429: Input Parameter:
1430: . tao - the Tao context
1432: Output Parameter:
1433: . type - the type of spectral penalty update routine
1435: Level: intermediate
1437: .seealso: TaoADMMSetUpdateType(), TaoADMMUpdateType, TAOADMM
1438: @*/
1439: PetscErrorCode TaoADMMGetUpdateType(Tao tao, TaoADMMUpdateType *type)
1440: {
1445: PetscUseMethod(tao,"TaoADMMGetUpdateType_C",(Tao,TaoADMMUpdateType*),(tao,type));
1446: return(0);
1447: }