1 #include <../src/tao/leastsquares/impls/brgn/brgn.h> 2 3 static PetscErrorCode GNHessianProd(Mat H,Vec in,Vec out) 4 { 5 TAO_BRGN *gn; 6 PetscErrorCode ierr; 7 8 PetscFunctionBegin; 9 ierr = MatShellGetContext(H,&gn);CHKERRQ(ierr); 10 ierr = MatMult(gn->subsolver->ls_jac,in,gn->r_work);CHKERRQ(ierr); 11 ierr = MatMultTranspose(gn->subsolver->ls_jac,gn->r_work,out);CHKERRQ(ierr); 12 /* out = out + lambda*D'*(diag.*(D*in)) */ 13 ierr = MatMult(gn->D,in,gn->y);CHKERRQ(ierr); /* y = D*in */ 14 ierr = VecPointwiseMult(gn->y_work,gn->diag,gn->y);CHKERRQ(ierr); /* y_work = diag.*(D*in), where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */ 15 ierr = MatMultTranspose(gn->D,gn->y_work,gn->x_work);CHKERRQ(ierr); /* x_work = D'*(diag.*(D*in)) */ 16 ierr = VecAXPY(out,gn->lambda,gn->x_work);CHKERRQ(ierr); 17 18 PetscFunctionReturn(0); 19 } 20 21 static PetscErrorCode GNObjectiveGradientEval(Tao tao,Vec X,PetscReal *fcn,Vec G,void *ptr) 22 { 23 TAO_BRGN *gn = (TAO_BRGN *)ptr; 24 PetscInt K; /* dimension of D*X */ 25 PetscScalar yESum; 26 PetscErrorCode ierr; 27 28 PetscFunctionBegin; 29 /* compute objective *fcn*/ 30 /* compute first term ||ls_res||^2 */ 31 ierr = TaoComputeResidual(tao,X,tao->ls_res);CHKERRQ(ierr); 32 ierr = VecDotBegin(tao->ls_res,tao->ls_res,fcn);CHKERRQ(ierr); 33 ierr = VecDotEnd(tao->ls_res,tao->ls_res,fcn);CHKERRQ(ierr); 34 /* add the second term lambda*sum(sqrt(y.^2+epsilon^2) - epsilon), where y = D*x*/ 35 ierr = MatMult(gn->D,X,gn->y);CHKERRQ(ierr); /* y = D*x */ 36 ierr = VecPointwiseMult(gn->y_work,gn->y,gn->y);CHKERRQ(ierr); 37 ierr = VecShift(gn->y_work,gn->epsilon*gn->epsilon);CHKERRQ(ierr); 38 ierr = VecSqrtAbs(gn->y_work);CHKERRQ(ierr); /* gn->y_work = sqrt(y.^2+epsilon^2) */ 39 ierr = VecSum(gn->y_work,&yESum);CHKERRQ(ierr);CHKERRQ(ierr); 40 ierr = VecGetSize(gn->y,&K);CHKERRQ(ierr); 41 *fcn = 0.5*(*fcn) + gn->lambda*(yESum - K*gn->epsilon); 42 43 /* compute gradient G */ 44 ierr = TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);CHKERRQ(ierr); 45 ierr = MatMultTranspose(tao->ls_jac,tao->ls_res,G);CHKERRQ(ierr); 46 /* compute G = G + lambda*D'*(y./sqrt(y.^2+epsilon^2)),where y = D*x */ 47 ierr = VecPointwiseDivide(gn->y_work,gn->y,gn->y_work);CHKERRQ(ierr); /* reuse y_work = y./sqrt(y.^2+epsilon^2) */ 48 ierr = MatMultTranspose(gn->D,gn->y_work,gn->x_work);CHKERRQ(ierr); 49 ierr = VecAXPY(G,gn->lambda,gn->x_work);CHKERRQ(ierr); 50 51 PetscFunctionReturn(0); 52 } 53 54 55 static PetscErrorCode GNComputeHessian(Tao tao,Vec X,Mat H,Mat Hpre,void *ptr) 56 { 57 TAO_BRGN *gn = (TAO_BRGN *)ptr; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 ierr = TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);CHKERRQ(ierr); 62 63 /* calculate and store diagonal matrix as a vector: diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3* --> diag = epsilon^2 ./ sqrt(y.^2+epsilon^2).^3,where y = D*x */ 64 ierr = MatMult(gn->D,X,gn->y);CHKERRQ(ierr); /* y = D*x */ 65 ierr = VecPointwiseMult(gn->y_work,gn->y,gn->y);CHKERRQ(ierr); 66 ierr = VecShift(gn->y_work,gn->epsilon*gn->epsilon);CHKERRQ(ierr); 67 ierr = VecCopy(gn->y_work,gn->diag);CHKERRQ(ierr); /* gn->diag = y.^2+epsilon^2 */ 68 ierr = VecSqrtAbs(gn->y_work);CHKERRQ(ierr); /* gn->y_work = sqrt(y.^2+epsilon^2) */ 69 ierr = VecPointwiseMult(gn->diag,gn->y_work,gn->diag);CHKERRQ(ierr); /* gn->diag = sqrt(y.^2+epsilon^2).^3 */ 70 ierr = VecReciprocal(gn->diag);CHKERRQ(ierr); 71 ierr = VecScale(gn->diag,gn->epsilon*gn->epsilon);CHKERRQ(ierr); 72 73 PetscFunctionReturn(0); 74 } 75 76 static PetscErrorCode GNHookFunction(Tao tao,PetscInt iter) 77 { 78 TAO_BRGN *gn = (TAO_BRGN *)tao->user_update; 79 PetscErrorCode ierr; 80 81 PetscFunctionBegin; 82 /* Update basic tao information from the subsolver */ 83 gn->parent->nfuncs = tao->nfuncs; 84 gn->parent->ngrads = tao->ngrads; 85 gn->parent->nfuncgrads = tao->nfuncgrads; 86 gn->parent->nhess = tao->nhess; 87 gn->parent->niter = tao->niter; 88 gn->parent->ksp_its = tao->ksp_its; 89 gn->parent->ksp_tot_its = tao->ksp_tot_its; 90 ierr = TaoGetConvergedReason(tao,&gn->parent->reason);CHKERRQ(ierr); 91 /* Update the solution vectors */ 92 if (iter == 0) { 93 ierr = VecSet(gn->x_old,0.0);CHKERRQ(ierr); 94 } else { 95 ierr = VecCopy(tao->solution,gn->x_old);CHKERRQ(ierr); 96 ierr = VecCopy(tao->solution,gn->parent->solution);CHKERRQ(ierr); 97 } 98 /* Update the gradient */ 99 ierr = VecCopy(tao->gradient,gn->parent->gradient);CHKERRQ(ierr); 100 /* Call general purpose update function */ 101 if (gn->parent->ops->update) { 102 ierr = (*gn->parent->ops->update)(gn->parent,gn->parent->niter);CHKERRQ(ierr); 103 } 104 PetscFunctionReturn(0); 105 } 106 107 static PetscErrorCode TaoSolve_BRGN(Tao tao) 108 { 109 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 110 PetscErrorCode ierr; 111 112 PetscFunctionBegin; 113 ierr = TaoSolve(gn->subsolver);CHKERRQ(ierr); 114 /* Update basic tao information from the subsolver */ 115 tao->nfuncs = gn->subsolver->nfuncs; 116 tao->ngrads = gn->subsolver->ngrads; 117 tao->nfuncgrads = gn->subsolver->nfuncgrads; 118 tao->nhess = gn->subsolver->nhess; 119 tao->niter = gn->subsolver->niter; 120 tao->ksp_its = gn->subsolver->ksp_its; 121 tao->ksp_tot_its = gn->subsolver->ksp_tot_its; 122 ierr = TaoGetConvergedReason(gn->subsolver,&tao->reason);CHKERRQ(ierr); 123 /* Update vectors */ 124 ierr = VecCopy(gn->subsolver->solution,tao->solution);CHKERRQ(ierr); 125 ierr = VecCopy(gn->subsolver->gradient,tao->gradient);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 129 static PetscErrorCode TaoSetFromOptions_BRGN(PetscOptionItems *PetscOptionsObject,Tao tao) 130 { 131 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 ierr = PetscOptionsHead(PetscOptionsObject,"least-squares problems with L1 regularizer: ||f(x)||^2 + lambda*||x||_1. Currently L1-norm is approximated with smooth form");CHKERRQ(ierr); 136 ierr = PetscOptionsReal("-tao_brgn_lambda","L1-norm regularizer weight","",gn->lambda,&gn->lambda,NULL);CHKERRQ(ierr); 137 ierr = PetscOptionsReal("-tao_brgn_epsilon","L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon)","",gn->epsilon,&gn->epsilon,NULL);CHKERRQ(ierr); 138 ierr = PetscOptionsTail();CHKERRQ(ierr); 139 ierr = TaoSetFromOptions(gn->subsolver);CHKERRQ(ierr); 140 PetscFunctionReturn(0); 141 } 142 143 static PetscErrorCode TaoView_BRGN(Tao tao,PetscViewer viewer) 144 { 145 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 146 PetscErrorCode ierr; 147 148 PetscFunctionBegin; 149 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 150 ierr = TaoView(gn->subsolver,viewer);CHKERRQ(ierr); 151 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 152 PetscFunctionReturn(0); 153 } 154 155 static PetscErrorCode TaoSetUp_BRGN(Tao tao) 156 { 157 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 158 PetscErrorCode ierr; 159 PetscBool is_bnls,is_bntr,is_bntl; 160 PetscInt i,n,N,K; /* dict has size K*N*/ 161 /*PetscScalar v; */ /* XH: hack to set value of matrix */ 162 163 PetscFunctionBegin; 164 if (!tao->ls_res) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualRoutine() must be called before setup!"); 165 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNLS,&is_bnls);CHKERRQ(ierr); 166 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTR,&is_bntr);CHKERRQ(ierr); 167 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTL,&is_bntl);CHKERRQ(ierr); 168 if ((is_bnls || is_bntr || is_bntl) && !tao->ls_jac) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualJacobianRoutine() must be called before setup!"); 169 if (!tao->gradient){ 170 ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr); 171 } 172 if (!gn->x_work){ 173 ierr = VecDuplicate(tao->solution,&gn->x_work);CHKERRQ(ierr); 174 } 175 if (!gn->r_work){ 176 ierr = VecDuplicate(tao->ls_res,&gn->r_work);CHKERRQ(ierr); 177 } 178 if (!gn->x_old) { 179 ierr = VecDuplicate(tao->solution,&gn->x_old);CHKERRQ(ierr); 180 ierr = VecSet(gn->x_old,0.0);CHKERRQ(ierr); 181 } 182 183 /*ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr);*/ 184 /* TODO: Safeguard against NULL matrix */ 185 /*if (!gn->D)*/ 186 ierr = MatGetSize(gn->D,&K,&N);CHKERRQ(ierr); /* Shell matrices still must have sizes defined */ 187 /* K = N for identity matrix, K=N-1 or N for gradient matrix */ 188 if (!gn->y){ 189 ierr = VecCreate(PETSC_COMM_SELF,&gn->y);CHKERRQ(ierr); 190 ierr = VecSetSizes(gn->y,PETSC_DECIDE,K);CHKERRQ(ierr); 191 ierr = VecSetFromOptions(gn->y);CHKERRQ(ierr); 192 ierr = VecSet(gn->y,0.0);CHKERRQ(ierr); 193 194 } 195 if (!gn->y_work){ 196 ierr = VecDuplicate(gn->y,&gn->y_work);CHKERRQ(ierr); 197 } 198 if (!gn->diag){ 199 ierr = VecDuplicate(gn->y,&gn->diag);CHKERRQ(ierr); 200 ierr = VecSet(gn->diag,0.0);CHKERRQ(ierr); 201 } 202 203 /* XH: debug: check matrix */ 204 #if 0 205 ierr = PetscPrintf(PETSC_COMM_SELF,"-------- Check D matrix: -------- \n"); CHKERRQ(ierr); 206 ierr = MatView(gn->D,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 207 #endif 208 209 if (!tao->setupcalled) { 210 /* Hessian setup */ 211 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 212 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 213 ierr = MatSetSizes(gn->H,n,n,N,N);CHKERRQ(ierr); 214 ierr = MatSetType(gn->H,MATSHELL);CHKERRQ(ierr); 215 ierr = MatSetUp(gn->H);CHKERRQ(ierr); 216 ierr = MatShellSetOperation(gn->H,MATOP_MULT,(void (*)(void))GNHessianProd);CHKERRQ(ierr); 217 ierr = MatShellSetContext(gn->H,(void*)gn);CHKERRQ(ierr); 218 /* Subsolver setup,include initial vector and dicttionary D */ 219 ierr = TaoSetUpdate(gn->subsolver,GNHookFunction,(void*)gn);CHKERRQ(ierr); 220 ierr = TaoSetInitialVector(gn->subsolver,tao->solution);CHKERRQ(ierr); 221 if (tao->bounded) { 222 ierr = TaoSetVariableBounds(gn->subsolver,tao->XL,tao->XU);CHKERRQ(ierr); 223 } 224 ierr = TaoSetResidualRoutine(gn->subsolver,tao->ls_res,tao->ops->computeresidual,tao->user_lsresP);CHKERRQ(ierr); 225 ierr = TaoSetJacobianResidualRoutine(gn->subsolver,tao->ls_jac,tao->ls_jac,tao->ops->computeresidualjacobian,tao->user_lsjacP);CHKERRQ(ierr); 226 ierr = TaoSetObjectiveAndGradientRoutine(gn->subsolver,GNObjectiveGradientEval,(void*)gn);CHKERRQ(ierr); 227 ierr = TaoSetHessianRoutine(gn->subsolver,gn->H,gn->H,GNComputeHessian,(void*)gn);CHKERRQ(ierr); 228 /* Propagate some options down */ 229 ierr = TaoSetTolerances(gn->subsolver,tao->gatol,tao->grtol,tao->gttol);CHKERRQ(ierr); 230 ierr = TaoSetMaximumIterations(gn->subsolver,tao->max_it);CHKERRQ(ierr); 231 ierr = TaoSetMaximumFunctionEvaluations(gn->subsolver,tao->max_funcs);CHKERRQ(ierr); 232 for (i=0; i<tao->numbermonitors; ++i) { 233 ierr = TaoSetMonitor(gn->subsolver,tao->monitor[i],tao->monitorcontext[i],tao->monitordestroy[i]);CHKERRQ(ierr); 234 ierr = PetscObjectReference((PetscObject)(tao->monitorcontext[i]));CHKERRQ(ierr); 235 } 236 ierr = TaoSetUp(gn->subsolver);CHKERRQ(ierr); 237 } 238 PetscFunctionReturn(0); 239 } 240 241 static PetscErrorCode TaoDestroy_BRGN(Tao tao) 242 { 243 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 244 PetscErrorCode ierr; 245 246 PetscFunctionBegin; 247 if (tao->setupcalled) { 248 ierr = VecDestroy(&tao->gradient);CHKERRQ(ierr); 249 ierr = VecDestroy(&gn->x_work);CHKERRQ(ierr); 250 ierr = VecDestroy(&gn->r_work);CHKERRQ(ierr); 251 ierr = VecDestroy(&gn->x_old);CHKERRQ(ierr); 252 ierr = VecDestroy(&gn->diag);CHKERRQ(ierr); 253 ierr = VecDestroy(&gn->y);CHKERRQ(ierr); 254 ierr = VecDestroy(&gn->y_work);CHKERRQ(ierr); 255 } 256 ierr = MatDestroy(&gn->H);CHKERRQ(ierr); 257 ierr = MatDestroy(&gn->D);CHKERRQ(ierr); 258 ierr = TaoDestroy(&gn->subsolver);CHKERRQ(ierr); 259 gn->parent = NULL; 260 ierr = PetscFree(tao->data);CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 264 /*MC 265 TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares 266 problems with bound constraints. This algorithm is a thin wrapper around TAOBNTL 267 that constructs the Guass-Newton problem with the user-provided least-squares 268 residual and Jacobian. The problem is regularized with an L2-norm proximal point 269 term. 270 271 Options Database Keys: 272 + -tao_bqnk_max_cg_its - maximum number of bounded conjugate-gradient iterations taken in each Newton loop 273 . -tao_bqnk_init_type - trust radius initialization method ("constant", "direction", "interpolation") 274 . -tao_bqnk_update_type - trust radius update method ("step", "direction", "interpolation") 275 - -tao_bqnk_as_type - active-set estimation method ("none", "bertsekas") 276 277 Level: beginner 278 M*/ 279 PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao) 280 { 281 TAO_BRGN *gn; 282 PetscErrorCode ierr; 283 284 PetscFunctionBegin; 285 ierr = PetscNewLog(tao,&gn);CHKERRQ(ierr); 286 287 tao->ops->destroy = TaoDestroy_BRGN; 288 tao->ops->setup = TaoSetUp_BRGN; 289 tao->ops->setfromoptions = TaoSetFromOptions_BRGN; 290 tao->ops->view = TaoView_BRGN; 291 tao->ops->solve = TaoSolve_BRGN; 292 293 tao->data = (void*)gn; 294 gn->lambda = 1e-4; 295 gn->epsilon = 1e-6; 296 gn->parent = tao; 297 298 ierr = MatCreate(PetscObjectComm((PetscObject)tao),&gn->H);CHKERRQ(ierr); 299 ierr = MatSetOptionsPrefix(gn->H,"tao_brgn_hessian_");CHKERRQ(ierr); 300 301 ierr = TaoCreate(PetscObjectComm((PetscObject)tao),&gn->subsolver);CHKERRQ(ierr); 302 ierr = TaoSetType(gn->subsolver,TAOBNLS);CHKERRQ(ierr); 303 ierr = TaoSetOptionsPrefix(gn->subsolver,"tao_brgn_subsolver_");CHKERRQ(ierr); 304 PetscFunctionReturn(0); 305 } 306 307 /*@C 308 TaoBRGNGetSubsolver - Get the pointer to the subsolver inside BRGN 309 310 Collective on Tao 311 312 Level: developer 313 314 Input Parameters: 315 + tao - the Tao solver context 316 - subsolver - the Tao sub-solver context 317 @*/ 318 PetscErrorCode TaoBRGNGetSubsolver(Tao tao,Tao *subsolver) 319 { 320 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 321 322 PetscFunctionBegin; 323 *subsolver = gn->subsolver; 324 PetscFunctionReturn(0); 325 } 326 327 /*@C 328 TaoBRGNSetL1RegularizerWeight - Set the L1-norm regularizer weight for the Gauss-Newton least-squares algorithm 329 330 Collective on Tao 331 332 Level: developer 333 334 Input Parameters: 335 + tao - the Tao solver context 336 - lambda - L1-norm regularizer weight 337 @*/ 338 PetscErrorCode TaoBRGNSetL1RegularizerWeight(Tao tao,PetscReal lambda) 339 { 340 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 341 342 /* Initialize lambda here */ 343 344 PetscFunctionBegin; 345 gn->lambda = lambda; 346 PetscFunctionReturn(0); 347 } 348 349 /*@C 350 TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm 351 352 Collective on Tao 353 354 Level: developer 355 356 Input Parameters: 357 + tao - the Tao solver context 358 - epsilon - L1-norm smooth approximation parameter 359 @*/ 360 PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao,PetscReal epsilon) 361 { 362 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 363 364 /* Initialize epsilon here */ 365 366 PetscFunctionBegin; 367 gn->epsilon = epsilon; 368 PetscFunctionReturn(0); 369 } 370 371 /*@C 372 TaoBRGNSetDictionaryMatrix - bind the dictionary matrix from user application context to gn->D, for compressed sensing (with least-squares problem) 373 374 Input Parameters: 375 + tao - the Tao context 376 . dict - the user specified dictionary matrix 377 378 Level: developer 379 @*/ 380 PetscErrorCode TaoBRGNSetDictionaryMatrix(Tao tao,Mat dict) 381 { 382 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 383 PetscErrorCode ierr; 384 PetscFunctionBegin; 385 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 386 if (dict) { 387 PetscValidHeaderSpecific(dict,MAT_CLASSID,2); 388 /*PetscCheckSameComm(tao,1,dict,2);*/ 389 ierr = PetscObjectReference((PetscObject)dict);CHKERRQ(ierr); 390 } 391 ierr = MatDestroy(&gn->D);CHKERRQ(ierr); 392 gn->D = dict; /* We allow to set a null dictionary, which means we just use default identity matrix? */ 393 PetscFunctionReturn(0); 394 } 395 396 /* XH: 397 Changed TaoBRGNSetTikhonovLambda --> TaoBRGNSetL1RegularizerWeight in brgn.c, peststao.h, and zbrgnf.c. 398 Added TaoBRGNSetL1SmoothEpsilon by following TaoBRGNSetL1RegularizerWeight. 399 Added TaoBRGNSetDictionaryMatrix by following TaoBRGNSetL1RegularizerWeight 400 Maybe change D*x to D(x), and A*x to A(x) as function handle 401 Maybe need to also keep y = D*x, to avoid duplicate frequent computation of D*x 402 */