5
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.