Actual source code: cs1.c

  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    - system routines        petscmat.h - matrices
  7:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
  8:      petscviewer.h - viewers               petscpc.h  - preconditioners

 10: */

 12: #include <petsctao.h>

 14: /*
 15: Description:   Compressive sensing test example 1.
 16:                0.5*||Ax-b||^2 + lambda*||D*x||_1
 17:                Xiang Huang: Nov 19, 2018

 19: Reference:     None
 20: */

 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";

 27: #define M 3
 28: #define N 5
 29: #define K 4

 31: /* User-defined application context */
 32: typedef struct {
 33:   /* Working space. linear least square:  f(x) = A*x - b */
 34:   PetscReal A[M][N]; /* array of coefficients */
 35:   PetscReal b[M];    /* array of observations */
 36:   PetscReal xGT[M];  /* array of ground truth object, which can be used to compare the reconstruction result */
 37:   PetscReal D[K][N]; /* array of coefficients for 0.5*||Ax-b||^2 + lambda*||D*x||_1 */
 38:   PetscReal J[M][N]; /* dense jacobian matrix array. For linear least square, J = A. For nonlinear least square, it is different from A */
 39:   PetscInt  idm[M];  /* Matrix row, column indices for jacobian and dictionary */
 40:   PetscInt  idn[N];
 41:   PetscInt  idk[K];
 42: } AppCtx;

 44: /* User provided Routines */
 45: PetscErrorCode InitializeUserData(AppCtx *);
 46: PetscErrorCode FormStartingPoint(Vec);
 47: PetscErrorCode FormDictionaryMatrix(Mat, AppCtx *);
 48: PetscErrorCode EvaluateFunction(Tao, Vec, Vec, void *);
 49: PetscErrorCode EvaluateJacobian(Tao, Vec, Mat, Mat, void *);

 51: /*--------------------------------------------------------------------*/
 52: int main(int argc, char **argv)
 53: {
 54:   Vec       x, f; /* solution, function f(x) = A*x-b */
 55:   Mat       J, D; /* Jacobian matrix, Transform matrix */
 56:   Tao       tao;  /* Tao solver context */
 57:   PetscInt  i;    /* iteration information */
 58:   PetscReal hist[100], resid[100];
 59:   PetscInt  lits[100];
 60:   AppCtx    user; /* user-defined work context */

 63:   PetscInitialize(&argc, &argv, (char *)0, help);

 65:   /* Allocate solution and vector function vectors */
 66:   VecCreateSeq(PETSC_COMM_SELF, N, &x);
 67:   VecCreateSeq(PETSC_COMM_SELF, M, &f);

 69:   /* Allocate Jacobian and Dictionary matrix. */
 70:   MatCreateSeqDense(PETSC_COMM_SELF, M, N, NULL, &J);
 71:   MatCreateSeqDense(PETSC_COMM_SELF, K, N, NULL, &D); /* XH: TODO: dense -> sparse/dense/shell etc, do it on fly  */

 73:   for (i = 0; i < M; i++) user.idm[i] = i;
 74:   for (i = 0; i < N; i++) user.idn[i] = i;
 75:   for (i = 0; i < K; i++) user.idk[i] = i;

 77:   /* Create TAO solver and set desired solution method */
 78:   TaoCreate(PETSC_COMM_SELF, &tao);
 79:   TaoSetType(tao, TAOBRGN);

 81:   /* User set application context: A, D matrice, and b vector. */
 82:   InitializeUserData(&user);

 84:   /* Set initial guess */
 85:   FormStartingPoint(x);

 87:   /* Fill the content of matrix D from user application Context */
 88:   FormDictionaryMatrix(D, &user);

 90:   /* Bind x to tao->solution. */
 91:   TaoSetSolution(tao, x);
 92:   /* Bind D to tao->data->D */
 93:   TaoBRGNSetDictionaryMatrix(tao, D);

 95:   /* Set the function and Jacobian routines. */
 96:   TaoSetResidualRoutine(tao, f, EvaluateFunction, (void *)&user);
 97:   TaoSetJacobianResidualRoutine(tao, J, J, EvaluateJacobian, (void *)&user);

 99:   /* Check for any TAO command line arguments */
100:   TaoSetFromOptions(tao);

102:   TaoSetConvergenceHistory(tao, hist, resid, 0, lits, 100, PETSC_TRUE);

104:   /* Perform the Solve */
105:   TaoSolve(tao);

107:   /* XH: Debug: View the result, function and Jacobian.  */
108:   PetscPrintf(PETSC_COMM_SELF, "-------- result x, residual f=A*x-b, and Jacobian=A. -------- \n");
109:   VecView(x, PETSC_VIEWER_STDOUT_SELF);
110:   VecView(f, PETSC_VIEWER_STDOUT_SELF);
111:   MatView(J, PETSC_VIEWER_STDOUT_SELF);
112:   MatView(D, PETSC_VIEWER_STDOUT_SELF);

114:   /* Free TAO data structures */
115:   TaoDestroy(&tao);

117:   /* Free PETSc data structures */
118:   VecDestroy(&x);
119:   VecDestroy(&f);
120:   MatDestroy(&J);
121:   MatDestroy(&D);

123:   PetscFinalize();
124:   return 0;
125: }

