1 /* XH: todo add cs1f.F90 and asjust makefile */ 2 /* 3 Include "petsctao.h" so that we can use TAO solvers. Note that this 4 file automatically includes libraries such as: 5 petsc.h - base PETSc routines petscvec.h - vectors 6 petscsys.h - sysem routines petscmat.h - matrices 7 petscis.h - index sets petscksp.h - Krylov subspace methods 8 petscviewer.h - viewers petscpc.h - preconditioners 9 10 */ 11 12 #include <petsctao.h> 13 14 /* 15 Description: Compressive sensing test example 1. 16 0.5*||Ax-b||^2 + lambda*||D*x||_1 17 Xiang Huang: Nov 19, 2018 18 19 Reference: None 20 */ 21 22 static char help[] = "Finds the least-squares solution to the under constraint linear model Ax = b, with L1-norm regularizer. \n\ 23 A is a M*N real matrix (M<N), x is sparse. \n\ 24 We find the sparse solution by solving 0.5*||Ax-b||^2 + lambda*||D*x||_1, where lambda (by default 1e-4) is a user specified weight.\n\ 25 D is the K*N transform matrix so that D*x is sparse. By default D is identity matrix, so that D*x = x.\n"; 26 /*T 27 Concepts: TAO^Solving a system of nonlinear equations, nonlinear least squares 28 Routines: TaoCreate(); 29 Routines: TaoSetType(); 30 Routines: TaoSetSeparableObjectiveRoutine(); 31 Routines: TaoSetJacobianRoutine(); 32 Routines: TaoSetInitialVector(); 33 Routines: TaoSetFromOptions(); 34 Routines: TaoSetConvergenceHistory(); TaoGetConvergenceHistory(); 35 Routines: TaoSolve(); 36 Routines: TaoView(); TaoDestroy(); 37 Processors: 1 38 T*/ 39 40 #define M 3 41 #define N 5 42 #define K 4 43 44 /* User-defined application context */ 45 typedef struct { 46 /* Working space. linear least square: f(x) = A*x - b */ 47 PetscReal A[M][N]; /* array of coefficients */ 48 PetscReal b[M]; /* array of observations */ 49 PetscReal xGT[M]; /* array of ground truth object, which can be used to compare the reconstruction result */ 50 PetscReal D[K][N]; /* array of coefficients for 0.5*||Ax-b||^2 + lambda*||D*x||_1 */ 51 PetscReal J[M][N]; /* dense jacobian matrix array. For linear least square, J = A. For nonlinear least square, it is different from A */ 52 PetscInt idm[M]; /* Matrix row, column indices for jacobian and dictionary */ 53 PetscInt idn[N]; 54 PetscInt idk[K]; 55 } AppCtx; 56 57 /* User provided Routines */ 58 PetscErrorCode InitializeUserData(AppCtx *); 59 PetscErrorCode FormStartingPoint(Vec); 60 PetscErrorCode FormDictionaryMatrix(Mat,AppCtx *); 61 PetscErrorCode EvaluateFunction(Tao,Vec,Vec,void *); 62 PetscErrorCode EvaluateJacobian(Tao,Vec,Mat,Mat,void *); 63 64 /*--------------------------------------------------------------------*/ 65 int main(int argc,char **argv) 66 { 67 PetscErrorCode ierr; /* used to check for functions returning nonzeros */ 68 Vec x,f; /* solution, function f(x) = A*x-b */ 69 Mat J,D; /* Jacobian matrix, Transform matrix */ 70 Tao tao; /* Tao solver context */ 71 PetscInt i; /* iteration information */ 72 PetscReal hist[100],resid[100]; 73 PetscInt lits[100]; 74 AppCtx user; /* user-defined work context */ 75 76 ierr = PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr; 77 78 /* Allocate solution and vector function vectors */ 79 ierr = VecCreateSeq(PETSC_COMM_SELF,N,&x);CHKERRQ(ierr); 80 ierr = VecCreateSeq(PETSC_COMM_SELF,M,&f);CHKERRQ(ierr); 81 82 /* Allocate Jacobian and Dictionary matrix. */ 83 ierr = MatCreateSeqDense(PETSC_COMM_SELF,M,N,NULL,&J);CHKERRQ(ierr); 84 ierr = MatCreateSeqDense(PETSC_COMM_SELF,K,N,NULL,&D);CHKERRQ(ierr); /* XH: TODO: dense -> sparse/dense/shell etc, do it on fly */ 85 86 for (i=0;i<M;i++) user.idm[i] = i; 87 for (i=0;i<N;i++) user.idn[i] = i; 88 for (i=0;i<K;i++) user.idk[i] = i; 89 90 /* Create TAO solver and set desired solution method */ 91 ierr = TaoCreate(PETSC_COMM_SELF,&tao);CHKERRQ(ierr); 92 ierr = TaoSetType(tao,TAOBRGN);CHKERRQ(ierr); 93 94 /* User set application context: A, D matrice, and b vector. */ 95 ierr = InitializeUserData(&user);CHKERRQ(ierr); 96 97 /* Set initial guess */ 98 ierr = FormStartingPoint(x);CHKERRQ(ierr); 99 100 /* Fill the content of matrix D from user application Context */ 101 ierr = FormDictionaryMatrix(D,&user);CHKERRQ(ierr); 102 103 /* Bind x to tao->solution. */ 104 ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr); 105 /* Bind D to tao->data->D */ 106 ierr = TaoBRGNSetDictionaryMatrix(tao,D);CHKERRQ(ierr); 107 108 /* Set the function and Jacobian routines. */ 109 ierr = TaoSetResidualRoutine(tao,f,EvaluateFunction,(void*)&user);CHKERRQ(ierr); 110 ierr = TaoSetJacobianResidualRoutine(tao,J,J,EvaluateJacobian,(void*)&user);CHKERRQ(ierr); 111 112 /* Check for any TAO command line arguments */ 113 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 114 115 ierr = TaoSetConvergenceHistory(tao,hist,resid,0,lits,100,PETSC_TRUE);CHKERRQ(ierr); 116 117 /* Perform the Solve */ 118 ierr = TaoSolve(tao);CHKERRQ(ierr); 119 120 /* XH: Debug: View the result, function and Jacobian. */ 121 ierr = PetscPrintf(PETSC_COMM_SELF, "-------- result x, residual f=A*x-b, and Jacobian=A. -------- \n");CHKERRQ(ierr); 122 ierr = VecView(x,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 123 ierr = VecView(f,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 124 ierr = MatView(J,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 125 ierr = MatView(D,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 126 127 /* Free TAO data structures */ 128 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 129 130 /* Free PETSc data structures */ 131 ierr = VecDestroy(&x);CHKERRQ(ierr); 132 ierr = VecDestroy(&f);CHKERRQ(ierr); 133 ierr = MatDestroy(&J);CHKERRQ(ierr); 134 ierr = MatDestroy(&D);CHKERRQ(ierr); 135 136 ierr = PetscFinalize(); 137 return ierr; 138 } 139 140 /*--------------------------------------------------------------------*/ 141 PetscErrorCode EvaluateFunction(Tao tao, Vec X, Vec F, void *ptr) 142 { 143 AppCtx *user = (AppCtx *)ptr; 144 PetscInt m,n; 145 const PetscReal *x; 146 PetscReal *b=user->b,*f; 147 PetscErrorCode ierr; 148 149 PetscFunctionBegin; 150 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 151 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 152 153 /* Even for linear least square, we do not direct use matrix operation f = A*x - b now, just for future modification and compatability for nonlinear least square */ 154 for (m=0;m<M;m++) { 155 f[m] = -b[m]; 156 for (n=0;n<N;n++) { 157 f[m] += user->A[m][n]*x[n]; 158 } 159 } 160 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 161 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 162 PetscLogFlops(M*N*2); 163 PetscFunctionReturn(0); 164 } 165 166 /*------------------------------------------------------------*/ 167 /* J[m][n] = df[m]/dx[n] */ 168 PetscErrorCode EvaluateJacobian(Tao tao, Vec X, Mat J, Mat Jpre, void *ptr) 169 { 170 AppCtx *user = (AppCtx *)ptr; 171 PetscInt m,n; 172 const PetscReal *x; 173 PetscErrorCode ierr; 174 175 PetscFunctionBegin; 176 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); /* not used for linear least square, but keep for future nonlinear least square) */ 177 /* XH: TODO: For linear least square, we can just set J=A fixed once, instead of keep update it! Maybe just create a function getFixedJacobian? 178 For nonlinear least square, we require x to compute J, keep codes here for future nonlinear least square*/ 179 for (m=0; m<M; ++m) { 180 for (n=0; n<N; ++n) { 181 user->J[m][n] = user->A[m][n]; 182 } 183 } 184 185 ierr = MatSetValues(J,M,user->idm,N,user->idn,(PetscReal *)user->J,INSERT_VALUES);CHKERRQ(ierr); 186 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 187 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 188 189 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);/* not used for linear least square, but keep for future nonlinear least square) */ 190 PetscLogFlops(0); /* 0 for linear least square, >0 for nonlinear least square */ 191 PetscFunctionReturn(0); 192 } 193 194 /* ------------------------------------------------------------ */ 195 /* Currently fixed matrix, in future may be dynamic for D(x)? */ 196 PetscErrorCode FormDictionaryMatrix(Mat D,AppCtx *user) 197 { 198 PetscErrorCode ierr; 199 200 PetscFunctionBegin; 201 ierr = MatSetValues(D,K,user->idk,N,user->idn,(PetscReal *)user->D,INSERT_VALUES);CHKERRQ(ierr); 202 ierr = MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 203 ierr = MatAssemblyEnd(D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 204 205 PetscLogFlops(0); /* 0 for fixed dictionary matrix, >0 for varying dictionary matrix */ 206 PetscFunctionReturn(0); 207 } 208 209 /* ------------------------------------------------------------ */ 210 PetscErrorCode FormStartingPoint(Vec X) 211 { 212 PetscErrorCode ierr; 213 PetscFunctionBegin; 214 ierr = VecSet(X,0.0);CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 218 /* ---------------------------------------------------------------------- */ 219 PetscErrorCode InitializeUserData(AppCtx *user) 220 { 221 PetscReal *b=user->b; /* **A=user->A, but we don't kown the dimension of A in this way, how to fix? */ 222 PetscInt m,n,k; /* loop index for M,N,K dimension. */ 223 224 PetscFunctionBegin; 225 /* b = A*x while x = [0;0;1;0;0] here*/ 226 m = 0; 227 b[m++] = 0.28; 228 b[m++] = 0.55; 229 b[m++] = 0.96; 230 231 /* matlab generated random matrix, uniformly distributed in [0,1] with 2 digits accuracy. rng(0); A = rand(M, N); A = round(A*100)/100; 232 A = [0.81 0.91 0.28 0.96 0.96 233 0.91 0.63 0.55 0.16 0.49 234 0.13 0.10 0.96 0.97 0.80] 235 */ 236 m=0; n=0; user->A[m][n++] = 0.81; user->A[m][n++] = 0.91; user->A[m][n++] = 0.28; user->A[m][n++] = 0.96; user->A[m][n++] = 0.96; 237 ++m; n=0; user->A[m][n++] = 0.91; user->A[m][n++] = 0.63; user->A[m][n++] = 0.55; user->A[m][n++] = 0.16; user->A[m][n++] = 0.49; 238 ++m; n=0; user->A[m][n++] = 0.13; user->A[m][n++] = 0.10; user->A[m][n++] = 0.96; user->A[m][n++] = 0.97; user->A[m][n++] = 0.80; 239 240 /* initialize to 0 */ 241 for (k=0; k<K; k++) { 242 for (n=0; n<N; n++) { 243 user->D[k][n] = 0.0; 244 } 245 } 246 /* Choice I: set D to identity matrix of size N*N for testing */ 247 /* for (k=0; k<K; k++) user->D[k][k] = 1.0; */ 248 /* Choice II: set D to Backward difference matrix of size (N-1)*N, with zero extended boundary assumption */ 249 for (k=0;k<K;k++) { 250 user->D[k][k] = -1.0; 251 user->D[k][k+1] = 1.0; 252 } 253 254 PetscFunctionReturn(0); 255 } 256 257 /*TEST 258 259 build: 260 requires: !complex !single !quad !define(PETSC_USE_64BIT_INDICES) 261 262 test: 263 localrunfiles: cs1Data_A_b_xGT 264 args: -tao_smonitor -tao_max_it 100 -tao_type pounders -tao_gatol 1.e-6 265 266 test: 267 suffix: 2 268 localrunfiles: cs1Data_A_b_xGT 269 args: -tao_monitor -tao_max_it 100 -tao_type brgn -tao_brgn_regularization_type l2prox -tao_brgn_regularizer_weight 1e-8 -tao_gatol 1.e-6 270 271 test: 272 suffix: 3 273 localrunfiles: cs1Data_A_b_xGT 274 args: -tao_monitor -tao_max_it 100 -tao_type brgn -tao_brgn_regularization_type l1dict -tao_brgn_regularizer_weight 1e-8 -tao_brgn_l1_smooth_epsilon 1e-6 -tao_gatol 1.e-6 275 276 test: 277 suffix: 4 278 localrunfiles: cs1Data_A_b_xGT 279 args: -tao_monitor -tao_max_it 100 -tao_type brgn -tao_brgn_regularization_type l2pure -tao_brgn_regularizer_weight 1e-8 -tao_gatol 1.e-6 280 281 TEST*/ 282