Actual source code: cs1.c
petsc-3.14.6 2021-03-30
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
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";
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*/
40: #define M 3
41: #define N 5
42: #define K 4
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;
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 *);
64: /*--------------------------------------------------------------------*/
65: int main(int argc,char **argv)
66: {
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 */
76: PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr;
78: /* Allocate solution and vector function vectors */
79: VecCreateSeq(PETSC_COMM_SELF,N,&x);
80: VecCreateSeq(PETSC_COMM_SELF,M,&f);
82: /* Allocate Jacobian and Dictionary matrix. */
83: MatCreateSeqDense(PETSC_COMM_SELF,M,N,NULL,&J);
84: MatCreateSeqDense(PETSC_COMM_SELF,K,N,NULL,&D); /* XH: TODO: dense -> sparse/dense/shell etc, do it on fly */
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;
90: /* Create TAO solver and set desired solution method */
91: TaoCreate(PETSC_COMM_SELF,&tao);
92: TaoSetType(tao,TAOBRGN);
94: /* User set application context: A, D matrice, and b vector. */
95: InitializeUserData(&user);
97: /* Set initial guess */
98: FormStartingPoint(x);
100: /* Fill the content of matrix D from user application Context */
101: FormDictionaryMatrix(D,&user);
103: /* Bind x to tao->solution. */
104: TaoSetInitialVector(tao,x);
105: /* Bind D to tao->data->D */
106: TaoBRGNSetDictionaryMatrix(tao,D);
108: /* Set the function and Jacobian routines. */
109: TaoSetResidualRoutine(tao,f,EvaluateFunction,(void*)&user);
110: TaoSetJacobianResidualRoutine(tao,J,J,EvaluateJacobian,(void*)&user);
112: /* Check for any TAO command line arguments */
113: TaoSetFromOptions(tao);
115: TaoSetConvergenceHistory(tao,hist,resid,0,lits,100,PETSC_TRUE);
117: /* Perform the Solve */
118: TaoSolve(tao);
120: /* XH: Debug: View the result, function and Jacobian. */
121: PetscPrintf(PETSC_COMM_SELF, "-------- result x, residual f=A*x-b, and Jacobian=A. -------- \n");
122: VecView(x,PETSC_VIEWER_STDOUT_SELF);
123: VecView(f,PETSC_VIEWER_STDOUT_SELF);
124: MatView(J,PETSC_VIEWER_STDOUT_SELF);
125: MatView(D,PETSC_VIEWER_STDOUT_SELF);
127: /* Free TAO data structures */
128: TaoDestroy(&tao);
130: /* Free PETSc data structures */
131: VecDestroy(&x);
132: VecDestroy(&f);
133: MatDestroy(&J);
134: MatDestroy(&D);
136: PetscFinalize();
137: return ierr;
138: }
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;
150: VecGetArrayRead(X,&x);
151: VecGetArray(F,&f);
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: VecRestoreArrayRead(X,&x);
161: VecRestoreArray(F,&f);
162: PetscLogFlops(2.0*M*N);
163: return(0);
164: }
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;
176: VecGetArrayRead(X,&x); /* 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: }
185: MatSetValues(J,M,user->idm,N,user->idn,(PetscReal *)user->J,INSERT_VALUES);
186: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
187: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
189: VecRestoreArrayRead(X,&x);/* 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: return(0);
192: }
194: /* ------------------------------------------------------------ */
195: /* Currently fixed matrix, in future may be dynamic for D(x)? */
196: PetscErrorCode FormDictionaryMatrix(Mat D,AppCtx *user)
197: {
201: MatSetValues(D,K,user->idk,N,user->idn,(PetscReal *)user->D,INSERT_VALUES);
202: MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY);
203: MatAssemblyEnd(D,MAT_FINAL_ASSEMBLY);
205: PetscLogFlops(0); /* 0 for fixed dictionary matrix, >0 for varying dictionary matrix */
206: return(0);
207: }
209: /* ------------------------------------------------------------ */
210: PetscErrorCode FormStartingPoint(Vec X)
211: {
214: VecSet(X,0.0);
215: return(0);
216: }
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. */
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;
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;
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: }
254: return(0);
255: }
257: /*TEST
259: build:
260: requires: !complex !single !quad !define(PETSC_USE_64BIT_INDICES)
262: test:
263: localrunfiles: cs1Data_A_b_xGT
264: args: -tao_smonitor -tao_max_it 100 -tao_type pounders -tao_gatol 1.e-6
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 -tao_brgn_subsolver_ksp_converged_reason
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
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
281: test:
282: suffix: 5
283: localrunfiles: cs1Data_A_b_xGT
284: 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
286: TEST*/