Actual source code: taoterml1.c

      1: #include <petsc/private/taoimpl.h>

      3: typedef struct _n_TaoTerm_L1 TaoTerm_L1;

      5: struct _n_TaoTerm_L1 {
      6:   PetscReal        epsilon;
      7:   PetscBool        epsilon_warning;
      8:   Vec              diff; // caches $x - p$
      9:   Vec              d;    // caches $d_i = sqrt((diff_i)^2 + epsilon^2)
     10:   PetscReal        d_epsilon;
     11:   PetscReal        d_sum;
     12:   PetscObjectId    d_x_id, d_p_id;
     13:   PetscObjectState d_x_state, d_p_state;
     14:   Vec              diag;
     15:   PetscObjectId    diag_x_id, diag_p_id;
     16:   PetscObjectState diag_x_state, diag_p_state;
     17:   PetscReal        diag_epsilon;
     18: };

     20: static PetscErrorCode TaoTermDestroy_L1(TaoTerm term)
     21: {
     22:   TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;

     24:   PetscFunctionBegin;
     25:   PetscCall(VecDestroy(&l1->diff));
     26:   PetscCall(VecDestroy(&l1->d));
     27:   PetscCall(VecDestroy(&l1->diag));
     28:   PetscCall(PetscFree(l1));
     29:   term->data = NULL;
     30:   PetscCall(TaoTermDestroy_ElementwiseDivergence_Internal(term));
     31:   PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermL1SetEpsilon_C", NULL));
     32:   PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermL1GetEpsilon_C", NULL));
     33:   PetscFunctionReturn(PETSC_SUCCESS);
     34: }

     36: static PetscErrorCode TaoTermL1ComputeData(TaoTerm term, Vec x, Vec params, Vec *_diff, Vec *d)
     37: {
     38:   TaoTerm_L1      *l1 = (TaoTerm_L1 *)term->data;
     39:   PetscObjectId    x_id, p_id       = 0;
     40:   PetscObjectState x_state, p_state = 0;
     41:   Vec              diff = x;

     43:   PetscFunctionBegin;
     44:   PetscCall(PetscObjectGetId((PetscObject)x, &x_id));
     45:   PetscCall(PetscObjectStateGet((PetscObject)x, &x_state));
     46:   if (params) {
     47:     PetscCall(PetscObjectGetId((PetscObject)params, &p_id));
     48:     PetscCall(PetscObjectStateGet((PetscObject)params, &p_state));
     49:   }
     50:   if (l1->d_x_id != x_id || l1->d_x_state != x_state || l1->d_p_id != p_id || l1->d_p_state != p_state || l1->d_epsilon != l1->epsilon) {
     51:     l1->d_x_id    = x_id;
     52:     l1->d_x_state = x_state;
     53:     l1->d_p_id    = p_id;
     54:     l1->d_p_state = p_state;
     55:     l1->d_epsilon = l1->epsilon;
     56:     diff          = x;
     57:     if (params) {
     58:       PetscCall(VecIfNotCongruentGetSameLayoutVec(x, &l1->diff));
     59:       PetscCall(VecWAXPY(l1->diff, -1.0, params, x));
     60:       diff = l1->diff;
     61:     }
     62:     if (l1->epsilon != 0.0) {
     63:       PetscCall(VecIfNotCongruentGetSameLayoutVec(x, &l1->d));
     64:       PetscCall(VecPointwiseMult(l1->d, diff, diff));
     65:       PetscCall(VecShift(l1->d, l1->epsilon * l1->epsilon));
     66:       PetscCall(VecSqrtAbs(l1->d));
     67:     }
     68:   }
     69:   if (params) diff = l1->diff;
     70:   *_diff = diff;
     71:   *d     = l1->d;
     72:   PetscFunctionReturn(PETSC_SUCCESS);
     73: }

     75: static PetscErrorCode TaoTermL1ComputeDiag(TaoTerm term, Vec x, Vec params, Vec *diag)
     76: {
     77:   Vec              diff, d;
     78:   TaoTerm_L1      *l1 = (TaoTerm_L1 *)term->data;
     79:   PetscObjectId    x_id, p_id       = 0;
     80:   PetscObjectState x_state, p_state = 0;

     82:   PetscFunctionBegin;
     83:   PetscCall(TaoTermL1ComputeData(term, x, params, &diff, &d));
     84:   PetscCall(PetscObjectGetId((PetscObject)x, &x_id));
     85:   PetscCall(PetscObjectStateGet((PetscObject)x, &x_state));
     86:   if (params) {
     87:     PetscCall(PetscObjectGetId((PetscObject)params, &p_id));
     88:     PetscCall(PetscObjectStateGet((PetscObject)params, &p_state));
     89:   }
     90:   if (l1->diag_x_id != x_id || l1->diag_x_state != x_state || l1->diag_p_id != p_id || l1->diag_p_state != p_state || l1->diag_epsilon != l1->epsilon) {
     91:     l1->diag_x_id    = x_id;
     92:     l1->diag_x_state = x_state;
     93:     l1->diag_p_id    = p_id;
     94:     l1->diag_p_state = p_state;
     95:     l1->diag_epsilon = l1->epsilon;
     96:     if (l1->epsilon != 0.0) {
     97:       PetscCall(VecIfNotCongruentGetSameLayoutVec(x, &l1->diag));
     98:       PetscCall(VecCopy(d, l1->diag));
     99:       PetscCall(VecPointwiseMult(l1->diag, l1->diag, d));
    100:       PetscCall(VecPointwiseMult(l1->diag, l1->diag, d));
    101:       PetscCall(VecReciprocal(l1->diag));
    102:       PetscCall(VecScale(l1->diag, l1->epsilon * l1->epsilon));
    103:     }
    104:   }
    105:   *diag = l1->diag;
    106:   PetscFunctionReturn(PETSC_SUCCESS);
    107: }

    109: static PetscErrorCode TaoTermComputeObjective_L1_Internal(TaoTerm term, Vec diff, Vec d, PetscReal *value)
    110: {
    111:   TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;

    113:   PetscFunctionBegin;
    114:   if (l1->epsilon == 0.0) {
    115:     PetscCall(VecNorm(diff, NORM_1, value));
    116:   } else {
    117:     PetscScalar sum;
    118:     PetscInt    n;

    120:     PetscCall(VecGetSize(d, &n));
    121:     PetscCall(VecSum(d, &sum));
    122:     *value = PetscRealPart(sum) - n * l1->epsilon;
    123:   }
    124:   PetscFunctionReturn(PETSC_SUCCESS);
    125: }

    127: static PetscErrorCode TaoTermComputeObjective_L1(TaoTerm term, Vec x, Vec params, PetscReal *value)
    128: {
    129:   Vec diff, d;

    131:   PetscFunctionBegin;
    132:   PetscCall(TaoTermL1ComputeData(term, x, params, &diff, &d));
    133:   PetscCall(TaoTermComputeObjective_L1_Internal(term, diff, d, value));
    134:   PetscFunctionReturn(PETSC_SUCCESS);
    135: }

    137: static PetscErrorCode TaoTermL1DerivativeCheck(TaoTerm term)
    138: {
    139:   TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;

    141:   PetscFunctionBegin;
    142:   if (l1->epsilon_warning == PETSC_FALSE) {
    143:     l1->epsilon_warning = PETSC_TRUE;
    144:     PetscCall(PetscInfo(term, "%s: Asking for derivatives of l1 norm, which is not smooth.  Consider smoothing the TaoTerm with TaoTermL1SetEpsilon() or using a derivative-free Tao solver\n", ((PetscObject)term)->prefix));
    145:   }
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

    149: static PetscErrorCode TaoTermComputeGradient_L1_Internal(TaoTerm term, Vec diff, Vec d, Vec g)
    150: {
    151:   TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;

    153:   PetscFunctionBegin;
    154:   if (l1->epsilon == 0.0) {
    155:     PetscCall(TaoTermL1DerivativeCheck(term));
    156:     PetscCall(VecPointwiseSign(g, diff, VEC_SIGN_ZERO_TO_ZERO));
    157:   } else {
    158:     PetscCall(VecPointwiseDivide(g, diff, d));
    159:   }
    160:   PetscFunctionReturn(PETSC_SUCCESS);
    161: }

    163: static PetscErrorCode TaoTermComputeGradient_L1(TaoTerm term, Vec x, Vec params, Vec g)
    164: {
    165:   Vec diff, d;

    167:   PetscFunctionBegin;
    168:   PetscCall(TaoTermL1ComputeData(term, x, params, &diff, &d));
    169:   PetscCall(TaoTermComputeGradient_L1_Internal(term, diff, d, g));
    170:   PetscFunctionReturn(PETSC_SUCCESS);
    171: }

    173: static PetscErrorCode TaoTermComputeObjectiveAndGradient_L1(TaoTerm term, Vec x, Vec params, PetscReal *value, Vec g)
    174: {
    175:   Vec diff, d;

    177:   PetscFunctionBegin;
    178:   PetscCall(TaoTermL1ComputeData(term, x, params, &diff, &d));
    179:   PetscCall(TaoTermComputeObjective_L1_Internal(term, diff, d, value));
    180:   PetscCall(TaoTermComputeGradient_L1_Internal(term, diff, d, g));
    181:   PetscFunctionReturn(PETSC_SUCCESS);
    182: }

    184: static PetscErrorCode TaoTermComputeHessian_L1_Internal(TaoTerm term, Vec diag, Mat H)
    185: {
    186:   TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;

    188:   PetscFunctionBegin;
    189:   if (l1->epsilon == 0.0) {
    190:     PetscCall(TaoTermL1DerivativeCheck(term));
    191:     PetscCall(MatZeroEntries(H));
    192:     PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
    193:     PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
    194:   } else {
    195:     PetscCall(MatDiagonalSet(H, diag, INSERT_VALUES));
    196:   }
    197:   PetscFunctionReturn(PETSC_SUCCESS);
    198: }

    200: static PetscErrorCode TaoTermComputeHessian_L1(TaoTerm term, Vec x, Vec params, Mat H, Mat Hpre)
    201: {
    202:   Vec diag = NULL; /* Appease -Wmaybe-uninitialized */

    204:   PetscFunctionBegin;
    205:   if (H == NULL && Hpre == NULL) PetscFunctionReturn(PETSC_SUCCESS);
    206:   PetscCall(TaoTermL1ComputeDiag(term, x, params, &diag));
    207:   if (H) PetscCall(TaoTermComputeHessian_L1_Internal(term, diag, H));
    208:   if (Hpre && Hpre != H) PetscCall(TaoTermComputeHessian_L1_Internal(term, diag, Hpre));
    209:   PetscFunctionReturn(PETSC_SUCCESS);
    210: }

    212: static PetscErrorCode TaoTermCreateHessianMatrices_L1(TaoTerm term, Mat *H, Mat *Hpre)
    213: {
    214:   PetscBool is_hdiag, is_hprediag;

    216:   PetscFunctionBegin;
    217:   PetscCall(PetscInfo(term, "%s: Creating TAOTERML1 Hessian Matrices. TAOTERML1 only accepts MATDIAGONAL for MatType, overriding any user-set MatType.\n", ((PetscObject)term)->prefix));
    218:   PetscCall(PetscStrcmp(term->H_mattype, MATDIAGONAL, &is_hdiag));
    219:   PetscCall(PetscStrcmp(term->Hpre_mattype, MATDIAGONAL, &is_hprediag));
    220:   if (!is_hdiag) {
    221:     PetscCall(PetscFree(term->H_mattype));
    222:     PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->H_mattype));
    223:   }
    224:   if (!is_hprediag) {
    225:     PetscCall(PetscFree(term->Hpre_mattype));
    226:     PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->Hpre_mattype));
    227:   }
    228:   PetscCall(TaoTermCreateHessianMatricesDefault(term, H, Hpre));
    229:   PetscFunctionReturn(PETSC_SUCCESS);
    230: }

    232: /*@
    233:   TaoTermL1SetEpsilon - Set an $\epsilon$ smoothing parameter.

    235:   Logically collective

    237:   Input Parameters:
    238: + term    - a `TaoTerm` of type `TAOTERML1`
    239: - epsilon - a real number $\geq 0$

    241:   Options Database Keys:
    242: . -tao_term_l1_epsilon <real> - $\epsilon$

    244:   Level: advanced

    246:   If $\epsilon = 0$ (the default), then `term` computes $\|x - p\|_1$, but if $\epsilon > 0$, then it computes
    247:   $\sum_{i=0}^{n-1} \left(\sqrt{(x_i-p_i)^2 + \epsilon^2} - \epsilon\right)$.

    249: .seealso: [](sec_tao_term),
    250:           `TaoTerm`,
    251:           `TAOTERML1`,
    252:           `TaoTermL1GetEpsilon()`
    253: @*/
    254: PetscErrorCode TaoTermL1SetEpsilon(TaoTerm term, PetscReal epsilon)
    255: {
    256:   PetscFunctionBegin;
    259:   PetscTryMethod(term, "TaoTermL1SetEpsilon_C", (TaoTerm, PetscReal), (term, epsilon));
    260:   PetscFunctionReturn(PETSC_SUCCESS);
    261: }

    263: static PetscErrorCode TaoTermL1SetEpsilon_L1(TaoTerm term, PetscReal epsilon)
    264: {
    265:   TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;

    267:   PetscFunctionBegin;
    268:   PetscCheck(epsilon >= 0, PetscObjectComm((PetscObject)term), PETSC_ERR_ARG_OUTOFRANGE, "L1 epsilon (%g) cannot be < 0.0", (double)epsilon);
    269:   l1->epsilon = epsilon;
    270:   PetscFunctionReturn(PETSC_SUCCESS);
    271: }

    273: /*@
    274:   TaoTermL1GetEpsilon - Get the $\epsilon$ smoothing parameter set by `TaoTermL1SetEpsilon()`.

    276:   Not collective

    278:   Input Parameter:
    279: . term - a `TaoTerm` of type `TAOTERML1`

    281:   Output Parameter:
    282: . epsilon - the smoothing parameter

    284:   Level: advanced

    286: .seealso: [](sec_tao_term),
    287:           `TaoTerm`,
    288:           `TAOTERML1`,
    289:           `TaoTermL1SetEpsilon()`
    290: @*/
    291: PetscErrorCode TaoTermL1GetEpsilon(TaoTerm term, PetscReal *epsilon)
    292: {
    293:   PetscFunctionBegin;
    295:   PetscAssertPointer(epsilon, 2);
    296:   PetscUseMethod(term, "TaoTermL1GetEpsilon_C", (TaoTerm, PetscReal *), (term, epsilon));
    297:   PetscFunctionReturn(PETSC_SUCCESS);
    298: }

    300: static PetscErrorCode TaoTermL1GetEpsilon_L1(TaoTerm term, PetscReal *epsilon)
    301: {
    302:   TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;

    304:   PetscFunctionBegin;
    305:   *epsilon = l1->epsilon;
    306:   PetscFunctionReturn(PETSC_SUCCESS);
    307: }

    309: static PetscErrorCode TaoTermView_L1(TaoTerm term, PetscViewer viewer)
    310: {
    311:   TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
    312:   PetscBool   is_ascii;

    314:   PetscFunctionBegin;
    315:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &is_ascii));
    316:   if (is_ascii) PetscCall(PetscViewerASCIIPrintf(viewer, "epsilon (tao_term_l1_epsilon): %g\n", (double)l1->epsilon));
    317:   PetscFunctionReturn(PETSC_SUCCESS);
    318: }

    320: static PetscErrorCode TaoTermSetFromOptions_L1(TaoTerm term, PetscOptionItems PetscOptionsObject)
    321: {
    322:   TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
    323:   PetscBool   is_hdiag, is_hprediag;

    325:   PetscFunctionBegin;
    326:   PetscOptionsHeadBegin(PetscOptionsObject, "TaoTerm l1 options");
    327:   PetscCall(PetscOptionsBoundedReal("-tao_term_l1_epsilon", "smoothing parameter", "TaoTermL1SetEpsilon", l1->epsilon, &l1->epsilon, NULL, 0.0));
    328:   PetscOptionsHeadEnd();
    329:   PetscCall(PetscStrcmp(term->H_mattype, MATDIAGONAL, &is_hdiag));
    330:   PetscCall(PetscStrcmp(term->Hpre_mattype, MATDIAGONAL, &is_hprediag));
    331:   if (!is_hdiag) {
    332:     PetscCall(PetscFree(term->H_mattype));
    333:     PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->H_mattype));
    334:   }
    335:   if (!is_hprediag) {
    336:     PetscCall(PetscFree(term->Hpre_mattype));
    337:     PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->Hpre_mattype));
    338:   }
    339:   PetscFunctionReturn(PETSC_SUCCESS);
    340: }

    342: static PetscErrorCode TaoTermIsComputeHessianFDPossible_L1(TaoTerm term, PetscBool3 *ispossible)
    343: {
    344:   PetscFunctionBegin;
    345:   *ispossible = PETSC_BOOL3_FALSE;
    346:   PetscFunctionReturn(PETSC_SUCCESS);
    347: }

    349: /*MC
    350:   TAOTERML1 - A `TaoTerm` that computes $\|x - p\|_1$, for solution $x$ and parameters $p$.

    352:   Level: intermediate

    354:   Options Database Keys:
    355: . -tao_term_l1_epsilon <real> - (default 0.0) a smoothing parameter (see `TaoTermL1SetEpsilon()`)

    357:   Notes:
    358:   This term is `TAOTERM_PARAMETERS_OPTIONAL`.  If the parameters argument is `NULL` for
    359:   evaluation routines the term computes $\|x\|_1$.

    361:   This term has a smoothing parameter $\epsilon$ that defaults to 0: if $\epsilon > 0$,
    362:   the term computes a smooth approximation of $\|x - p\|_1$, see `TaoTermL1SetEpsilon()`.

    364:   The default Hessian creation mode (see `TaoTermGetCreateHessianMode()`) is `H == Hpre` and `TaoTermCreateHessianMatrices()`
    365:   will create a `MATDIAGONAL` for the Hessian.

    367: .seealso: [](sec_tao_term),
    368:           `TaoTerm`,
    369:           `TaoTermType`,
    370:           `TaoTermCreateL1()`,
    371:           `TaoTermL1GetEpsilon()`,
    372:           `TaoTermL1SetEpsilon()`,
    373:           `TAOTERMHALFL2SQUARED`,
    374:           `TAOTERMQUADRATIC`
    375: M*/
    376: PETSC_INTERN PetscErrorCode TaoTermCreate_L1(TaoTerm term)
    377: {
    378:   TaoTerm_L1 *l1;

    380:   PetscFunctionBegin;
    381:   PetscCall(TaoTermCreate_ElementwiseDivergence_Internal(term));
    382:   PetscCall(PetscNew(&l1));
    383:   term->data = (void *)l1;

    385:   PetscCall(PetscFree(term->H_mattype));
    386:   PetscCall(PetscFree(term->Hpre_mattype));

    388:   PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->H_mattype));
    389:   PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->Hpre_mattype));

    391:   term->ops->destroy                    = TaoTermDestroy_L1;
    392:   term->ops->view                       = TaoTermView_L1;
    393:   term->ops->setfromoptions             = TaoTermSetFromOptions_L1;
    394:   term->ops->objective                  = TaoTermComputeObjective_L1;
    395:   term->ops->gradient                   = TaoTermComputeGradient_L1;
    396:   term->ops->objectiveandgradient       = TaoTermComputeObjectiveAndGradient_L1;
    397:   term->ops->hessian                    = TaoTermComputeHessian_L1;
    398:   term->ops->createhessianmatrices      = TaoTermCreateHessianMatrices_L1;
    399:   term->ops->iscomputehessianfdpossible = TaoTermIsComputeHessianFDPossible_L1;

    401:   PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermL1SetEpsilon_C", TaoTermL1SetEpsilon_L1));
    402:   PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermL1GetEpsilon_C", TaoTermL1GetEpsilon_L1));
    403:   PetscFunctionReturn(PETSC_SUCCESS);
    404: }

    406: /*@
    407:   TaoTermCreateL1 - Create a `TaoTerm` for the objective function term $\|x - p\|_1$.

    409:   Collective

    411:   Input Parameters:
    412: + comm    - the MPI communicator where the term will be computed
    413: . n       - the local size of the $x$ and $p$ vectors (or `PETSC_DECIDE`)
    414: . N       - the global size of the $x$ and $p$ vectors (or `PETSC_DECIDE`)
    415: - epsilon - a non-negative smoothing parameter (see `TaoTermL1SetEpsilon()`)

    417:   Output Parameter:
    418: . term - the `TaoTerm`

    420:   Level: beginner

    422:   Note:
    423:   If you would like to add an L1 regularization term $\alpha \|x\|_1$ to the objective function of a `Tao`, do the following\:
    424: .vb
    425:   VecGetLocalSize(x, &n);
    426:   VecGetSize(x, &N);
    427:   TaoTermCreateL1(PetscObjectComm((PetscObject)x), n, N, 0.0, &term);
    428:   TaoAddTerm(tao, "reg_", alpha, term, NULL, NULL);
    429:   TaoTermDestroy(&term);
    430: .ve
    431:   If you would like to have a dictionary matrix term $\alpha \|D x\|_1$, do the same but pass `D` as the map of the term\:
    432: .vb
    433:   MatGetLocalSize(D, &m, NULL);
    434:   MatGetSize(D, &M, NULL);
    435:   TaoTermCreateL1(PetscObjectComm((PetscObject)D), m, M, 0.0, &term);
    436:   TaoAddTerm(tao, "reg_", alpha, term, NULL, D);
    437:   TaoTermDestroy(&term);
    438: .ve

    440: .seealso: [](sec_tao_term),
    441:           `TaoTerm`,
    442:           `TAOTERML1`,
    443:           `TaoTermL1GetEpsilon()`,
    444:           `TaoTermL1SetEpsilon()`,
    445:           `TaoTermCreateHalfL2Squared()`,
    446:           `TaoTermCreateQuadratic()`
    447: @*/
    448: PetscErrorCode TaoTermCreateL1(MPI_Comm comm, PetscInt n, PetscInt N, PetscReal epsilon, TaoTerm *term)
    449: {
    450:   TaoTerm _term;

    452:   PetscFunctionBegin;
    453:   PetscAssertPointer(term, 5);
    454:   PetscCheck(epsilon >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "L1 epsilon (%g) cannot be < 0.0", (double)epsilon);
    455:   PetscCall(TaoTermCreate(comm, &_term));
    456:   PetscCall(TaoTermSetType(_term, TAOTERML1));
    457:   PetscCall(TaoTermSetSolutionSizes(_term, n, N, 1));
    458:   PetscCall(TaoTermL1SetEpsilon(_term, epsilon));
    459:   *term = _term;
    460:   PetscFunctionReturn(PETSC_SUCCESS);
    461: }