October 23, 2013

Working on Matlab integration.

Quick test shows we can get about 20 evaluations per second. Not bad, considering each evaluation consists of 9 image likleihoods. There's some room for improvement here, but probably not worth pursuing at the moment:

- leapfrog rendering: render (a) while evaluating on (b).
- cleaning up geometry and fragment shaders (fewer branches, less storage)
- try other GPU blurring routines.

Scenario: we have a decent likelihood that is linear-gaussian (and a gaussian prior), so we can compute the marginal likelihood in closed-form. However, we'd like to incorporate additional sources of evidence whose likelihoods aren't gaussian. We'll see how we can estimate the joint marginal likelihood with simple Monte-Carlo sampling (no MCMC needed, no gradient).

The old marginal likelihood looked like this:

\[
\begin{align}
p(D_1) &= \int p(x) p(D_1 | x) dx
\end{align}
\]

After introducing the extra likelihood term, the joint probability is no longer linear-gaussian, so the exact marginal likelihood involves an intractible integral. However, by re-arranging, we see we can get a good monte-carlo approximation:

\[
\begin{align}
p(D_1, D_2) &= \int p(x) p(D_1 | x) p(D_2 | x) dx \\
p(D_1, D_2) &= \int p(x | D_1) p(D_1) p(D_2 | x) dx & \left(\text{Bayes thm (see below)}\right) \\
p(D_1, D_2) &= p(D_1) \int p(x | D_1) p(D_2 | x) dx \\
p(D_1, D_2) &= p(D_1) \frac{1}{N} \sum p\left(D_2 | x^{(*)}\right) & \text{(Monte Carlo)}
\end{align}
\]

In the second line, I've replaced the first two terms using Bayes theorem. In the last line, the x-stars are samples from \(p(x | D_1)\), which we have in closed-form due to linear-Gaussian prior and likelihood for \(D_1\).

Thus, we see if at least one source of data yields a linear-gaussian likelihood, we can incorporate additional data with arbitrary likelihoods in a principled way. In many cases, \(p(x | D_1) \) has low variance, so a small number of Monte-Carlo samples are sufficient for a good estimate -- even a single sample could suffice. Even if the estimates are bad, they are unbiased, so any MCMC involving the marginal likelihood will converge to the target distribution.

We can alternatively derive this in terms of importance sampling, setting the proposal probability q(x) to \(p(x | D_1) \):

\[
\begin{align}
p(D_1, D_2) &\approx 1/N \sum p(x) p(D_1 | x)p(D_2 | x) \frac{1}{q(x)} \\
&= 1/N \sum p(x) p(D_1 | x)p(D_2 | x) \frac{p(D_1)}{p(x) p(D_1 | x)} \\
&= 1/N \sum p(D_2 | x) p(D_1) \\
&= p(D_1) 1/N \sum p(D_2 | x)
\end{align}
\]

**Implementation**

**File:** curve_tree_ml_2.m

Basic idea

- construct a thinned output index set (optional, but smart)
- construct a posterior distribution over the thinned set
- Add perturbation variance to the posterior
- take average over n trials: sample curveset and evaluate pixel likelihood

Step 2 required updating to construct_attachment_covariance, which only constructs symmetric covariance matrices. We need the covariance between indices with observations and the desired output indices. Fully refactored that function into `construct_attacment_covariance_2.m`

; confirmed correctenss in the case of the self-covariance by using an existing test for version 1 of that function. If non-symmetric, the upper-triangular blocks are processed in a second pass, swapping indices so we can re-use existing code.

Need to try view-specific sampling, i.e. sample 9 different curves from 9 different views. This refactor affects likelihood server, client, and message format. I'm worried about the performance hit, but probably not worth worrying about (or futile). Coarse sampling of indices could mitigate. In any case, view-specific sampling is probably necessary, because we're using such low blurring levels in the Bd_likelihood, so the reconstruction needs to fall near the data. We've seen how plant motion and camera miscalibration cause "good" 3D curves to reproject up to 10 pixels away from the data in some views.

**Side-note** - the tcp connection is working very reliably so far! Even after the machine sleeping/resuming several times, and suspending the server job for serveral hours, the socket is still valid and communicating flawlessly!

need to refactor send_curves.m to send `num_views`

curves instead of one.

Now receive as num_curves x num_views cell array.

Coded vector^{3} to/from message.

todo: rewrite likelihood server to receive vector3, somehow pass per-view models to likelihood

```
multi-view likelihood is out - it's one-model, multi-view. We need multi-model, multi-view.
no, MV likelihood is okay, just add an extra operator() that receives a sequence of renderables
whose size is equal to the number of views
```

Some protocol for starting and loading the likelihood server from matlab code Stress test likelihood server

Posted by Kyle Simek