Due: Friday, Dec. 4.
You have been using prior models since Lab 3. A prior is a probability distribution over images, and you have written down its density, differentiated it, and minimized things containing it. You have never once looked at an image drawn from one.
This lab draws them. You will sample the GMRF prior of Lab 3 and the QGGMRF prior of Lab 4, and see what those models actually believe an image looks like. Then you will sample a posterior, which gives something the MAP estimate can never give: a pixel by pixel map of how sure the answer is.
The background is Chapter 15, on the Metropolis sampler, the Hastings-Metropolis sampler, and the Gibbs sampler. The notation here matches the book.
Only one file, and only in Module 3. Modules 1 and 2 need no data at all, because a sampler makes its own images out of random numbers.
Source: image 23 of the Kodak Lossless True Color Image Suite, released for unrestricted use.
wget https://cabouman.github.io/grad_labs/labs/resources/kodim23.pgm
lab7.py, holding functions and nothing else. Run your
experiments from a separate script.float64. No for loop
over pixels. This rule is harder to keep here than in any
earlier lab, because samplers are usually written one pixel at a
time. Step 3 shows how to keep it.The GMRF prior of Lab 3 is the density
$$ p(x) = \frac{1}{z} \exp \left\{ - u(x) \right\} \ , \qquad u(x) = \frac{1}{2\sigma_x^2} \sum_{ \{s,r\} \in \mathcal{P} } g_{s-r} \, ( x_s - x_r )^2 \ , \tag{1} $$with the same \(3 \times 3\) array \(g\) and the same 8-point neighborhood as Lab 3. The function \(u\) is called the energy. Low energy means high probability.
Sampling from (1) directly is out of reach. The normalizing constant \(z\) is an integral over every image there is, and nobody can compute it. What is easy is the conditional distribution of one pixel given all the others, because everything except that one pixel is a constant. The samplers in this chapter are built entirely out of that observation.
Collect the terms of (1) that contain \(x_s\). They are the pairs that join \(s\) to one of its neighbors, so
$$ u(x) = \frac{1}{2\sigma_x^2} \sum_{ r \in \partial s } g_{s-r} ( x_s - x_r )^2 + \mbox{terms without } x_s \ . $$That is a quadratic in \(x_s\), so the conditional density is Gaussian. Completing the square gives
$$ X_s \, | \, \{ X_r \}_{r \not= s} \ \sim \ N \! \left( \mu_s , \ \frac{\sigma_x^2}{c_s} \right) \ , \qquad \mu_s = \frac{1}{c_s} \sum_{ r \in \partial s } g_{s-r} \, x_r \ , \tag{2} $$where \( c_s = \sum_{r \in \partial s} g_{s-r} \) is the same normalization you computed in Lab 3, equal to 1 for an interior pixel and less along an edge.
The Gibbs sampler is now just this: visit the pixels in some order, and replace each one with a draw from (2). Repeat. The chapter proves that the distribution of the whole image converges to \(p(x)\), whatever you started from.
Written literally, the Gibbs sampler updates one pixel at a time, and a 256 by 256 image needs 65,536 updates per sweep. In Python that is unusably slow, and it breaks the rule this course has kept since Lab 2.
There is a way out, and it comes from the neighborhood. Two pixels can be updated at the same time if neither is a neighbor of the other, because then neither appears in the other's conditional. Color the lattice by the parity of the row and column index:
| 0 | 1 | 0 | 1 |
| 2 | 3 | 2 | 3 |
| 0 | 1 | 0 | 1 |
| 2 | 3 | 2 | 3 |
Two pixels of the same color differ by an even number of rows and an even number of columns, so they are never 8-neighbors. Every pixel of one color can therefore be drawn simultaneously, as one array operation. Four of those makes a full sweep.
def neighbor_stats(x):
"""Return the neighbor sums and neighbor weights of the GMRF prior.
Uses the fixed 3 x 3 array g and the 8-point neighborhood of Lab 3,
with zero padding, so that pixels on the boundary have fewer
neighbors.
x (H, W) float64 tensor
returns (nbr_sum, c), each (H, W) float64, where
nbr_sum[s] = sum over r in ds of g[s-r] x[r]
c[s] = sum over r in ds of g[s-r]
"""
def gibbs_sweep_gmrf(x, sigma_x):
"""Run one Gibbs sweep for the GMRF prior of equation (1).
Updates all four color classes in turn, each as one array operation,
drawing each pixel from the conditional of equation (2).
x (H, W) float64 tensor, the current sample
sigma_x float
returns (H, W) float64 tensor, the sample after one sweep
"""
Task 3a. Draw from the prior. Start from \(x = 0\) on a 256 by 256 array and run 200 sweeps with \(\sigma_x = 0.05\). Display the result. Then repeat for \(\sigma_x = 0.01\) and \(\sigma_x = 0.2\). Display each with its own gray scale, chosen to fill the range, and state the range you used.
The three sampled images, one for each \(\sigma_x\), with their display ranges. Add two sentences: what changed with \(\sigma_x\), and one specific way in which these images do not resemble a photograph.
Task 3b. Watch the brightness wander. Record the mean of the whole image after every sweep, for 500 sweeps, and plot it. Then plot the standard deviation of \(x - \mathrm{mean}(x)\) over the same sweeps.
The two plots from Task 3b, and two sentences saying which quantity converged, which did not, and which part of equation (1) is responsible. Say what you did about it before displaying the images in D1.
Now change the prior to the QGGMRF of Lab 4. The energy becomes
$$ u(x) = \sum_{ \{s,r\} \in \mathcal{P} } b_{s,r} \, \rho( x_s - x_r ) \ , \qquad b_{s,r} = g_{s-r} \ , \tag{3} $$with the QGGMRF potential from Lab 4, repeated here:
$$ \rho( \Delta ) = \frac{ |\Delta|^p }{ p\,\sigma_x^p } \left( \frac{ \left| \frac{\Delta}{T \sigma_x} \right|^{q-p} } { 1 + \left| \frac{\Delta}{T \sigma_x} \right|^{q-p} } \right) \ . \tag{4} $$Repeat Step 2 with this energy. Collect the terms containing \(x_s\) and you get a sum of potentials, not a quadratic. It is a perfectly good 1-D density, but it is not one of the distributions your random number generator knows how to draw from, and there is no formula for its inverse CDF. The Gibbs sampler has nothing to draw.
Use \(q = 2\), \(p = 1.2\), and \(T = 0.1\), the baseline values from Lab 4. With those, equations (3) and (4) are exactly the prior term of the Lab 4 cost function.
The Metropolis sampler needs far less than the Gibbs sampler does. It never draws from the conditional. It only ever compares two energies.
Pick a pixel \(s\). Propose a new value by adding noise,
$$ W_s = X_s + Z \ , \qquad Z \sim N( 0, \sigma_{\mathrm{prop}}^2 ) \ , \tag{5} $$compute the change in energy \( \Delta E = u(W) - u(X) \), and accept the proposal with probability
$$ \alpha = \min \left\{ 1 , \ \exp( - \Delta E ) \right\} \ . \tag{6} $$If the proposal lowers the energy it is always accepted. If it raises the energy it is sometimes accepted, and how often depends on how much it raised it. Those occasional uphill moves are the whole reason this produces samples from \(p(x)\) rather than sliding into a minimum of \(u\) and stopping.
Because the proposal (5) is symmetric, meaning \( q(w|x) = q(x|w) \), the acceptance rule needs no correction term and the plain Metropolis form (6) applies. The change in energy is local:
$$ \Delta E = \sum_{ r \in \partial s } b_{s,r} \left[ \rho( w_s - x_r ) - \rho( x_s - x_r ) \right] \ . \tag{7} $$The four-color trick from Step 3 works here too, and for the same reason. Pixels of one color are not neighbors, so their \(\Delta E\) values do not depend on each other and all of them can be proposed, scored, and accepted or rejected in one array operation.
def qggmrf_delta_energy(x, w, sigma_x, p, q, T):
"""Return the per-pixel energy change of equation (7).
Computes, for every pixel s at once, the change in the QGGMRF energy
of equation (3) that would result from replacing x[s] by w[s] with
every other pixel held at its value in x.
x, w (H, W) float64 tensors, current and proposed values
returns (H, W) float64 tensor
"""
def metropolis_sweep_qggmrf(x, sigma_x, p, q, T, sigma_prop):
"""Run one coordinatewise Metropolis sweep for the QGGMRF prior.
Updates all four color classes in turn. For each class, proposes
with equation (5), scores with equation (7), and accepts with
equation (6).
x (H, W) float64 tensor, the current sample
sigma_prop float, the standard deviation of the proposal
returns (x_new, accept_rate), the sample after one sweep and the
fraction of the pixels whose proposal was accepted
"""
The proposal width \(\sigma_{\mathrm{prop}}\) does not change the distribution being sampled. Equation (6) guarantees that for any width. What it changes is how fast the sampler explores, and it fails in both directions. Too small, and almost every proposal is accepted but each one moves almost nowhere. Too large, and the proposals are ambitious but nearly all are rejected, so the image sits still.
Task 6a. Starting from \(x = 0\) on a 256 by 256 array with \(\sigma_x = 0.05\), run 300 sweeps for each of \( \sigma_{\mathrm{prop}} \in \{ 0.002,\ 0.02,\ 0.2 \} \). For each run, record the acceptance rate and the energy \(u(x)\) after every sweep.
The three sampled images side by side, with the acceptance rate of each printed underneath.
One plot of energy versus sweep with all three runs on it, and three sentences: which width worked best, what went wrong in each of the other two, and what acceptance rate you would aim for next time.
Task 6b. Using your best \(\sigma_{\mathrm{prop}}\), draw one QGGMRF sample and place it beside the GMRF sample from D1 at the same \(\sigma_x = 0.05\).
The two samples side by side, GMRF and QGGMRF, with two sentences on how they differ. Relate the difference to the shape of \(\rho\) in equation (4) for large \(|\Delta|\), and to the reason Lab 4 gave for preferring the QGGMRF.
So far the samplers have drawn from priors, which is to say from imagination. Now condition on a measurement. Take the denoising problem of Lab 3, with \(A = I\), and form \( y = x + w \) with \(\sigma_w = 0.05\). The posterior energy is the MAP cost:
$$ u(x) = \frac{1}{2\sigma_w^2} \| y - x \|^2 + \frac{1}{2\sigma_x^2} \sum_{ \{s,r\} \in \mathcal{P} } g_{s-r} ( x_s - x_r )^2 \ . \tag{8} $$Both terms are quadratic in \(x_s\), so the conditional is Gaussian again and the Gibbs sampler is back in business. Completing the square as in Step 2 gives
$$ X_s \, | \, \{ X_r \}_{r \not= s} \ \sim \ N \! \left( \frac{ \frac{y_s}{\sigma_w^2} + \frac{1}{\sigma_x^2} \sum_{r \in \partial s} g_{s-r} x_r } { \frac{1}{\sigma_w^2} + \frac{c_s}{\sigma_x^2} } \ , \ \left( \frac{1}{\sigma_w^2} + \frac{c_s}{\sigma_x^2} \right)^{\!\!-1} \right) \ . \tag{9} $$def gibbs_sweep_posterior(x, y, sigma_w, sigma_x):
"""Run one Gibbs sweep for the denoising posterior of equation (8).
Updates all four color classes in turn, drawing each pixel from the
conditional of equation (9).
x (H, W) float64 tensor, the current sample
y (H, W) float64 tensor, the noisy data
returns (H, W) float64 tensor
"""
def posterior_statistics(y, sigma_w, sigma_x, num_burn, num_keep):
"""Sample the posterior and accumulate its first two moments.
Runs num_burn sweeps and discards them, then runs num_keep more,
accumulating the pixelwise mean and standard deviation over the kept
sweeps.
returns (mean, std, energy_history), two (H, W) float64 tensors and
a list of length num_burn + num_keep
"""
Task 7a. Use a 256 by 256 crop of kodim23 as the clean image, with \(\sigma_w = 0.05\) and \(\sigma_x = 0.05\). Run 100 burn-in sweeps and 400 kept sweeps. Produce the posterior mean image and the pixelwise standard deviation image. Also compute the MAP estimate for the same problem with your Lab 3 code.
Five images: the clean crop, the noisy \(y\), one single posterior sample, the posterior mean, and the MAP estimate, all over \([0,1]\). Report the RMSE of the last three against the clean image.
The pixelwise standard deviation image, displayed with a gray scale chosen to show its structure, with the range stated. Add two sentences saying where in the image the posterior is least certain, and why that is where you would expect it to be.
The theory says the distribution converges. It does not say when. Every sweep before convergence is contaminated by wherever you started, which is why Task 7a threw the first hundred away. A hundred was a guess, and guesses should be checked.
Task 8a. Run the posterior sampler twice more, for 500 sweeps each, once starting from \(x = y\) and once starting from \(x = 0\). Plot the energy \(u(x)\) versus sweep for both on one figure.
The energy plot for the two starting points, and three sentences: roughly how many sweeps the two took to become indistinguishable, why agreement between two different starting points is better evidence than either curve flattening on its own, and whether your choice of 100 burn-in sweeps was enough.
Answer each in a few sentences, using your own images and numbers as evidence.
Your answers to the five questions above.
Submit two files through Brightspace.
lab7.py, as a plain .py file. Do
not paste the code into the report and do not send a zip
file.
Also commit your code to your image-processing-labs
repository under lab7 and push it before you submit.
Concise and clear beats long.
You have now seen the models you have been using since Lab 3, and you have produced something no optimizer can produce: an estimate that says how sure it is. Every reconstruction you have made in this course was a single point picked out of a distribution. This lab is the one where you looked at the distribution.