127: /*--------------------------------------------------------------------*/
128: PetscErrorCode EvaluateFunction(Tao tao, Vec X, Vec F, void *ptr)
129: {
130:   AppCtx          *user = (AppCtx *)ptr;
131:   PetscInt         m, n;
132:   const PetscReal *x;
133:   PetscReal       *b = user->b, *f;

135:   VecGetArrayRead(X, &x);
136:   VecGetArray(F, &f);

138:   /* Even for linear least square, we do not direct use matrix operation f = A*x - b now, just for future modification and compatibility for nonlinear least square */
139:   for (m = 0; m < M; m++) {
140:     f[m] = -b[m];
141:     for (n = 0; n < N; n++) f[m] += user->A[m][n] * x[n];
142:   }
143:   VecRestoreArrayRead(X, &x);
144:   VecRestoreArray(F, &f);
145:   PetscLogFlops(2.0 * M * N);
146:   return 0;
147: }

149: /*------------------------------------------------------------*/
150: /* J[m][n] = df[m]/dx[n] */
151: PetscErrorCode EvaluateJacobian(Tao tao, Vec X, Mat J, Mat Jpre, void *ptr)
152: {
153:   AppCtx          *user = (AppCtx *)ptr;
154:   PetscInt         m, n;
155:   const PetscReal *x;

157:   VecGetArrayRead(X, &x); /* not used for linear least square, but keep for future nonlinear least square) */
158:   /* 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?
159:     For nonlinear least square, we require x to compute J, keep codes here for future nonlinear least square*/
160:   for (m = 0; m < M; ++m) {
161:     for (n = 0; n < N; ++n) user->J[m][n] = user->A[m][n];
162:   }

164:   MatSetValues(J, M, user->idm, N, user->idn, (PetscReal *)user->J, INSERT_VALUES);
165:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
166:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);

168:   VecRestoreArrayRead(X, &x); /* not used for linear least square, but keep for future nonlinear least square) */
169:   PetscLogFlops(0);                      /* 0 for linear least square, >0 for nonlinear least square */
170:   return 0;
171: }

173: /* ------------------------------------------------------------ */
174: /* Currently fixed matrix, in future may be dynamic for D(x)? */
175: PetscErrorCode FormDictionaryMatrix(Mat D, AppCtx *user)
176: {
177:   MatSetValues(D, K, user->idk, N, user->idn, (PetscReal *)user->D, INSERT_VALUES);
178:   MatAssemblyBegin(D, MAT_FINAL_ASSEMBLY);
179:   MatAssemblyEnd(D, MAT_FINAL_ASSEMBLY);

181:   PetscLogFlops(0); /* 0 for fixed dictionary matrix, >0 for varying dictionary matrix */
182:   return 0;
183: }

