Multi-image processing with PyMC3

Question

I have an image processing problem I thought I could use to experiment with learning more about PyMC3. I have spent a good amount of time fiddling with non-linear solvers and brute-force methods, and so far nothing makes me happy.

The problem I have involves a complicated method for co-registering two images of the same scene but recorded in different modalities. Think of trying to match a regular black and white visible image with a thermal infrared image. Or, from a medical-imaging perspective, think of trying to match MRI data with X-ray data.

Just to keep things simple, I can represent my data processing work flow with the following function:

def process_and_compare(image_src, image_dst, parameters):
    """
    Parameters
    ----------
    image_src : 2D array
    image_trg : 2D array
    parameters : sequence of 7 scalars defining image transform

    Output
    ------
    metric : Scalar value indicating how well the transformed source image
             matches up with the target image.
    """

    image_src_warp = image_warper(image_src, parameters)
    metric = compare_two_images(image_src_warp, image_trg)

    return metric

This function takes as input two images and a vector of model parameters. Some complicated number crunching happens on the inside. When it's done a scalar is returned that indicates how well the model (defined solely by the parameter vector) aligns the two images. The details of how the source image is warped or how two images are compared are black boxes for now. In the end I ultimately the primary result I want is the warped image corresponding to the model resulting in the best match. But right now while I'm still playing with my algorithms, I think I can learn a lot by visualizing the posterior distributions of my model parameters for some simple test case images. I initially thought PyMC would make this easy, but once I started looking into the actual implementation detail I got a bit confused.

I have looked over recent PyMC3 presentations by Thomas Wiecki and also read through much of the great online book by Cam Davidson-Pilon. So far it looks to me that the great new features of PyMC3 versus PyMC2 are (in part) a snazzy model specification syntax and automatic use of Theano for processing acceleration.

In the examples I've seen so far, it looks like the data model is now often fully specified using the new syntax system. But in my case I have this more complicated function.

Here are my questions:

1. Can someone point me to an existing PyMC example involving a black box data model implemented as a user function? Either PyMC2 or PyMC3 would great!

2. Will I be able to realize the benefits of Theano deep in my Python data model function once I figure out how to make this work with PyMC3?

Answer 1

The function that you specified above would be incorporated into a PyMC model as a Deterministic node, where it is calculated based on some (presumably) stochastic parent nodes (your parameters). This node would then be connected downstream to a likelihood (observed stochastic node) that provides the information for fitting the parameters. For example, you might have some parametric distribution that describes the distribution of the error corresponding to the metric output by process_and_compare.

The PyMC wiki has several model examples from a range of domains for PyMC 2. There are PyMC 3 examples in the pymc/examples folder in the master branch.

As far as Theano goes, the motivation behind using it as a dependency for PyMC is due the fact that the current state-of-the-art in MCMC involves using gradient information, so we needed the ability to calculate gradients for arbitrary models. We hope to eventually benefit from its GPU capabilities, but for now its just for the gradients. All PyMC's objects are Theano tensors in version 3, so if you have other plans for Theano in the context of building Bayesian models, then chances are it can be made to work. For example, we might eventually want to implement probabilistic graphical models in PyMC, so Theano will probably facilitate that as well.