# Making sense of standard deviation after sampling using Cholesky

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:

Figure 2: