Lines Matching defs:gn
7 TAO_BRGN *gn;
10 PetscCall(MatShellGetContext(H, &gn));
11 PetscCall(MatMult(gn->subsolver->ls_jac, in, gn->r_work));
12 PetscCall(MatMultTranspose(gn->subsolver->ls_jac, gn->r_work, out));
13 switch (gn->reg_type) {
15 PetscCall(MatMult(gn->Hreg, in, gn->x_work));
16 PetscCall(VecAXPY(out, gn->lambda, gn->x_work));
19 PetscCall(VecAXPY(out, gn->lambda, in));
22 PetscCall(VecAXPY(out, gn->lambda, in));
26 if (gn->D) {
27 PetscCall(MatMult(gn->D, in, gn->y)); /* y = D*in */
29 PetscCall(VecCopy(in, gn->y));
31 PetscCall(VecPointwiseMult(gn->y_work, gn->diag, gn->y)); /* y_work = diag.*(D*in), where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */
32 if (gn->D) {
33 PetscCall(MatMultTranspose(gn->D, gn->y_work, gn->x_work)); /* x_work = D'*(diag.*(D*in)) */
35 PetscCall(VecCopy(gn->y_work, gn->x_work));
37 PetscCall(VecAXPY(out, gn->lambda, gn->x_work));
40 PetscCall(VecPointwiseMult(gn->x_work, gn->damping, in));
41 PetscCall(VecAXPY(out, 1, gn->x_work));
46 static PetscErrorCode ComputeDamping(TAO_BRGN *gn)
54 PetscCall(VecGetArray(gn->damping, &damping_ary));
55 PetscCall(VecGetArrayRead(gn->diag, &diag_ary));
56 PetscCall(VecGetLocalSize(gn->damping, &n));
58 PetscCall(VecScale(gn->damping, gn->lambda));
59 PetscCall(VecRestoreArray(gn->damping, &damping_ary));
60 PetscCall(VecRestoreArrayRead(gn->diag, &diag_ary));
89 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
92 PetscCheck(gn->reg_type == TAOBRGN_REGULARIZATION_LM, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "Damping vector is only available if regularization type is lm.");
93 *d = gn->damping;
99 TAO_BRGN *gn = (TAO_BRGN *)ptr;
114 switch (gn->reg_type) {
116 PetscCall((*gn->regularizerobjandgrad)(tao, X, &f_reg, gn->x_work, gn->reg_obj_ctx));
117 *fcn += gn->lambda * f_reg;
118 PetscCall(VecAXPY(G, gn->lambda, gn->x_work));
123 *fcn += gn->lambda * 0.5 * f_reg;
125 PetscCall(VecAXPY(G, gn->lambda, X));
129 PetscCall(VecAXPBYPCZ(gn->x_work, 1.0, -1.0, 0.0, X, gn->x_old));
130 PetscCall(VecDot(gn->x_work, gn->x_work, &f_reg));
131 *fcn += gn->lambda * 0.5 * f_reg;
133 PetscCall(VecAXPBYPCZ(G, gn->lambda, -gn->lambda, 1.0, X, gn->x_old));
137 if (gn->D) {
138 PetscCall(MatMult(gn->D, X, gn->y)); /* y = D*x */
140 PetscCall(VecCopy(X, gn->y));
142 PetscCall(VecPointwiseMult(gn->y_work, gn->y, gn->y));
143 PetscCall(VecShift(gn->y_work, gn->epsilon * gn->epsilon));
144 PetscCall(VecSqrtAbs(gn->y_work)); /* gn->y_work = sqrt(y.^2+epsilon^2) */
145 PetscCall(VecSum(gn->y_work, &yESum));
146 PetscCall(VecGetSize(gn->y, &K));
147 *fcn += gn->lambda * (yESum - K * gn->epsilon);
149 PetscCall(VecPointwiseDivide(gn->y_work, gn->y, gn->y_work)); /* reuse y_work = y./sqrt(y.^2+epsilon^2) */
150 if (gn->D) {
151 PetscCall(MatMultTranspose(gn->D, gn->y_work, gn->x_work));
153 PetscCall(VecCopy(gn->y_work, gn->x_work));
155 PetscCall(VecAXPY(G, gn->lambda, gn->x_work));
167 TAO_BRGN *gn = (TAO_BRGN *)ptr;
173 if (gn->mat_explicit) PetscCall(MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_REUSE_MATRIX, PETSC_DETERMINE, &gn->H));
175 switch (gn->reg_type) {
177 PetscCall((*gn->regularizerhessian)(tao, X, gn->Hreg, gn->reg_hess_ctx));
178 if (gn->mat_explicit) PetscCall(MatAXPY(gn->H, 1.0, gn->Hreg, DIFFERENT_NONZERO_PATTERN));
181 if (gn->mat_explicit) PetscCall(MatShift(gn->H, gn->lambda));
184 if (gn->mat_explicit) PetscCall(MatShift(gn->H, gn->lambda));
188 if (gn->D) {
189 PetscCall(MatMult(gn->D, X, gn->y)); /* y = D*x */
191 PetscCall(VecCopy(X, gn->y));
193 PetscCall(VecPointwiseMult(gn->y_work, gn->y, gn->y));
194 PetscCall(VecShift(gn->y_work, gn->epsilon * gn->epsilon));
195 PetscCall(VecCopy(gn->y_work, gn->diag)); /* gn->diag = y.^2+epsilon^2 */
196 PetscCall(VecSqrtAbs(gn->y_work)); /* gn->y_work = sqrt(y.^2+epsilon^2) */
197 PetscCall(VecPointwiseMult(gn->diag, gn->y_work, gn->diag)); /* gn->diag = sqrt(y.^2+epsilon^2).^3 */
198 PetscCall(VecReciprocal(gn->diag));
199 PetscCall(VecScale(gn->diag, gn->epsilon * gn->epsilon));
200 if (gn->mat_explicit) PetscCall(MatDiagonalSet(gn->H, gn->diag, ADD_VALUES));
204 PetscCall(MatGetSize(gn->parent->ls_jac, NULL, &n));
206 PetscCall(MatGetColumnNorms(gn->parent->ls_jac, NORM_2, cnorms));
207 PetscCall(MatGetOwnershipRangeColumn(gn->parent->ls_jac, &cstart, &cend));
208 PetscCall(VecGetArray(gn->diag, &diag_ary));
210 PetscCall(VecRestoreArray(gn->diag, &diag_ary));
212 PetscCall(ComputeDamping(gn));
213 if (gn->mat_explicit) PetscCall(MatDiagonalSet(gn->H, gn->damping, ADD_VALUES));
223 TAO_BRGN *gn = (TAO_BRGN *)ctx;
227 gn->parent->nfuncs = tao->nfuncs;
228 gn->parent->ngrads = tao->ngrads;
229 gn->parent->nfuncgrads = tao->nfuncgrads;
230 gn->parent->nhess = tao->nhess;
231 gn->parent->niter = tao->niter;
232 gn->parent->ksp_its = tao->ksp_its;
233 gn->parent->ksp_tot_its = tao->ksp_tot_its;
234 gn->parent->fc = tao->fc;
235 PetscCall(TaoGetConvergedReason(tao, &gn->parent->reason));
238 PetscCall(VecSet(gn->x_old, 0.0));
240 PetscCall(VecCopy(tao->solution, gn->x_old));
241 PetscCall(VecCopy(tao->solution, gn->parent->solution));
244 PetscCall(VecCopy(tao->gradient, gn->parent->gradient));
247 if (gn->reg_type == TAOBRGN_REGULARIZATION_LM) {
249 if (gn->fc_old > tao->fc) {
250 gn->lambda = gn->lambda * gn->downhill_lambda_change;
253 gn->lambda = gn->lambda * gn->uphill_lambda_change;
256 gn->fc_old = tao->fc;
260 if (gn->parent->ops->update) PetscCall((*gn->parent->ops->update)(gn->parent, gn->parent->niter, gn->parent->user_update));
266 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
269 *type = gn->reg_type;
299 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
302 gn->reg_type = type;
330 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
333 PetscCall(TaoSolve(gn->subsolver));
335 tao->nfuncs = gn->subsolver->nfuncs;
336 tao->ngrads = gn->subsolver->ngrads;
337 tao->nfuncgrads = gn->subsolver->nfuncgrads;
338 tao->nhess = gn->subsolver->nhess;
339 tao->niter = gn->subsolver->niter;
340 tao->ksp_its = gn->subsolver->ksp_its;
341 tao->ksp_tot_its = gn->subsolver->ksp_tot_its;
342 PetscCall(TaoGetConvergedReason(gn->subsolver, &tao->reason));
344 PetscCall(VecCopy(gn->subsolver->solution, tao->solution));
345 PetscCall(VecCopy(gn->subsolver->gradient, tao->gradient));
351 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
356 PetscCall(PetscOptionsBool("-tao_brgn_mat_explicit", "switches the Hessian construction to be an explicit matrix rather than MATSHELL", "", gn->mat_explicit, &gn->mat_explicit, NULL));
357 PetscCall(PetscOptionsReal("-tao_brgn_regularizer_weight", "regularizer weight (default 1e-4)", "", gn->lambda, &gn->lambda, NULL));
358 PetscCall(PetscOptionsReal("-tao_brgn_l1_smooth_epsilon", "L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)", "", gn->epsilon, &gn->epsilon, NULL));
359 PetscCall(PetscOptionsReal("-tao_brgn_lm_downhill_lambda_change", "Factor to decrease trust region by on downhill steps", "", gn->downhill_lambda_change, &gn->downhill_lambda_change, NULL));
360 PetscCall(PetscOptionsReal("-tao_brgn_lm_uphill_lambda_change", "Factor to increase trust region by on uphill steps", "", gn->uphill_lambda_change, &gn->uphill_lambda_change, NULL));
361 PetscCall(PetscOptionsEnum("-tao_brgn_regularization_type", "regularization type", "", TaoBRGNRegularizationTypes, (PetscEnum)gn->reg_type, (PetscEnum *)&gn->reg_type, NULL));
364 if (gn->reg_type == TAOBRGN_REGULARIZATION_LM) {
365 PetscCall(TaoGetLineSearch(gn->subsolver, &ls));
368 PetscCall(TaoSetFromOptions(gn->subsolver));
374 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
381 PetscCall(PetscViewerASCIIPrintf(viewer, "Regularizer weight: %g\n", (double)gn->lambda));
382 PetscCall(PetscViewerASCIIPrintf(viewer, "BRGN Regularization Type: %s\n", TaoBRGNRegularizationTypes[gn->reg_type]));
383 switch (gn->reg_type) {
385 PetscCall(PetscViewerASCIIPrintf(viewer, "L1 smooth epsilon: %g\n", (double)gn->epsilon));
388 PetscCall(PetscViewerASCIIPrintf(viewer, "Downhill trust region decrease factor:: %g\n", (double)gn->downhill_lambda_change));
389 PetscCall(PetscViewerASCIIPrintf(viewer, "Uphill trust region increase factor:: %g\n", (double)gn->uphill_lambda_change));
400 PetscCall(TaoView(gn->subsolver, viewer));
407 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
413 PetscCall(PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNLS, &is_bnls));
414 PetscCall(PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTR, &is_bntr));
415 PetscCall(PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTL, &is_bntl));
418 if (!gn->x_work) PetscCall(VecDuplicate(tao->solution, &gn->x_work));
419 if (!gn->r_work) PetscCall(VecDuplicate(tao->ls_res, &gn->r_work));
420 if (!gn->x_old) {
421 PetscCall(VecDuplicate(tao->solution, &gn->x_old));
422 PetscCall(VecSet(gn->x_old, 0.0));
425 if (TAOBRGN_REGULARIZATION_L1DICT == gn->reg_type) {
426 if (!gn->y) {
427 if (gn->D) {
428 PetscCall(MatGetSize(gn->D, &K, &N)); /* Shell matrices still must have sizes defined. K = N for identity matrix, K=N-1 or N for gradient matrix */
429 PetscCall(MatCreateVecs(gn->D, NULL, &gn->y));
431 PetscCall(VecDuplicate(tao->solution, &gn->y)); /* If user does not setup dict matrix, use identity matrix, K=N */
433 PetscCall(VecSet(gn->y, 0.0));
435 if (!gn->y_work) PetscCall(VecDuplicate(gn->y, &gn->y_work));
436 if (!gn->diag) {
437 PetscCall(VecDuplicate(gn->y, &gn->diag));
438 PetscCall(VecSet(gn->diag, 0.0));
441 if (TAOBRGN_REGULARIZATION_LM == gn->reg_type) {
442 if (!gn->diag) PetscCall(MatCreateVecs(tao->ls_jac, &gn->diag, NULL));
443 if (!gn->damping) PetscCall(MatCreateVecs(tao->ls_jac, &gn->damping, NULL));
448 if (gn->mat_explicit) {
450 PetscCall(MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &gn->H));
454 PetscCall(MatCreate(PetscObjectComm((PetscObject)tao), &gn->H));
455 PetscCall(MatSetSizes(gn->H, n, n, N, N));
456 PetscCall(MatSetType(gn->H, MATSHELL));
457 PetscCall(MatSetOption(gn->H, MAT_SYMMETRIC, PETSC_TRUE));
458 PetscCall(MatShellSetOperation(gn->H, MATOP_MULT, (PetscErrorCodeFn *)GNHessianProd));
459 PetscCall(MatShellSetContext(gn->H, gn));
461 PetscCall(MatSetUp(gn->H));
463 PetscCall(TaoSetUpdate(gn->subsolver, GNHookFunction, gn));
464 PetscCall(TaoSetSolution(gn->subsolver, tao->solution));
465 if (tao->bounded) PetscCall(TaoSetVariableBounds(gn->subsolver, tao->XL, tao->XU));
466 PetscCall(TaoSetResidualRoutine(gn->subsolver, tao->ls_res, tao->ops->computeresidual, tao->user_lsresP));
467 PetscCall(TaoSetJacobianResidualRoutine(gn->subsolver, tao->ls_jac, tao->ls_jac, tao->ops->computeresidualjacobian, tao->user_lsjacP));
468 PetscCall(TaoSetObjectiveAndGradient(gn->subsolver, NULL, GNObjectiveGradientEval, gn));
469 PetscCall(TaoSetHessian(gn->subsolver, gn->H, gn->H, GNComputeHessian, gn));
471 PetscCall(TaoSetTolerances(gn->subsolver, tao->gatol, tao->grtol, tao->gttol));
472 PetscCall(TaoSetMaximumIterations(gn->subsolver, tao->max_it));
473 PetscCall(TaoSetMaximumFunctionEvaluations(gn->subsolver, tao->max_funcs));
474 PetscCall(TaoSetUp(gn->subsolver));
481 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
486 PetscCall(VecDestroy(&gn->x_work));
487 PetscCall(VecDestroy(&gn->r_work));
488 PetscCall(VecDestroy(&gn->x_old));
489 PetscCall(VecDestroy(&gn->diag));
490 PetscCall(VecDestroy(&gn->y));
491 PetscCall(VecDestroy(&gn->y_work));
493 PetscCall(VecDestroy(&gn->damping));
494 PetscCall(VecDestroy(&gn->diag));
495 PetscCall(MatDestroy(&gn->H));
496 PetscCall(MatDestroy(&gn->D));
497 PetscCall(MatDestroy(&gn->Hreg));
498 PetscCall(TaoDestroy(&gn->subsolver));
499 gn->parent = NULL;
536 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
539 *subsolver = gn->subsolver;
567 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
570 gn->lambda = lambda;
598 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
601 gn->epsilon = epsilon;
606 TaoBRGNSetDictionaryMatrix - bind the dictionary matrix from user application context to gn->D, for compressed sensing (with least-squares problem)
626 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
634 PetscCall(MatDestroy(&gn->D));
635 gn->D = dict;
669 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
672 if (ctx) gn->reg_obj_ctx = ctx;
673 if (func) gn->regularizerobjandgrad = func;
707 TAO_BRGN *gn = (TAO_BRGN *)tao->data;
714 if (ctx) gn->reg_hess_ctx = ctx;
715 if (func) gn->regularizerhessian = func;
718 PetscCall(MatDestroy(&gn->Hreg));
719 gn->Hreg = Hreg;
747 TAO_BRGN *gn;
750 PetscCall(PetscNew(&gn));
760 tao->data = gn;
761 gn->reg_type = TAOBRGN_REGULARIZATION_L2PROX;
762 gn->lambda = 1e-4;
763 gn->epsilon = 1e-6;
764 gn->downhill_lambda_change = 1. / 5.;
765 gn->uphill_lambda_change = 1.5;
766 gn->parent = tao;
768 PetscCall(TaoCreate(PetscObjectComm((PetscObject)tao), &gn->subsolver));
769 PetscCall(TaoSetType(gn->subsolver, TAOBNLS));
770 PetscCall(TaoSetOptionsPrefix(gn->subsolver, "tao_brgn_subsolver_"));