Codes for "Sequential Monte Carlo with Gaussian Mixture Approximation for Infinite-Dimensional Statistical Inverse Problems"

Overview

The program relies on FEniCSx (Version 2024), sklearn (Version 1.5), and mpi4py.

Directory core contains the main functions and classes that are useful for implementing the algorithms. Specifically,

1. probability.py: This file contains classes of GaussianElliptic2[The Gaussian measure implemented by finite element methods based on solving elliptic differential equations used for generating samples and also contains the functionality of evaluate the gradient and Hessian operators]; GaussianFiniteRank[The Gaussian measure implemented by finite element methods and eigensystem decomposition].

• noise.py: This file contains the class NoiseGaussianIID.
• model.py: This file contains the class Domain and two classes Domain2D and Domain1D inherit from the parent class Domain; contains the parent class ModelBase of the model classes employed in specific examples, the class ModelBase incorporates the components of the domain, prior, equation solver, and noise.
• linear_eq_solver.py: contains the function cg_my which is our implementation of the conjugate gradient algorithm for solving linear equations.
• eigensystem.py: This file contains the function double_pass which is our implementation of an algorithm for calculating eigensystem.
• approximate_sample.py: This file contains the class LaplaceApproximate, which can be used to compute the Laplace approximation of the posterior measures.
• optimizer.py: This file contains the class OptimBase[incorporate an implementation of armijo_line_search can be employed for each optimizer]; the class GradientDescent[an implementation of the gradient descent algorithm]; the class NewtonCG[an implementation of the Newton conjugate gradient algorithm].
• sample.py: This file contains the class pCN, which is a type of discrete invariant Markov chain Monte Carlo sampling algorithm.
• Plotting.py: This file contains some functions that can draw functions generated by FEniCS.
• misc.py: This file contains functions of trans2spnumpy, trans2sptorch, spnumpy2sptorch, sptorch2spnumpy, and sptensor2cude, which will be useful for transferring sparse matrixes to different forms required for doing calculations in numpy, pytorch, and FEniCS. This file also contains the function construct_measurement_matrix, which will be used for generating a sparse matrix S. The matrix S times a function generated by FEniCS to get the values at the measurement points.

Directory DarcyFlow contains the main functions and classes for the Darcy flow problem.

1. common.py: This file contains the class EquSolver and the class ModelDarcyFlow. The class EquSolver contains the implementations of solvers of forward, adjoint, incremental forward, and incremental adjoint equations. These equations are necessary for implementing gradient and Newton-type optimization algorithms. The class ModelDarcyFlow contains the function of calculating the loss, the gradient, and the Hessian operator.

• generate_data.py The file first samples from the prior and saves them. Then, it solves the PDE with the function as its parameter and saves the measurements of the solution function at specified locations. For the convenience of the readers, we have prepared the data in the folders ./SequentialMonteCarlo/data and ./SequentialMonteCarlo/data_sparse_measure.

The subdirectory SequentialMonteCarlo contains the main files for our numerical results.

1. example1SMCGM: Used as the first example in the paper, solving the inverse problem with GMM approximation.

• example1SMCpCN: Used as a comparative example for the first example in the paper, solving the inverse problem with preconditioned Crank-Nicolson (pCN) method.
• example2SMCGM: Used as the second example in the paper, solving the Darcy flow inverse problem with SMCGM method.
• example2SMCpCN: Used as a comparative example for the second example in the paper, solving the Darcy flow inverse problem with SMC-pCN method.
• example2ResultAnalysis: Used for result analysis of the second example in the paper, plotting the densities of various Fourier coefficients for comparison of the posterior.
• example2_4_1MeshIndpdt.py: Solves the Darcy flow inverse problem using SMC-GM.
• example2_4_2MeshDpdt.py: Solves the Darcy flow inverse problem using SMC-RW.
• example2_4_3PlotMeshIndpdt.py: Reads the results from 2_4_1 and 2_4_2 and plots them for comparison.
• example2_4_0MeshIndpdt.sh: For the reader's convenience, we have packaged 2_4_1 to 2_4_3 into a .sh file, which will automatically execute the above three commands, thereby completing the verification of dimension independence.
• example3_1SMCGM.py: Solves the Darcy flow inverse problem using SMCGM under sparse measurement conditions.
• example3_2SMCpCN.py: Solves the Darcy flow inverse problem using SMCpCN under sparse measurement conditions.
• example3_3ResultAnalysis.py: Reads the data from 3_1 and 3_2, and plots them for comparison.
• Mix_Gaussian.py: Contains the class Mix_Gauss_IS_proposal, which can estimate a mixture of Gaussians from data and sample from it, used in the SMCGM algorithm.

Workflows

The first example

Run the following command sequentially for the first example:

python example1SMCGM.py

python example1SMCpCN.py

When all of the commands are executed, you will find a directory named "RESULTS". In the directories ./results/example1SMCGM and ./results/example1SMCpCN, there are images of the 1000 samples spanning from the first to the last layer of the SMCGM and SMCpCN algorithms. These images fully illustrate the execution process of the algorithms, with the final layer images being used in the paper.

The second example

Run the following command sequentially for the second example:

mpiexec -n 20 python example2_1SMCGM.py

mpiexec -n 20 python example2_2SMCpCN.py

mpiexec -n 20 python example2_3ResultAnalysis.py

Run the following code to verify the dimension independence.

./example2_4_0MeshIndpdt.sh

In paths ./results/example2SMCGM and ./results/example2SMCpCN, there are the mean values of each layer for SMCGM and SMCpCN solutions of the Darcy flow inverse problem, respectively. The mean value of the last layer is displayed in the paper. Additionally, in path ./results/example2Analysis, there are comparisons of the density functions of the Fourier coefficients estimated from the samples using both methods, and the results of the dimension independence verification are also in this directory.

The third example

Run the following command sequentially for the third example:

mpiexec -n 20 python example3SMCGM.py

mpiexec -n 20 python example3SMCpCN.py

In paths ./results/example3SMCGM and ./results/example3SMCpCN, there are the mean values of each layer for SMCGM and SMCpCN solutions of the Darcy flow inverse problem, respectively. The mean values of the last layer is displayed in the paper. Additionally, in path ./results/example3Analysis, we present a comparison of the solutions obtained by substituting these means into the forward problem with the true data.