Due: Friday, Oct. 2.
Your Lab 3 reconstruction had one clear defect. The prior was quadratic, so it pulled hardest on the pixel pairs that differed the most, and the pairs that differ the most are the edges. The restoration worked, and it blurred every edge in the picture.
In this lab you replace the quadratic prior with a non-Gaussian one, the QGGMRF, which leaves edges alone. That fix has a price. The cost function stops being quadratic, and the gradient descent you have been using stops working well on it.
The way out is majorization. At each outer iteration you build a quadratic function that sits on top of the true cost, minimize that instead, and repeat. Every outer iteration is then the quadratic problem you already solved in Lab 3, and the true cost is guaranteed to go down at every step without ever being evaluated. Module 2 is where that happens, and it is the part worth reading twice.
The background is Chapter 6 for non-Gaussian MRF priors and the QGGMRF potential, and Chapter 8 for surrogate functions and majorization. The notation here matches the book.
The same image and blur kernel as Lab 3, so that you can compare directly with what you already have.
Source: image 23 of the Kodak Lossless True Color Image Suite, released for unrestricted use.
Source: A. Levin, Y. Weiss, F. Durand, and W. T. Freeman, "Understanding and Evaluating Blind Deconvolution Algorithms," IEEE CVPR 2009.
wget https://cabouman.github.io/grad_labs/labs/resources/kodim23.pgm wget https://cabouman.github.io/grad_labs/labs/resources/levin09_kernels.zip unzip levin09_kernels.zip
lab4.py, holding functions and nothing else. Run
your experiments from a separate script.forward, adjoint,
load_pgm, and load_kernel from your
lab3.py. Those four are done and you should not write
them again.float64, so that the cost differences
you plot are not lost in rounding. No for loop over
pixels.Write the MAP cost function in the general pairwise form
$$ f(x) = \frac{1}{2\sigma_w^2} \| y - A x \|^2 + \sum_{\{s,r\} \in \mathcal{P}} b_{s,r}\, \rho\!\left( x_s - x_r \right) \ , \tag{1} $$where \( \mathcal{P} \) is the set of neighboring pixel pairs, \( b_{s,r} \) are the neighbor weights, and \( \rho \) is the potential function. In Lab 3 you used the quadratic potential
$$ \rho( \Delta ) = \frac{\Delta^2}{2\sigma_x^2} \ , \qquad \rho^\prime( \Delta ) = \frac{\Delta}{\sigma_x^2} \ . \tag{2} $$The derivative \( \rho^\prime \) is called the influence function, because it says how hard the prior pulls two neighboring pixels together when their difference is \( \Delta \).
Now look at what (2) does. Its influence grows without bound. The larger the difference between two neighboring pixels, the harder the prior pulls them together. But a real edge is exactly a place where two neighboring pixels ought to differ a lot. The quadratic prior therefore fights hardest precisely where it should not, and the edges lose.
What we want instead is a potential that behaves like a quadratic for small differences, so that noise in flat regions is still suppressed, but whose influence flattens out for large differences, so that edges are left alone.
The QGGMRF potential does exactly that. In the scaled form from Chapter 6,
$$ \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{3} $$with influence function
$$ \rho^\prime( \Delta ) = \frac{ |\Delta|^{p-1} }{ \sigma_x^p }\, \frac{ \left| \frac{\Delta}{T\sigma_x} \right|^{q-p} \left( \frac{q}{p} + \left| \frac{\Delta}{T\sigma_x} \right|^{q-p} \right) } { \left( 1 + \left| \frac{\Delta}{T\sigma_x} \right|^{q-p} \right)^2 } \,\mathrm{sign}( \Delta ) \ . \tag{4} $$There are three parameters plus the scale \( \sigma_x \), and each one does a single job.
Requiring \( 1 \le p \lt q \) keeps \( \rho \) convex. With \( q = 2 \) the potential has a continuous, bounded second derivative at \( \Delta = 0 \), which Module 2 needs. The plain GGMRF potential \( |\Delta|^p / (p \sigma_x^p) \) does not have that property, since for \( p \lt 2 \) its second derivative blows up at zero.
Write these two functions.
def rho(delta, sigma_x, p, q, T):
"""Return the QGGMRF potential of equation (3), elementwise.
delta float64 tensor of any shape
returns float64 tensor of the same shape
"""
def rho_prime(delta, sigma_x, p, q, T):
"""Return the QGGMRF influence function of equation (4), elementwise.
delta float64 tensor of any shape
returns float64 tensor of the same shape
"""
Task 2a. Plot \( \rho^\prime(\Delta) \) for the quadratic potential (2) and for the QGGMRF (4) on the same axes over \( \Delta \in [-0.3, 0.3] \), using \( p = 1.2 \), \( q = 2 \), \( T = 0.1 \), and \( \sigma_x = 0.2 \).
The plot of the two influence functions, and two sentences using it to say what each prior does to a pixel pair straddling a strong edge.
Put the QGGMRF potential into the cost (1) and two things go wrong for plain gradient descent.
First, there is no longer one curvature. The second derivative of \( \rho \) is about \( 1/(p\,T^{2-p}\sigma_x^2) \) near \( \Delta = 0 \) and falls toward zero for large \( |\Delta| \). Across an image the neighboring pixel differences span both regimes at once, and gradient descent has one step size for the whole image. It must be small enough for the stiffest part of the cost, which makes it slow everywhere else.
Second, the safe step size depends on the current image, so it changes from iteration to iteration. There is no fixed step size you can compute once and trust.
Majorization fixes both. Instead of minimizing the hard function, we repeatedly minimize an easy one built to stand in for it.
Let \( f(x) \) be the function we want to minimize and let \( x^\prime \) be the current image, called the point of approximation. A function \( q(x; x^\prime) \) is a surrogate function for \( f \) at \( x^\prime \) if
$$ f( x^\prime ) = q( x^\prime; x^\prime ) \ , \qquad f( x ) \le q( x; x^\prime ) \quad \text{for all } x \ . \tag{5} $$In words, the surrogate touches the true cost at the current image and lies above it everywhere else. Picture two curves that touch at one point, with the surrogate sitting on top.
In practice we build a \( Q(x; x^\prime) \) that differs from \( q(x; x^\prime) \) by a term depending on \( x^\prime \) but not on \( x \). That changes nothing, because a constant does not move a minimum. Stated in terms of \( Q \), the condition is
$$ f( x ) \ \le \ Q( x; x^\prime ) - Q( x^\prime; x^\prime ) + f( x^\prime ) \quad \text{for all } x, x^\prime \ . \tag{6} $$This is the inequality you will check numerically in Step 6.
The majorization-minimization algorithm is the obvious one. From the current image, build the surrogate, minimize it, and use the result as the new current image:
$$ x^{(k+1)} = \arg\min_{x} \ Q\!\left( x; x^{(k)} \right) \ . \tag{7} $$Here is the whole argument for why it works, in one line.
$$ f\!\left( x^{(k+1)} \right) \ \le \ q\!\left( x^{(k+1)}; x^{(k)} \right) \ \le \ q\!\left( x^{(k)}; x^{(k)} \right) \ = \ f\!\left( x^{(k)} \right) \ . \tag{8} $$Read it left to right. The first inequality is the upper bound in (5). The second holds because \( x^{(k+1)} \) minimizes the surrogate. The equality is the touching condition in (5). Together they give \( f(x^{(k+1)}) \le f(x^{(k)}) \), so the true cost never goes up.
Three consequences are worth stating plainly.
The data term \( \frac{1}{2\sigma_w^2}\|y - Ax\|^2 \) is already quadratic, so leave it alone. Only the prior term needs replacing.
Two properties from Chapter 8 make this legitimate. Additivity says that if you have a surrogate for each term of a sum, the sum of the surrogates is a surrogate for the sum. Composition says that if \( Q(\Delta; \Delta^\prime) \) is a surrogate for \( \rho(\Delta) \), then substituting \( \Delta = x_s - x_r \) into both gives a surrogate for \( \rho( x_s - x_r ) \) as a function of \(x\). So the whole job reduces to one scalar problem: find a quadratic surrogate for the scalar function \( \rho(\Delta) \). Solve it once and it applies to every clique in the image.
We want the surrogate to be quadratic in \( \Delta \), and we choose it symmetric about zero with no linear term:
$$ \rho( \Delta; \Delta^\prime ) = \frac{a_2}{2}\, \Delta^2 \ , \tag{9} $$where \( a_2 \) depends on the current difference \( \Delta^\prime = x_s^\prime - x_r^\prime \). This form is chosen on purpose: the surrogate prior term then looks exactly like a Gaussian prior term, which is what lets you reuse Lab 3.
There is one free number and one condition that fixes it. The surrogate must be tangent to \( \rho \) at \( \Delta = \Delta^\prime \). Matching derivatives, \( a_2 \Delta^\prime = \rho^\prime(\Delta^\prime) \), gives
$$ a_2 = \frac{ \rho^\prime( \Delta^\prime ) }{ \Delta^\prime } \ , \tag{10} $$and therefore the symmetric bound surrogate
$$ \rho( \Delta; \Delta^\prime ) = \left\{ \begin{array}{ll} \displaystyle \frac{ \rho^\prime( \Delta^\prime ) }{ 2\Delta^\prime }\, \Delta^2 & \text{if } \Delta^\prime \ne 0 \ , \\[10pt] \displaystyle \frac{ \rho^{\prime\prime}( 0 ) }{ 2 }\, \Delta^2 & \text{if } \Delta^\prime = 0 \ . \end{array} \right. \tag{11} $$The second case is the limit of the first as \( \Delta^\prime \to 0 \). This is where \( q = 2 \) earns its keep. For the QGGMRF, \( \rho^{\prime\prime}(0) \) is finite and equals \( 2 / (p\,\sigma_x^2 T^{2-p}) \). For the total variation potential \( |\Delta| \), and for the GGMRF with \( p \lt 2 \), that quantity is infinite and the whole method fails.
Because it is symmetric, this quadratic touches \( \rho \) at two points, \( \Delta = \pm\Delta^\prime \), not one. That is why it is a much tighter bound than the alternative maximum curvature surrogate, which uses the largest second derivative of \( \rho \) anywhere and is therefore too conservative. The tighter the bound, the more each iteration moves.
Now substitute (11) into the prior term of (1). Each clique contributes a weight times \( (x_s - x_r)^2 \), so collect the weights into one coefficient per clique:
$$ \tilde{b}_{s,r} \ \leftarrow \ b_{s,r}\, \frac{ \rho^\prime( x_s^\prime - x_r^\prime ) } { 2\,( x_s^\prime - x_r^\prime ) } \ . \tag{12} $$The surrogate MAP cost is then
$$ Q( x; x^\prime ) = \frac{1}{2\sigma_w^2} \| y - A x \|^2 + \sum_{\{s,r\} \in \mathcal{P}} \tilde{b}_{s,r}\, \left( x_s - x_r \right)^2 \ . \tag{13} $$Compare (13) with the Lab 3 cost. They are the same function. The only difference is that the neighbor weights are no longer fixed constants. They are recomputed from the current image at the start of each outer iteration.
For the QGGMRF, substituting (4) into (12) gives the coefficient you will implement:
$$ \tilde{b}_{s,r} \ \leftarrow \ b_{s,r}\, \frac{ |\Delta^\prime|^{p-2} }{ 2\sigma_x^p }\; \frac{ \left| \frac{\Delta^\prime}{T\sigma_x} \right|^{q-p} \left( \frac{q}{p} + \left| \frac{\Delta^\prime}{T\sigma_x} \right|^{q-p} \right) } { \left( 1 + \left| \frac{\Delta^\prime}{T\sigma_x} \right|^{q-p} \right)^2 } \ , \qquad \Delta^\prime = x_s^\prime - x_r^\prime \ , \tag{14} $$and, for \( q = 2 \), the limiting value at \( \Delta^\prime = 0 \):
$$ \tilde{b}_{s,r} \ \leftarrow \ \frac{ b_{s,r} }{ p\, \sigma_x^2\, T^{2-p} } \ . \tag{15} $$For the Gaussian prior of (2), the same formula (12) gives \( \tilde{b}_{s,r} = b_{s,r} / (2\sigma_x^2) \), a constant that never changes. The Gaussian case is the special case of this algorithm in which the surrogate is exact and the re-weighting does nothing. Your code handles both with the same outer loop.
torch.where evaluates both branches, so clamp the
denominator before dividing rather than relying on the
where to protect you.
Task 4a. Derive equation (15) from equation (14). Take the limit \( \Delta^\prime \to 0 \) with \( q = 2 \), showing the step where the two powers of \( \Delta^\prime \) cancel. This is a paper-and-pencil derivation, and it is the fastest way to understand why the code needs the threshold.
Your derivation of equation (15) from equation (14).
Task 4b. Plot \( \rho(\Delta) \) from (3) over \( \Delta \in [-0.3, 0.3] \) at \( p = 1.2 \), \( q = 2 \), \( T = 0.1 \), \( \sigma_x = 0.2 \). On the same axes plot the surrogate (11) for \( \Delta^\prime = 0.02 \) and for \( \Delta^\prime = 0.15 \), and mark the points \( \Delta = \pm\Delta^\prime \).
The plot of the potential with its two surrogates, and one sentence stating whether each surrogate touches the potential at the two marked points and lies above it everywhere else.
The outer loop re-majorizes and the inner loop runs gradient descent on the quadratic surrogate. Gradient descent is the only optimizer in this lab, as in Lab 3. The gradient of the surrogate (13) at pixel \( s \) is
$$ \left[ \nabla Q( x; x^\prime ) \right]_s = - \frac{1}{\sigma_w^2} \left[ A^{t} ( y - A x ) \right]_s + 2 \sum_{r \in \partial s} \tilde{b}_{s,r} \, ( x_s - x_r ) \ . \tag{16} $$Both terms are cheap. The first is the forward and adjoint pair you built in Lab 3. The second is a local filtering operation over the 8-point neighborhood, with weights that change each outer iteration.
Gradient descent on a quadratic is stable and monotone for a step size \( \alpha \le 1 / \lambda_{\max} \), where \( \lambda_{\max} \) is the largest eigenvalue of the Hessian of \( Q \). You can bound it without computing it. For the data term \( \| A \|^2 \le 1 \), because the blur kernel is nonnegative and sums to one. For the prior term the Hessian is twice a weighted graph Laplacian, whose eigenvalues are bounded by twice the largest row sum. So use
$$ \alpha = \left( \frac{1}{\sigma_w^2} + 4 \max_{s} \sum_{r \in \partial s} \tilde{b}_{s,r} \right)^{-1} \ . \tag{17} $$Recompute \( \alpha \) each time you recompute the \( \tilde{b}_{s,r} \). It is a conservative bound, so it is safe but not the fastest possible step.
Majorization-Minimization with Gradient Descent
Given: y, A, sigma_w, sigma_x, p, q, T, K outer, M inner
Initialize x = y
For k = 0 to K-1: # outer loop
b_tilde = surrogate weights from x, using (14) and (15)
alpha = step size from b_tilde, using (17)
For m = 0 to M-1: # inner loop
g = gradient of Q(x ; x_prime) from (16)
x = x - alpha * g
Record the TRUE cost f(x) from (1) and (3)
Note that \( x^\prime \) is the image at the top of the outer iteration. The \( \tilde{b}_{s,r} \) are held fixed for all \( M \) inner steps, so the gradient (16) is the gradient of the surrogate built at \( x^\prime \), evaluated at the current \( x \).
The neighborhood is the 8-point one from Lab 3, with \( b_{s,r} = 1/6 \) for the four horizontal and vertical neighbors and \( b_{s,r} = 1/12 \) for the four diagonal neighbors, so that \( \sum_{r \in \partial s} b_{s,r} = 1 \) at an interior pixel. As in Lab 3, a pair counts only when both pixels lie inside the image.
To keep everybody's code comparable, the neighbor differences and the surrogate weights are stored as \( (8, H, W) \) tensors, one plane per neighbor offset, in this fixed order.
NEIGHBORS = [(-1, -1), (-1, 0), (-1, 1),
( 0, -1), ( 0, 1),
( 1, -1), ( 1, 0), ( 1, 1)]
B_WEIGHTS = [1/12, 1/6, 1/12,
1/6, 1/6,
1/12, 1/6, 1/12]
def neighbor_diffs(x):
"""Return the differences x[s] - x[r] for the 8 neighbor offsets.
Plane i holds x minus x shifted by NEIGHBORS[i]. Entries whose
neighbor falls outside the image are set to zero, and the matching
entries of the mask are zero, so that a pair counts only when both
pixels are inside.
x (H, W) float64 tensor
returns (diffs, valid), each (8, H, W), float64 and bool
"""
def surrogate_weights(x, sigma_x, p, q, T, prior):
"""Return the b_tilde coefficients of equations (14) and (15).
prior "qggmrf" or "gaussian"; for "gaussian" the result is the
constant b / (2 sigma_x^2) of equation (12)
returns (8, H, W) float64 tensor, zero where the pair is invalid
"""
def step_size(b_tilde, sigma_w):
"""Return the step size alpha of equation (17)."""
def surrogate_grad(x, b_tilde, y, a, sigma_w):
"""Return the gradient of the surrogate Q, equation (16).
returns (H, W) float64 tensor
"""
def true_cost(x, y, a, sigma_w, sigma_x, p, q, T, prior):
"""Return the true cost f(x) of equation (1).
Used only for reporting. It is never called inside the optimizer.
returns float64 scalar tensor
"""
def mm_reconstruct(y, a, sigma_w, sigma_x, p, q, T,
num_outer, num_inner, prior):
"""Run the majorization-minimization algorithm above, starting at y.
returns (x, cost_history), the final image and a list of the true
cost after each outer iteration, of length num_outer + 1
whose first entry is f(y)
"""
Give your AI assistant a specification, not a description. At a
minimum it needs equations (1), (3), (4), (14), (15), (16), and (17),
the neighborhood and its weights, the signatures above, and the
pseudocode in Step 5. Require float64, no pixel loops, and the
prior argument so that the Gaussian and QGGMRF runs go
through the identical outer loop, inner loop, and cost reporting.
The specification you gave your AI assistant, quoted, and a short paragraph on what it got wrong the first time and how you found out.
Task 6a. Three checks, before you look at a single picture. Each one catches a specific bug that this lab really does produce.
torch.autograd applied to the surrogate (13), on a
small random image with \( \tilde{b} \) held fixed. Report the
largest absolute difference. This catches a wrong sign or a
missing factor of 2.prior="gaussian",
\( \sigma_x = 0.2 \), and 2000 total gradient steps, and compare
against your Lab 3 deconvolution at \( \sigma_x = 0.2 \). The two
programs use different step sizes and follow different paths, but
the cost is strictly convex, so both must arrive at the same
minimizer. Report the largest absolute pixel difference.The three numbers from Task 6a: the largest gradient error, the number of majorization trials with the smallest slack, and the largest pixel difference against your Lab 3 result.
Everybody runs the same problem, so the results can be compared.
Use kodim23.pgm scaled to \( [0,1] \) as
\( x_{\text{true}} \), and form
where \(A\) is convolution with Levin kernel 1 and circular boundaries, exactly as in Lab 3. Do not clip \(y\). Use \( \sigma_w = 0.02 \) in the cost as well, that is, assume the noise level is known.
Unless a run says otherwise, use \( K = 20 \) outer iterations and \( M = 10 \) inner gradient steps, for 200 gradient evaluations in total. The baseline prior parameters are \( p = 1.2 \), \( q = 2 \), \( T = 0.1 \), \( \sigma_x = 0.2 \). Keep \( q = 2 \) everywhere.
Pick one \( 128 \times 128 \) crop containing a sharp, high contrast edge. Use the same crop in every image figure and state its coordinates.
| Run | Prior | Parameters |
|---|---|---|
| R1 | Gaussian | \( \sigma_x = 0.2 \), 200 gradient steps |
| R2 | QGGMRF | baseline: \( p=1.2 \), \( T=0.1 \), \( \sigma_x=0.2 \) |
| R3, R4 | QGGMRF | \( \sigma_x = 0.04 \) and \( \sigma_x = 1.0 \) |
| R5, R6 | QGGMRF | \( T = 0.01 \) and \( T = 1.0 \) |
| R7, R8 | QGGMRF | \( p = 1.0 \) and \( p = 1.6 \) |
| R9, R10 | QGGMRF | baseline, but \( M=1, K=200 \) and \( M=50, K=4 \) |
R3 through R8 change one parameter at a time from the baseline. R9 and R10 hold the total number of gradient evaluations at 200 and change only how often the surrogate is rebuilt.
Four images with the same gray scale: \( x_{\text{true}} \), \( y \), the R1 reconstruction, and the R2 reconstruction. Then the same four as your \( 128 \times 128 \) crop, with the crop coordinates stated.
An edge profile plot. Pick one horizontal line crossing a strong edge in your crop, and plot the pixel values along it for \( x_{\text{true}} \), R1, and R2 on the same axes.
The crops for R3 through R8, each labeled with its parameter value.
An image of \( \sum_{r \in \partial s} \tilde{b}_{s,r} \) at the final outer iteration of R2, scaled for visibility, with the display range stated.
Two convergence plots. First, the true cost for R1 and R2 against the number of gradient evaluations, plotted as \( f(x) - f_{\min} \) on a log axis, where \( f_{\min} \) is the lowest value you reached in a long run of that configuration. State whether the true cost decreased at every outer iteration of R2, and if it ever increased, by how much and where. Second, the same plot for R2, R9, and R10.
A table with one row per run, R1 through R10, giving the final true cost, the RMSE against \( x_{\text{true}} \) in \([0,1]\) units, and the wall clock time.
Answer each in a few sentences, pointing at your own figures and numbers. These are the point of the lab.
Your answers to the eight questions above.
Submit two files through Brightspace.
lab4.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 lab4 and push it before you submit.
Concise and clear beats long.
You now have the two pieces that make modern model-based reconstruction work: a prior that does not destroy edges, and a way to optimize a non-quadratic cost using nothing but the quadratic solver you already had. In Lab 5 the prior is replaced by a denoiser, and the same idea of solving a hard problem through a sequence of easy ones shows up again.