185: /* ------------------------------------------------------------ */
186: PetscErrorCode FormStartingPoint(Vec X)
187: {
188:   VecSet(X, 0.0);
189:   return 0;
190: }

192: /* ---------------------------------------------------------------------- */
193: PetscErrorCode InitializeUserData(AppCtx *user)
194: {
195:   PetscReal *b = user->b; /* **A=user->A, but we don't know the dimension of A in this way, how to fix? */
196:   PetscInt   m, n, k;     /* loop index for M,N,K dimension. */

198:   /* b = A*x while x = [0;0;1;0;0] here*/
199:   m      = 0;
200:   b[m++] = 0.28;
201:   b[m++] = 0.55;
202:   b[m++] = 0.96;

204:   /* matlab generated random matrix, uniformly distributed in [0,1] with 2 digits accuracy. rng(0); A = rand(M, N); A = round(A*100)/100;
205:   A = [0.81  0.91  0.28  0.96  0.96
206:        0.91  0.63  0.55  0.16  0.49
207:        0.13  0.10  0.96  0.97  0.80]
208:   */
209:   m               = 0;
210:   n               = 0;
211:   user->A[m][n++] = 0.81;
212:   user->A[m][n++] = 0.91;
213:   user->A[m][n++] = 0.28;
214:   user->A[m][n++] = 0.96;
215:   user->A[m][n++] = 0.96;
216:   ++m;
217:   n               = 0;
218:   user->A[m][n++] = 0.91;
219:   user->A[m][n++] = 0.63;
220:   user->A[m][n++] = 0.55;
221:   user->A[m][n++] = 0.16;
222:   user->A[m][n++] = 0.49;
223:   ++m;
224:   n               = 0;
225:   user->A[m][n++] = 0.13;
226:   user->A[m][n++] = 0.10;
227:   user->A[m][n++] = 0.96;
228:   user->A[m][n++] = 0.97;
229:   user->A[m][n++] = 0.80;

231:   /* initialize to 0 */
232:   for (k = 0; k < K; k++) {
233:     for (n = 0; n < N; n++) user->D[k][n] = 0.0;
234:   }
235:   /* Choice I: set D to identity matrix of size N*N for testing */
236:   /* for (k=0; k<K; k++) user->D[k][k] = 1.0; */
237:   /* Choice II: set D to Backward difference matrix of size (N-1)*N, with zero extended boundary assumption */
238:   for (k = 0; k < K; k++) {
239:     user->D[k][k]     = -1.0;
240:     user->D[k][k + 1] = 1.0;
241:   }

243:   return 0;
244: }

246: /*TEST

248:    build:
249:       requires: !complex !single !quad !defined(PETSC_USE_64BIT_INDICES)

251:    test:
252:       localrunfiles: cs1Data_A_b_xGT
253:       args: -tao_smonitor -tao_max_it 100 -tao_type pounders -tao_gatol 1.e-6

255:    test:
256:       suffix: 2
257:       localrunfiles: cs1Data_A_b_xGT
258:       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 -tao_brgn_subsolver_tao_bnk_ksp_converged_reason

260:    test:
261:       suffix: 3
262:       localrunfiles: cs1Data_A_b_xGT
263:       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

265:    test:
266:       suffix: 4
267:       localrunfiles: cs1Data_A_b_xGT
268:       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

270:    test:
271:       suffix: 5
272:       localrunfiles: cs1Data_A_b_xGT
273:       args: -tao_monitor -tao_max_it 100 -tao_type brgn -tao_brgn_regularization_type lm -tao_gatol 1.e-6 -tao_brgn_subsolver_tao_type bnls

275: TEST*/