Actual source code: ex4.c

      1: static char help[] = "Simple example to test separable objective optimizers.\n";

      3: #include <petsc.h>
      4: #include <petsctao.h>
      5: #include <petscvec.h>
      6: #include <petscmath.h>

      8: #define NWORKLEFT 4
      9: #define NWORKRIGHT 12

     11: typedef struct _UserCtx
     12: {
     13:   PetscInt    m;                     /* The row dimension of F */
     14:   PetscInt    n;                     /* The column dimension of F */
     15:   PetscInt    matops;                /* Matrix format. 0 for stencil, 1 for random */
     16:   PetscInt    iter;                  /* Numer of iterations for ADMM */
     17:   PetscReal   hStart;                /* Starting point for Taylor test */
     18:   PetscReal   hFactor;               /* Taylor test step factor */
     19:   PetscReal   hMin;                  /* Taylor test end goal */
     20:   PetscReal   alpha;                 /* regularization constant applied to || x ||_p */
     21:   PetscReal   eps;                   /* small constant for approximating gradient of || x ||_1 */
     22:   PetscReal   mu;                    /* the augmented Lagrangian term in ADMM */
     23:   PetscReal   abstol;
     24:   PetscReal   reltol;
     25:   Mat         F;                     /* matrix in least squares component $(1/2) * || F x - d ||_2^2$ */
     26:   Mat         W;                     /* Workspace matrix. ATA */
     27:   Mat         Hm;                    /* Hessian Misfit*/
     28:   Mat         Hr;                    /* Hessian Reg*/
     29:   Vec         d;                     /* RHS in least squares component $(1/2) * || F x - d ||_2^2$ */
     30:   Vec         workLeft[NWORKLEFT];   /* Workspace for temporary vec */
     31:   Vec         workRight[NWORKRIGHT]; /* Workspace for temporary vec */
     32:   NormType    p;
     33:   PetscRandom rctx;
     34:   PetscBool   taylor;                /* Flag to determine whether to run Taylor test or not */
     35:   PetscBool   use_admm;              /* Flag to determine whether to run Taylor test or not */
     36: }* UserCtx;

     38: static PetscErrorCode CreateRHS(UserCtx ctx)
     39: {

     43:   /* build the rhs d in ctx */
     44:   VecCreate(PETSC_COMM_WORLD,&(ctx->d));
     45:   VecSetSizes(ctx->d,PETSC_DECIDE,ctx->m);
     46:   VecSetFromOptions(ctx->d);
     47:   VecSetRandom(ctx->d,ctx->rctx);
     48:   return(0);
     49: }

     51: static PetscErrorCode CreateMatrix(UserCtx ctx)
     52: {
     53:   PetscInt       Istart,Iend,i,j,Ii,gridN,I_n, I_s, I_e, I_w;
     54: #if defined(PETSC_USE_LOG)
     55:   PetscLogStage  stage;
     56: #endif

     60:   /* build the matrix F in ctx */
     61:   MatCreate(PETSC_COMM_WORLD, &(ctx->F));
     62:   MatSetSizes(ctx->F,PETSC_DECIDE, PETSC_DECIDE, ctx->m, ctx->n);
     63:   MatSetType(ctx->F,MATAIJ); /* TODO: Decide specific SetType other than dummy*/
     64:   MatMPIAIJSetPreallocation(ctx->F, 5, NULL, 5, NULL); /*TODO: some number other than 5?*/
     65:   MatSeqAIJSetPreallocation(ctx->F, 5, NULL);
     66:   MatSetUp(ctx->F);
     67:   MatGetOwnershipRange(ctx->F,&Istart,&Iend);
     68:   PetscLogStageRegister("Assembly", &stage);
     69:   PetscLogStagePush(stage);

     71:   /* Set matrix elements in  2-D fiveopoint stencil format. */
     72:   if (!(ctx->matops)) {
     73:     if (ctx->m != ctx->n) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "Stencil matrix must be square");
     74:     gridN = (PetscInt) PetscSqrtReal((PetscReal) ctx->m);
     75:     if (gridN * gridN != ctx->m) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "Number of rows must be square");
     76:     for (Ii=Istart; Ii<Iend; Ii++) {
     77:       i   = Ii / gridN; j = Ii % gridN;
     78:       I_n = i * gridN + j + 1;
     79:       if (j + 1 >= gridN) I_n = -1;
     80:       I_s = i * gridN + j - 1;
     81:       if (j - 1 < 0) I_s = -1;
     82:       I_e = (i + 1) * gridN + j;
     83:       if (i + 1 >= gridN) I_e = -1;
     84:       I_w = (i - 1) * gridN + j;
     85:       if (i - 1 < 0) I_w = -1;
     86:       MatSetValue(ctx->F, Ii, Ii, 4., INSERT_VALUES);
     87:       MatSetValue(ctx->F, Ii, I_n, -1., INSERT_VALUES);
     88:       MatSetValue(ctx->F, Ii, I_s, -1., INSERT_VALUES);
     89:       MatSetValue(ctx->F, Ii, I_e, -1., INSERT_VALUES);
     90:       MatSetValue(ctx->F, Ii, I_w, -1., INSERT_VALUES);
     91:     }
     92:   } else {MatSetRandom(ctx->F, ctx->rctx);}
     93:   MatAssemblyBegin(ctx->F, MAT_FINAL_ASSEMBLY);
     94:   MatAssemblyEnd(ctx->F, MAT_FINAL_ASSEMBLY);
     95:   PetscLogStagePop();
     96:   /* Stencil matrix is symmetric. Setting symmetric flag for ICC/Cholesky preconditioner */
     97:   if (!(ctx->matops)) {
     98:     MatSetOption(ctx->F,MAT_SYMMETRIC,PETSC_TRUE);
     99:   }
    100:   MatTransposeMatMult(ctx->F,ctx->F, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &(ctx->W));
    101:   /* Setup Hessian Workspace in same shape as W */
    102:   MatDuplicate(ctx->W,MAT_DO_NOT_COPY_VALUES,&(ctx->Hm));
    103:   MatDuplicate(ctx->W,MAT_DO_NOT_COPY_VALUES,&(ctx->Hr));
    104:   return(0);
    105: }

    107: static PetscErrorCode SetupWorkspace(UserCtx ctx)
    108: {
    109:   PetscInt       i;

    113:   MatCreateVecs(ctx->F, &ctx->workLeft[0], &ctx->workRight[0]);
    114:   for (i=1; i<NWORKLEFT; i++) {
    115:     VecDuplicate(ctx->workLeft[0], &(ctx->workLeft[i]));
    116:   }
    117:   for (i=1; i<NWORKRIGHT; i++) {
    118:     VecDuplicate(ctx->workRight[0], &(ctx->workRight[i]));
    119:   }
    120:   return(0);
    121: }

    123: static PetscErrorCode ConfigureContext(UserCtx ctx)
    124: {

    128:   ctx->m        = 16;
    129:   ctx->n        = 16;
    130:   ctx->eps      = 1.e-3;
    131:   ctx->abstol   = 1.e-4;
    132:   ctx->reltol   = 1.e-2;
    133:   ctx->hStart   = 1.;
    134:   ctx->hMin     = 1.e-3;
    135:   ctx->hFactor  = 0.5;
    136:   ctx->alpha    = 1.;
    137:   ctx->mu       = 1.0;
    138:   ctx->matops   = 0;
    139:   ctx->iter     = 10;
    140:   ctx->p        = NORM_2;
    141:   ctx->taylor   = PETSC_TRUE;
    142:   ctx->use_admm = PETSC_FALSE;
    143:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Configure separable objection example", "ex4.c");
    144:   PetscOptionsInt("-m", "The row dimension of matrix F", "ex4.c", ctx->m, &(ctx->m), NULL);
    145:   PetscOptionsInt("-n", "The column dimension of matrix F", "ex4.c", ctx->n, &(ctx->n), NULL);
    146:   PetscOptionsInt("-matrix_format","Decide format of F matrix. 0 for stencil, 1 for random", "ex4.c", ctx->matops, &(ctx->matops), NULL);
    147:   PetscOptionsInt("-iter","Iteration number ADMM", "ex4.c", ctx->iter, &(ctx->iter), NULL);
    148:   PetscOptionsReal("-alpha", "The regularization multiplier. 1 default", "ex4.c", ctx->alpha, &(ctx->alpha), NULL);
    149:   PetscOptionsReal("-epsilon", "The small constant added to |x_i| in the denominator to approximate the gradient of ||x||_1", "ex4.c", ctx->eps, &(ctx->eps), NULL);
    150:   PetscOptionsReal("-mu", "The augmented lagrangian multiplier in ADMM", "ex4.c", ctx->mu, &(ctx->mu), NULL);
    151:   PetscOptionsReal("-hStart", "Taylor test starting point. 1 default.", "ex4.c", ctx->hStart, &(ctx->hStart), NULL);
    152:   PetscOptionsReal("-hFactor", "Taylor test multiplier factor. 0.5 default", "ex4.c", ctx->hFactor, &(ctx->hFactor), NULL);
    153:   PetscOptionsReal("-hMin", "Taylor test ending condition. 1.e-3 default", "ex4.c", ctx->hMin, &(ctx->hMin), NULL);
    154:   PetscOptionsReal("-abstol", "Absolute stopping criterion for ADMM", "ex4.c", ctx->abstol, &(ctx->abstol), NULL);
    155:   PetscOptionsReal("-reltol", "Relative stopping criterion for ADMM", "ex4.c", ctx->reltol, &(ctx->reltol), NULL);
    156:   PetscOptionsBool("-taylor","Flag for Taylor test. Default is true.", "ex4.c", ctx->taylor, &(ctx->taylor), NULL);
    157:   PetscOptionsBool("-use_admm","Use the ADMM solver in this example.", "ex4.c", ctx->use_admm, &(ctx->use_admm), NULL);
    158:   PetscOptionsEnum("-p","Norm type.", "ex4.c", NormTypes, (PetscEnum)ctx->p, (PetscEnum *) &(ctx->p), NULL);
    159:   PetscOptionsEnd();
    160:   /* Creating random ctx */
    161:   PetscRandomCreate(PETSC_COMM_WORLD,&(ctx->rctx));
    162:   PetscRandomSetFromOptions(ctx->rctx);
    163:   CreateMatrix(ctx);
    164:   CreateRHS(ctx);
    165:   SetupWorkspace(ctx);
    166:   return(0);
    167: }

    169: static PetscErrorCode DestroyContext(UserCtx *ctx)
    170: {
    171:   PetscInt       i;

    175:   MatDestroy(&((*ctx)->F));
    176:   MatDestroy(&((*ctx)->W));
    177:   MatDestroy(&((*ctx)->Hm));
    178:   MatDestroy(&((*ctx)->Hr));
    179:   VecDestroy(&((*ctx)->d));
    180:   for (i=0; i<NWORKLEFT; i++) {
    181:     VecDestroy(&((*ctx)->workLeft[i]));
    182:   }
    183:   for (i=0; i<NWORKRIGHT; i++) {
    184:     VecDestroy(&((*ctx)->workRight[i]));
    185:   }
    186:   PetscRandomDestroy(&((*ctx)->rctx));
    187:   PetscFree(*ctx);
    188:   return(0);
    189: }

    191: /* compute (1/2) * ||F x - d||^2 */
    192: static PetscErrorCode ObjectiveMisfit(Tao tao, Vec x, PetscReal *J, void *_ctx)
    193: {
    194:   UserCtx        ctx = (UserCtx) _ctx;
    196:   Vec            y;

    199:   y    = ctx->workLeft[0];
    200:   MatMult(ctx->F, x, y);
    201:   VecAXPY(y, -1., ctx->d);
    202:   VecDot(y, y, J);
    203:   *J  *= 0.5;
    204:   return(0);
    205: }

    207: /* compute V = FTFx - FTd */
    208: static PetscErrorCode GradientMisfit(Tao tao, Vec x, Vec V, void *_ctx)
    209: {
    210:   UserCtx        ctx = (UserCtx) _ctx;
    212:   Vec            FTFx, FTd;

    215:   /* work1 is A^T Ax, work2 is Ab, W is A^T A*/
    216:   FTFx = ctx->workRight[0];
    217:   FTd  = ctx->workRight[1];
    218:   MatMult(ctx->W,x,FTFx);
    219:   MatMultTranspose(ctx->F, ctx->d, FTd);
    220:   VecWAXPY(V, -1., FTd, FTFx);
    221:   return(0);
    222: }

    224: /* returns FTF */
    225: static PetscErrorCode HessianMisfit(Tao tao, Vec x, Mat H, Mat Hpre, void *_ctx)
    226: {
    227:   UserCtx        ctx = (UserCtx) _ctx;

    231:   if (H != ctx->W) {MatCopy(ctx->W, H, DIFFERENT_NONZERO_PATTERN);}
    232:   if (Hpre != ctx->W) {MatCopy(ctx->W, Hpre, DIFFERENT_NONZERO_PATTERN);}
    233:   return(0);
    234: }

    236: /* computes augment Lagrangian objective (with scaled dual):
    237:  * 0.5 * ||F x - d||^2  + 0.5 * mu ||x - z + u||^2 */
    238: static PetscErrorCode ObjectiveMisfitADMM(Tao tao, Vec x, PetscReal *J, void *_ctx)
    239: {
    240:   UserCtx        ctx = (UserCtx) _ctx;
    241:   PetscReal      mu, workNorm, misfit;
    242:   Vec            z, u, temp;

    246:   mu   = ctx->mu;
    247:   z    = ctx->workRight[5];
    248:   u    = ctx->workRight[6];
    249:   temp = ctx->workRight[10];
    250:   /* misfit = f(x) */
    251:   ObjectiveMisfit(tao, x, &misfit, _ctx);
    252:   VecCopy(x,temp);
    253:   /* temp = x - z + u */
    254:   VecAXPBYPCZ(temp,-1.,1.,1.,z,u);
    255:   /* workNorm = ||x - z + u||^2 */
    256:   VecDot(temp, temp, &workNorm);
    257:   /* augment Lagrangian objective (with scaled dual): f(x) + 0.5 * mu ||x -z + u||^2 */
    258:   *J = misfit + 0.5 * mu * workNorm;
    259:   return(0);
    260: }

    262: /* computes FTFx - FTd  mu*(x - z + u) */
    263: static PetscErrorCode GradientMisfitADMM(Tao tao, Vec x, Vec V, void *_ctx)
    264: {
    265:   UserCtx        ctx = (UserCtx) _ctx;
    266:   PetscReal      mu;
    267:   Vec            z, u, temp;

    271:   mu   = ctx->mu;
    272:   z    = ctx->workRight[5];
    273:   u    = ctx->workRight[6];
    274:   temp = ctx->workRight[10];
    275:   GradientMisfit(tao, x, V, _ctx);
    276:   VecCopy(x, temp);
    277:   /* temp = x - z + u */
    278:   VecAXPBYPCZ(temp,-1.,1.,1.,z,u);
    279:   /* V =  FTFx - FTd  mu*(x - z + u) */
    280:   VecAXPY(V, mu, temp);
    281:   return(0);
    282: }

    284: /* returns FTF + diag(mu) */
    285: static PetscErrorCode HessianMisfitADMM(Tao tao, Vec x, Mat H, Mat Hpre, void *_ctx)
    286: {
    287:   UserCtx        ctx = (UserCtx) _ctx;

    291:   MatCopy(ctx->W, H, DIFFERENT_NONZERO_PATTERN);
    292:   MatShift(H, ctx->mu);
    293:   if (Hpre != H) {
    294:     MatCopy(H, Hpre, DIFFERENT_NONZERO_PATTERN);
    295:   }
    296:   return(0);
    297: }

    299: /* computes || x ||_p (mult by 0.5 in case of NORM_2) */
    300: static PetscErrorCode ObjectiveRegularization(Tao tao, Vec x, PetscReal *J, void *_ctx)
    301: {
    302:   UserCtx        ctx = (UserCtx) _ctx;
    303:   PetscReal      norm;

    307:   *J = 0;
    308:   VecNorm (x, ctx->p, &norm);
    309:   if (ctx->p == NORM_2) norm = 0.5 * norm * norm;
    310:   *J = ctx->alpha * norm;
    311:   return(0);
    312: }

    314: /* NORM_2 Case: return x
    315:  * NORM_1 Case: x/(|x| + eps)
    316:  * Else: TODO */
    317: static PetscErrorCode GradientRegularization(Tao tao, Vec x, Vec V, void *_ctx)
    318: {
    319:   UserCtx        ctx = (UserCtx) _ctx;
    321:   PetscReal      eps = ctx->eps;

    324:   if (ctx->p == NORM_2) {
    325:     VecCopy(x, V);
    326:   } else if (ctx->p == NORM_1) {
    327:     VecCopy(x, ctx->workRight[1]);
    328:     VecAbs(ctx->workRight[1]);
    329:     VecShift(ctx->workRight[1], eps);
    330:     VecPointwiseDivide(V, x, ctx->workRight[1]);
    331:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_OUTOFRANGE, "Example only works for NORM_1 and NORM_2");
    332:   return(0);
    333: }

    335: /* NORM_2 Case: returns diag(mu)
    336:  * NORM_1 Case: diag(mu* 1/sqrt(x_i^2 + eps) * (1 - x_i^2/ABS(x_i^2+eps)))  */
    337: static PetscErrorCode HessianRegularization(Tao tao, Vec x, Mat H, Mat Hpre, void *_ctx)
    338: {
    339:   UserCtx        ctx = (UserCtx) _ctx;
    340:   PetscReal      eps = ctx->eps;
    341:   Vec            copy1,copy2,copy3;

    345:   if (ctx->p == NORM_2) {
    346:     /* Identity matrix scaled by mu */
    347:     MatZeroEntries(H);
    348:     MatShift(H,ctx->mu);
    349:     if (Hpre != H) {
    350:       MatZeroEntries(Hpre);
    351:       MatShift(Hpre,ctx->mu);
    352:     }
    353:   } else if (ctx->p == NORM_1) {
    354:     /* 1/sqrt(x_i^2 + eps) * (1 - x_i^2/ABS(x_i^2+eps)) */
    355:     copy1 = ctx->workRight[1];
    356:     copy2 = ctx->workRight[2];
    357:     copy3 = ctx->workRight[3];
    358:     /* copy1 : 1/sqrt(x_i^2 + eps) */
    359:     VecCopy(x, copy1);
    360:     VecPow(copy1,2);
    361:     VecShift(copy1, eps);
    362:     VecSqrtAbs(copy1);
    363:     VecReciprocal(copy1);
    364:     /* copy2:  x_i^2.*/
    365:     VecCopy(x,copy2);
    366:     VecPow(copy2,2);
    367:     /* copy3: abs(x_i^2 + eps) */
    368:     VecCopy(x,copy3);
    369:     VecPow(copy3,2);
    370:     VecShift(copy3, eps);
    371:     VecAbs(copy3);
    372:     /* copy2: 1 - x_i^2/abs(x_i^2 + eps) */
    373:     VecPointwiseDivide(copy2, copy2,copy3);
    374:     VecScale(copy2, -1.);
    375:     VecShift(copy2, 1.);
    376:     VecAXPY(copy1,1.,copy2);
    377:     VecScale(copy1, ctx->mu);
    378:     MatZeroEntries(H);
    379:     MatDiagonalSet(H, copy1,INSERT_VALUES);
    380:     if (Hpre != H) {
    381:       MatZeroEntries(Hpre);
    382:       MatDiagonalSet(Hpre, copy1,INSERT_VALUES);
    383:     }
    384:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_OUTOFRANGE, "Example only works for NORM_1 and NORM_2");
    385:   return(0);
    386: }

    388: /* NORM_2 Case: 0.5 || x ||_2 + 0.5 * mu * ||x + u - z||^2
    389:  * Else : || x ||_2 + 0.5 * mu * ||x + u - z||^2 */
    390: static PetscErrorCode ObjectiveRegularizationADMM(Tao tao, Vec z, PetscReal *J, void *_ctx)
    391: {
    392:   UserCtx        ctx = (UserCtx) _ctx;
    393:   PetscReal      mu, workNorm, reg;
    394:   Vec            x, u, temp;

    398:   mu   = ctx->mu;
    399:   x    = ctx->workRight[4];
    400:   u    = ctx->workRight[6];
    401:   temp = ctx->workRight[10];
    402:   ObjectiveRegularization(tao, z, &reg, _ctx);
    403:   VecCopy(z,temp);
    404:   /* temp = x + u -z */
    405:   VecAXPBYPCZ(temp,1.,1.,-1.,x,u);
    406:   /* workNorm = ||x + u - z ||^2 */
    407:   VecDot(temp, temp, &workNorm);
    408:   *J   = reg + 0.5 * mu * workNorm;
    409:   return(0);
    410: }

    413: /* NORM_2 Case: x - mu*(x + u - z)
    414:  * NORM_1 Case: x/(|x| + eps) - mu*(x + u - z)
    415:  * Else: TODO */
    416: static PetscErrorCode GradientRegularizationADMM(Tao tao, Vec z, Vec V, void *_ctx)
    417: {
    418:   UserCtx        ctx = (UserCtx) _ctx;
    419:   PetscReal      mu;
    420:   Vec            x, u, temp;

    424:   mu   = ctx->mu;
    425:   x    = ctx->workRight[4];
    426:   u    = ctx->workRight[6];
    427:   temp = ctx->workRight[10];
    428:   GradientRegularization(tao, z, V, _ctx);
    429:   VecCopy(z, temp);
    430:   /* temp = x + u -z */
    431:   VecAXPBYPCZ(temp,1.,1.,-1.,x,u);
    432:   VecAXPY(V, -mu, temp);
    433:   return(0);
    434: }

    436: /* NORM_2 Case: returns diag(mu)
    437:  * NORM_1 Case: FTF + diag(mu) */
    438: static PetscErrorCode HessianRegularizationADMM(Tao tao, Vec x, Mat H, Mat Hpre, void *_ctx)
    439: {
    440:   UserCtx        ctx = (UserCtx) _ctx;

    444:   if (ctx->p == NORM_2) {
    445:     /* Identity matrix scaled by mu */
    446:     MatZeroEntries(H);
    447:     MatShift(H,ctx->mu);
    448:     if (Hpre != H) {
    449:       MatZeroEntries(Hpre);
    450:       MatShift(Hpre,ctx->mu);
    451:     }
    452:   } else if (ctx->p == NORM_1) {
    453:     HessianMisfit(tao, x, H, Hpre, (void*) ctx);
    454:     MatShift(H, ctx->mu);
    455:     if (Hpre != H) {MatShift(Hpre, ctx->mu);}
    456:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_OUTOFRANGE, "Example only works for NORM_1 and NORM_2");
    457:   return(0);
    458: }

    460: /* NORM_2 Case : (1/2) * ||F x - d||^2 + 0.5 * || x ||_p
    461: *  NORM_1 Case : (1/2) * ||F x - d||^2 + || x ||_p */
    462: static PetscErrorCode ObjectiveComplete(Tao tao, Vec x, PetscReal *J, void *ctx)
    463: {
    464:   PetscReal      Jm, Jr;

    468:   ObjectiveMisfit(tao, x, &Jm, ctx);
    469:   ObjectiveRegularization(tao, x, &Jr, ctx);
    470:   *J   = Jm + Jr;
    471:   return(0);
    472: }

    474: /* NORM_2 Case: FTFx - FTd + x
    475:  * NORM_1 Case: FTFx - FTd + x/(|x| + eps) */
    476: static PetscErrorCode GradientComplete(Tao tao, Vec x, Vec V, void *ctx)
    477: {
    478:   UserCtx        cntx = (UserCtx) ctx;

    482:   GradientMisfit(tao, x, cntx->workRight[2], ctx);
    483:   GradientRegularization(tao, x, cntx->workRight[3], ctx);
    484:   VecWAXPY(V,1,cntx->workRight[2],cntx->workRight[3]);
    485:   return(0);
    486: }

    488: /* NORM_2 Case: diag(mu) + FTF
    489:  * NORM_1 Case: diag(mu* 1/sqrt(x_i^2 + eps) * (1 - x_i^2/ABS(x_i^2+eps))) + FTF  */
    490: static PetscErrorCode HessianComplete(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
    491: {
    492:   Mat            tempH;

    496:   MatDuplicate(H, MAT_SHARE_NONZERO_PATTERN, &tempH);
    497:   HessianMisfit(tao, x, H, H, ctx);
    498:   HessianRegularization(tao, x, tempH, tempH, ctx);
    499:   MatAXPY(H, 1., tempH, DIFFERENT_NONZERO_PATTERN);
    500:   if (Hpre != H) {
    501:     MatCopy(H, Hpre, DIFFERENT_NONZERO_PATTERN);
    502:   }
    503:   MatDestroy(&tempH);
    504:   return(0);
    505: }

    507: static PetscErrorCode TaoSolveADMM(UserCtx ctx,  Vec x)
    508: {
    510:   PetscInt       i;
    511:   PetscReal      u_norm, r_norm, s_norm, primal, dual, x_norm, z_norm;
    512:   Tao            tao1,tao2;
    513:   Vec            xk,z,u,diff,zold,zdiff,temp;
    514:   PetscReal      mu;

    517:   xk    = ctx->workRight[4];
    518:   z     = ctx->workRight[5];
    519:   u     = ctx->workRight[6];
    520:   diff  = ctx->workRight[7];
    521:   zold  = ctx->workRight[8];
    522:   zdiff = ctx->workRight[9];
    523:   temp  = ctx->workRight[11];
    524:   mu    = ctx->mu;
    525:   VecSet(u, 0.);
    526:   TaoCreate(PETSC_COMM_WORLD, &tao1);
    527:   TaoSetType(tao1,TAONLS);
    528:   TaoSetObjectiveRoutine(tao1, ObjectiveMisfitADMM, (void*) ctx);
    529:   TaoSetGradientRoutine(tao1, GradientMisfitADMM, (void*) ctx);
    530:   TaoSetHessianRoutine(tao1, ctx->Hm, ctx->Hm, HessianMisfitADMM, (void*) ctx);
    531:   VecSet(xk, 0.);
    532:   TaoSetInitialVector(tao1, xk);
    533:   TaoSetOptionsPrefix(tao1, "misfit_");
    534:   TaoSetFromOptions(tao1);
    535:   TaoCreate(PETSC_COMM_WORLD, &tao2);
    536:   if (ctx->p == NORM_2) {
    537:     TaoSetType(tao2,TAONLS);
    538:     TaoSetObjectiveRoutine(tao2, ObjectiveRegularizationADMM, (void*) ctx);
    539:     TaoSetGradientRoutine(tao2, GradientRegularizationADMM, (void*) ctx);
    540:     TaoSetHessianRoutine(tao2, ctx->Hr, ctx->Hr, HessianRegularizationADMM, (void*) ctx);
    541:   }
    542:   VecSet(z, 0.);
    543:   TaoSetInitialVector(tao2, z);
    544:   TaoSetOptionsPrefix(tao2, "reg_");
    545:   TaoSetFromOptions(tao2);

    547:   for (i=0; i<ctx->iter; i++) {
    548:     VecCopy(z,zold);
    549:     TaoSolve(tao1); /* Updates xk */
    550:     if (ctx->p == NORM_1){
    551:       VecWAXPY(temp,1.,xk,u);
    552:       TaoSoftThreshold(temp,-ctx->alpha/mu,ctx->alpha/mu,z);
    553:     } else {
    554:       TaoSolve(tao2); /* Update zk */
    555:     }
    556:     /* u = u + xk -z */
    557:     VecAXPBYPCZ(u,1.,-1.,1.,xk,z);
    558:     /* r_norm : norm(x-z) */
    559:     VecWAXPY(diff,-1.,z,xk);
    560:     VecNorm(diff,NORM_2,&r_norm);
    561:     /* s_norm : norm(-mu(z-zold)) */
    562:     VecWAXPY(zdiff, -1.,zold,z);
    563:     VecNorm(zdiff,NORM_2,&s_norm);
    564:     s_norm = s_norm * mu;
    565:     /* primal : sqrt(n)*ABSTOL + RELTOL*max(norm(x), norm(-z))*/
    566:     VecNorm(xk,NORM_2,&x_norm);
    567:     VecNorm(z,NORM_2,&z_norm);
    568:     primal = PetscSqrtReal(ctx->n)*ctx->abstol + ctx->reltol*PetscMax(x_norm,z_norm);
    569:     /* Duality : sqrt(n)*ABSTOL + RELTOL*norm(mu*u)*/
    570:     VecNorm(u,NORM_2,&u_norm);
    571:     dual = PetscSqrtReal(ctx->n)*ctx->abstol + ctx->reltol*u_norm*mu;
    572:     PetscPrintf(PetscObjectComm((PetscObject)tao1),"Iter %D : ||x-z||: %g, mu*||z-zold||: %g\n", i, (double) r_norm, (double) s_norm);
    573:     if (r_norm < primal && s_norm < dual) break;
    574:   }
    575:   VecCopy(xk, x);
    576:   TaoDestroy(&tao1);
    577:   TaoDestroy(&tao2);
    578:   return(0);
    579: }

    581: /* Second order Taylor remainder convergence test */
    582: static PetscErrorCode TaylorTest(UserCtx ctx, Tao tao, Vec x, PetscReal *C)
    583: {
    584:   PetscReal      h,J,temp;
    585:   PetscInt       i,j;
    586:   PetscInt       numValues;
    587:   PetscReal      Jx,Jxhat_comp,Jxhat_pred;
    588:   PetscReal      *Js, *hs;
    589:   PetscReal      gdotdx;
    590:   PetscReal      minrate = PETSC_MAX_REAL;
    591:   MPI_Comm       comm = PetscObjectComm((PetscObject)x);
    592:   Vec            g, dx, xhat;

    596:   VecDuplicate(x, &g);
    597:   VecDuplicate(x, &xhat);
    598:   /* choose a perturbation direction */
    599:   VecDuplicate(x, &dx);
    600:   VecSetRandom(dx,ctx->rctx);
    601:   /* evaluate objective at x: J(x) */
    602:   TaoComputeObjective(tao, x, &Jx);
    603:   /* evaluate gradient at x, save in vector g */
    604:   TaoComputeGradient(tao, x, g);
    605:   VecDot(g, dx, &gdotdx);

    607:   for (numValues=0, h=ctx->hStart; h>=ctx->hMin; h*=ctx->hFactor) numValues++;
    608:   PetscCalloc2(numValues, &Js, numValues, &hs);
    609:   for (i=0, h=ctx->hStart; h>=ctx->hMin; h*=ctx->hFactor, i++) {
    610:     VecWAXPY(xhat, h, dx, x);
    611:     TaoComputeObjective(tao, xhat, &Jxhat_comp);
    612:     /* J(\hat(x)) \approx J(x) + g^T (xhat - x) = J(x) + h * g^T dx */
    613:     Jxhat_pred = Jx + h * gdotdx;
    614:     /* Vector to dJdm scalar? Dot?*/
    615:     J     = PetscAbsReal(Jxhat_comp - Jxhat_pred);
    616:     PetscPrintf (comm, "J(xhat): %g, predicted: %g, diff %g\n", (double) Jxhat_comp,(double) Jxhat_pred, (double) J);
    617:     Js[i] = J;
    618:     hs[i] = h;
    619:   }
    620:   for (j=1; j<numValues; j++) {
    621:     temp    = PetscLogReal(Js[j]/Js[j-1]) / PetscLogReal (hs[j]/hs[j-1]);
    622:     PetscPrintf (comm, "Convergence rate step %D: %g\n", j-1, (double) temp);
    623:     minrate = PetscMin(minrate, temp);
    624:   }
    625:   /* If O is not ~2, then the test is wrong */
    626:   PetscFree2(Js, hs);
    627:   *C   = minrate;
    628:   VecDestroy(&dx);
    629:   VecDestroy(&xhat);
    630:   VecDestroy(&g);
    631:   return(0);
    632: }

    634: int main(int argc, char ** argv)
    635: {
    636:   UserCtx        ctx;
    637:   Tao            tao;
    638:   Vec            x;
    639:   Mat            H;

    642:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
    643:   PetscNew(&ctx);
    644:   ConfigureContext(ctx);
    645:   /* Define two functions that could pass as objectives to TaoSetObjectiveRoutine(): one
    646:    * for the misfit component, and one for the regularization component */
    647:   /* ObjectiveMisfit() and ObjectiveRegularization() */

    649:   /* Define a single function that calls both components adds them together: the complete objective,
    650:    * in the absence of a Tao implementation that handles separability */
    651:   /* ObjectiveComplete() */
    652:   TaoCreate(PETSC_COMM_WORLD, &tao);
    653:   TaoSetType(tao,TAONM);
    654:   TaoSetObjectiveRoutine(tao, ObjectiveComplete, (void*) ctx);
    655:   TaoSetGradientRoutine(tao, GradientComplete, (void*) ctx);
    656:   MatDuplicate(ctx->W, MAT_SHARE_NONZERO_PATTERN, &H);
    657:   TaoSetHessianRoutine(tao, H, H, HessianComplete, (void*) ctx);
    658:   MatCreateVecs(ctx->F, NULL, &x);
    659:   VecSet(x, 0.);
    660:   TaoSetInitialVector(tao, x);
    661:   TaoSetFromOptions(tao);
    662:   if (ctx->use_admm) {
    663:     TaoSolveADMM(ctx,x);
    664:   } else {TaoSolve(tao);}
    665:   /* examine solution */
    666:   VecViewFromOptions(x, NULL, "-view_sol");
    667:   if (ctx->taylor) {
    668:     PetscReal rate;
    669:     TaylorTest(ctx, tao, x, &rate);
    670:   }
    671:   MatDestroy(&H);
    672:   TaoDestroy(&tao);
    673:   VecDestroy(&x);
    674:   DestroyContext(&ctx);
    675:   PetscFinalize();
    676:   return ierr;
    677: }

    679: /*TEST

    681:   build:
    682:     requires: !complex

    684:   test:
    685:     suffix: 0
    686:     args:

    688:   test:
    689:     suffix: l1_1
    690:     args: -p 1 -tao_type lmvm -alpha 1. -epsilon 1.e-7 -m 64 -n 64 -view_sol -matrix_format 1

    692:   test:
    693:     suffix: hessian_1
    694:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 1 -tao_type nls -tao_nls_ksp_monitor

    696:   test:
    697:     suffix: hessian_2
    698:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 2 -tao_type nls -tao_nls_ksp_monitor

    700:   test:
    701:     suffix: nm_1
    702:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 1 -tao_type nm -tao_max_it 50

    704:   test:
    705:     suffix: nm_2
    706:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 2 -tao_type nm -tao_max_it 50

    708:   test:
    709:     suffix: lmvm_1
    710:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 1 -tao_type lmvm -tao_max_it 40

    712:   test:
    713:     suffix: lmvm_2
    714:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 2 -tao_type lmvm -tao_max_it 15

    716:   test:
    717:     suffix: soft_threshold_admm_1
    718:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 1 -use_admm

    720:   test:
    721:     suffix: hessian_admm_1
    722:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 1 -use_admm -reg_tao_type nls -misfit_tao_type nls

    724:   test:
    725:     suffix: hessian_admm_2
    726:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 2 -use_admm -reg_tao_type nls -misfit_tao_type nls

    728:   test:
    729:     suffix: nm_admm_1
    730:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 1 -use_admm -reg_tao_type nm -misfit_tao_type nm

    732:   test:
    733:     suffix: nm_admm_2
    734:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 2 -use_admm -reg_tao_type nm -misfit_tao_type nm -iter 7

    736:   test:
    737:     suffix: lmvm_admm_1
    738:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 1 -use_admm -reg_tao_type lmvm -misfit_tao_type lmvm

    740:   test:
    741:     suffix: lmvm_admm_2
    742:     args: -matrix_format 1 -m 100 -n 100 -tao_monitor -p 2 -use_admm -reg_tao_type lmvm -misfit_tao_type lmvm

    744: TEST*/