Actual source code: taotermhalfl2squared.c

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

      3: typedef struct _n_TaoTerm_HalfL2Squared TaoTerm_HalfL2Squared;

      5: struct _n_TaoTerm_HalfL2Squared {
      6:   Vec pdiff_work;
      7: };

      9: static PetscErrorCode TaoTermDestroy_Halfl2squared(TaoTerm term)
     10: {
     11:   TaoTerm_HalfL2Squared *l2 = (TaoTerm_HalfL2Squared *)term->data;

     13:   PetscFunctionBegin;
     14:   PetscCall(VecDestroy(&l2->pdiff_work));
     15:   PetscCall(PetscFree(l2));
     16:   term->data = NULL;
     17:   PetscCall(TaoTermDestroy_ElementwiseDivergence_Internal(term));
     18:   PetscFunctionReturn(PETSC_SUCCESS);
     19: }

     21: static PetscErrorCode TaoTermComputeObjective_Halfl2squared(TaoTerm term, Vec x, Vec params, PetscReal *value)
     22: {
     23:   PetscScalar            v;
     24:   TaoTerm_HalfL2Squared *l2 = (TaoTerm_HalfL2Squared *)term->data;

     26:   PetscFunctionBegin;
     27:   if (params) {
     28:     if (l2->pdiff_work == NULL) PetscCall(VecDuplicate(x, &l2->pdiff_work));

     30:     PetscCall(VecCopy(x, l2->pdiff_work));
     31:     PetscCall(VecAXPY(l2->pdiff_work, -1.0, params));
     32:     PetscCall(VecDot(l2->pdiff_work, l2->pdiff_work, &v));
     33:   } else PetscCall(VecDot(x, x, &v));
     34:   *value = 0.5 * PetscRealPart(v);
     35:   PetscFunctionReturn(PETSC_SUCCESS);
     36: }

     38: static PetscErrorCode TaoTermComputeObjectiveAndGradient_Halfl2squared(TaoTerm term, Vec x, Vec params, PetscReal *value, Vec g)
     39: {
     40:   PetscScalar v;

     42:   PetscFunctionBegin;
     43:   if (params) PetscCall(VecWAXPY(g, -1.0, params, x));
     44:   else PetscCall(VecCopy(x, g));
     45:   PetscCall(VecDot(g, g, &v));
     46:   *value = 0.5 * PetscRealPart(v);
     47:   PetscFunctionReturn(PETSC_SUCCESS);
     48: }

     50: static PetscErrorCode TaoTermComputeGradient_Halfl2squared(TaoTerm term, Vec x, Vec params, Vec g)
     51: {
     52:   PetscFunctionBegin;
     53:   if (params) PetscCall(VecWAXPY(g, -1.0, params, x));
     54:   else PetscCall(VecCopy(x, g));
     55:   PetscFunctionReturn(PETSC_SUCCESS);
     56: }

     58: static PetscErrorCode TaoTermComputeHessian_Halfl2squared(TaoTerm term, Vec x, Vec params, Mat H, Mat Hpre)
     59: {
     60:   PetscFunctionBegin;
     61:   if (H) {
     62:     PetscCall(MatZeroEntries(H));
     63:     PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
     64:     PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
     65:     PetscCall(MatShift(H, 1.0));
     66:   }
     67:   if (Hpre && Hpre != H) {
     68:     PetscCall(MatZeroEntries(Hpre));
     69:     PetscCall(MatAssemblyBegin(Hpre, MAT_FINAL_ASSEMBLY));
     70:     PetscCall(MatAssemblyEnd(Hpre, MAT_FINAL_ASSEMBLY));
     71:     PetscCall(MatShift(Hpre, 1.0));
     72:   }
     73:   PetscFunctionReturn(PETSC_SUCCESS);
     74: }

     76: static PetscErrorCode TaoTermCreateHessianMatrices_Halfl2squared(TaoTerm term, Mat *H, Mat *Hpre)
     77: {
     78:   PetscBool is_hdiag, is_hprediag;

     80:   PetscFunctionBegin;
     81:   PetscCall(PetscInfo(term, "%s: Creating TAOTERMHALFL2SQUARED Hessian Matrices. TAOTERMHALFL2SQUARED only accepts MATDIAGONAL for MatType, overriding any user-set MatType.\n", ((PetscObject)term)->prefix));
     82:   PetscCall(PetscStrcmp(term->H_mattype, MATDIAGONAL, &is_hdiag));
     83:   PetscCall(PetscStrcmp(term->Hpre_mattype, MATDIAGONAL, &is_hprediag));
     84:   if (!is_hdiag) {
     85:     PetscCall(PetscFree(term->H_mattype));
     86:     PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->H_mattype));
     87:   }
     88:   if (!is_hprediag) {
     89:     PetscCall(PetscFree(term->Hpre_mattype));
     90:     PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->Hpre_mattype));
     91:   }
     92:   PetscCall(TaoTermCreateHessianMatricesDefault(term, H, Hpre));
     93:   PetscFunctionReturn(PETSC_SUCCESS);
     94: }

     96: static PetscErrorCode TaoTermIsComputeHessianFDPossible_Halfl2squared(TaoTerm term, PetscBool3 *ispossible)
     97: {
     98:   PetscFunctionBegin;
     99:   *ispossible = PETSC_BOOL3_FALSE;
    100:   PetscFunctionReturn(PETSC_SUCCESS);
    101: }

    103: /*MC
    104:   TAOTERMHALFL2SQUARED - A `TaoTerm` that computes $\tfrac{1}{2}\|x - p\|_2^2$, for solution $x$ and parameters $p$.

    106:   Level: intermediate

    108:   Notes:
    109:   By default this term is `TAOTERM_PARAMETERS_OPTIONAL`.  If the parameters
    110:   argument is `NULL` in the evaluation routines (`TaoTermComputeObjective()`,
    111:   `TaoTermComputeGradient()`, etc.), then it is assumed $p = 0$ and the term computes
    112:   $\tfrac{1}{2}\|x\|_2^2$.

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

    117: .seealso: [](sec_tao_term),
    118:           `TaoTerm`,
    119:           `TaoTermType`,
    120:           `TaoTermCreateHalfL2Squared()`,
    121:           `TAOTERML1`,
    122:           `TAOTERMQUADRATIC`
    123: M*/
    124: PETSC_INTERN PetscErrorCode TaoTermCreate_Halfl2squared(TaoTerm term)
    125: {
    126:   TaoTerm_HalfL2Squared *l2;

    128:   PetscFunctionBegin;
    129:   PetscCall(TaoTermCreate_ElementwiseDivergence_Internal(term));
    130:   PetscCall(PetscNew(&l2));
    131:   term->data = (void *)l2;

    133:   PetscCall(PetscFree(term->H_mattype));
    134:   PetscCall(PetscFree(term->Hpre_mattype));

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

    139:   term->ops->destroy                    = TaoTermDestroy_Halfl2squared;
    140:   term->ops->objective                  = TaoTermComputeObjective_Halfl2squared;
    141:   term->ops->gradient                   = TaoTermComputeGradient_Halfl2squared;
    142:   term->ops->objectiveandgradient       = TaoTermComputeObjectiveAndGradient_Halfl2squared;
    143:   term->ops->hessian                    = TaoTermComputeHessian_Halfl2squared;
    144:   term->ops->createhessianmatrices      = TaoTermCreateHessianMatrices_Halfl2squared;
    145:   term->ops->iscomputehessianfdpossible = TaoTermIsComputeHessianFDPossible_Halfl2squared;
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

    149: /*@
    150:   TaoTermCreateHalfL2Squared - Create a `TaoTerm` for the objective term $\tfrac{1}{2}\|x - p\|_2^2$, for solution $x$ and parameters $p$.

    152:   Collective

    154:   Input Parameters:
    155: + comm - the MPI communicator where the `TaoTerm` will be computed
    156: . n    - the local size of the $x$ and $p$ vectors (or `PETSC_DECIDE`)
    157: - N    - the global size of the $x$ and $p$ vectors (or `PETSC_DECIDE`)

    159:   Output Parameter:
    160: . term - the `TaoTerm`

    162:   Level: beginner

    164:   Note:
    165:   If you would like to add a Tikhonov regularization term $\alpha \tfrac{1}{2}\|x\|_2^2$ to the objective function of a `Tao`, do the following\:
    166: .vb
    167:   VecGetSizes(x, &n, &N);
    168:   TaoTermCreateHalfL2Squared(PetscObjectComm((PetscObject)x), n, N, &term);
    169:   TaoAddTerm(tao, "reg_", alpha, term, NULL, NULL);
    170:   TaoTermDestroy(&term);
    171: .ve
    172:   If you would like to add a biased regularization term $\alpha \tfrac{1}{2}\|x - p \|_2^2$, do the same but pass `p` as the parameters of the term\:
    173: .vb
    174:   TaoAddTerm(tao, "reg_", alpha, term, p, NULL);
    175: .ve

    177: .seealso: [](sec_tao_term),
    178:           `TaoTerm`,
    179:           `TAOTERMHALFL2SQUARED`,
    180:           `TaoTermCreateL1()`,
    181:           `TaoTermCreateQuadratic()`
    182: @*/
    183: PetscErrorCode TaoTermCreateHalfL2Squared(MPI_Comm comm, PetscInt n, PetscInt N, TaoTerm *term)
    184: {
    185:   TaoTerm _term;

    187:   PetscFunctionBegin;
    188:   PetscAssertPointer(term, 4);
    189:   PetscCall(TaoTermCreate(comm, &_term));
    190:   PetscCall(TaoTermSetType(_term, TAOTERMHALFL2SQUARED));
    191:   PetscCall(TaoTermSetSolutionSizes(_term, n, N, 1));
    192:   *term = _term;
    193:   PetscFunctionReturn(PETSC_SUCCESS);
    194: }