I have an inverse problem with over 65,000 degrees of freedom. I am using Bayesian formulation to solve this problem. After using the optimization algorithm to obtain MAP solution, I want to calculate uncertainty associated with this solution.

For which, I calculate the covariance matrix by inverting the Hessian of the objective function. Since I cannot calculate the exact Hessian, I use an approximate Hessian given by $J^T J$ where $J$ is the Jacobian matrix. Then, the Hessian, after including the prior, is given by $H = J^J/\sigma^2 + I/\alpha^2$, where $I$ is an identity matrix and $\alpha$ is the standard deviation of prior, and $\sigma$ is the standard deviation of the data that I am using.

Now, after calculating the covariance I sample the outcomes by using the following:

$m = m_\mathrm{map} + R's$

Where $R$ is obtained by Cholesky decomposition of the covariance matrix and $s$ is supposed to be a vector containing normally distributed random numbers.

Now, the issue is that the outcome of the sampling does not seem to make any sense, Figure 1 shows the standard deviation. I would expect some kind of smooth spatial distribution as my $m$ varies locally. I am also attaching plot for $R' \mathrm{ones}(n,1)$ (essentially assuming all the elements of $s$ to be unity), Figure 2. I was expecting the standard deviation from the sampling to follow this trend. Am I missing something here?

Also, I am sure that 10000 samples are enough as all the m(i) have normal distribution.

I might have missed some vital information, please let me know if further description is required.

Note, that m = m(x,y) and the axis in the figures represent x and y coordinates.

Figure 1:

Standard deviation of m obtained after sampling for 10000 times

Figure 2:

R's, where s is a vector containing all ones.



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.