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*/