logTransitionProbability in ICP-proposal

[Paul Dickens]

Hi,

I was wondering if you could help me understand the logic in this code snippet. The function in question (logTransitionProbability) computes the probability of transitioning from one set of parameters (`from`) to another (`to`). As far as I can work out, firstly, the blending effect of the step length is removed, and then the logpdf of the shape parameters is computed.

What I can't understand is that the posterior distribution is calculated in line 79, and then used  to calculate the logpdf in line 83. However, the logpdf is simply the likelihood of drawing the coefficients from a multivariate normal distribution with mean zero and identity matrix covariance. This is equal to the likelihood of drawing the coefficients from the prior distribution, not the posterior.

In your accompanying paper, the transition probability is defined as "equal to the probability of sampling alpha0 from the analytically computed posterior model

computed in step 4 of the ICP-proposal."

It seems to me that the coefficients should be adjusted to account for the different mean and pcaBasis in the posterior distribution before the logpdf is calculated.

I would be grateful if you could confirm whether I'm correct in this or help correct any misunderstanding.

Kind regards,

Paul Dickens

[Dennis Madsen]

Hi Paul, 

Thanks a lot for your interest in our recent research. 

You are indeed right that there is a problem in the current public version. The parameters need to be projected into the space of the posterior model as well. 

Below is an example of a correct solution. I'll have the repository updated soon. Thanks for pointing this out. 

With the correct implementation, the femur experiment as shown in the paper will still give the same acceptance rate and "best-fit" performance. 

override def logTransitionProbability(from: ModelFittingParameters, to: ModelFittingParameters): Double = {
val pos = cashedPosterior(from)
val posterior = StatisticalMeshModel(referenceMesh, pos)

val compensatedTo = from.shapeParameters.parameters + ((to.shapeParameters.parameters - from.shapeParameters.parameters) / stepLength)
val toMesh = model.instance(compensatedTo)

val projectedTo = posterior.coefficients(toMesh)
val pdf = pos.logpdf(projectedTo)
pdf
}

Best,

Dennis

[Paul Dickens]

Hi Dennis,

Thanks for your reply and code. And I appreciate you making the MCMC IPC-based registration method available to the public. It's really great work. I'll edit my local copy with your modifications. Just out of interest, I'm assuming that this proposal generator is sampling from a symmetric distribution, so the log transition ratio (that is used in the Metropolis Hastings algorithm) will always be zero, so won't affect the outcome. Is this a correct assumption?

Also, I've found what I think are a couple of typos in the code at https://github.com/unibas-gravis/icp-proposal. If you don't mind I'll submit these as PRs for your consideration.

All the best,

Paul

[Dennis Madsen]

Hi Paul,

For the MH algorithm, the proposals do not need to be symmetric, this is only for the Metropolis algorithm. 

The log-transition probability is exactly there to adjust for the asymmetry in the proposal. So the ratio for the proposal will not always be zero. 

Contributions to our projects are always very welcome - so feel free to submit a PR for the project. 

Best,

Dennis

[Paul Dickens]

Hi Dennis,

Reviving an old thread about your icp-proposal method. I have implemented the code you posted above for logTransitionProbability in the IcpNonRigidProposal class, but I am finding that the transition ratio for the proposal tends to be very large (in other words it is is much more likely to transition in the forward direction than in the backward direction. This stands to reason, I guess, because of the way the ICP-based posterior distribution works. However, since the Metropolis-Hastings algorithm calculates the acceptance criterion as (proposedP - currentP - t), where t is the transition probability ratio, all log values), this calculation is getting swamped by the very high t, so (almost) all proposals end up being rejected. It also tends to scale with the rank of the Gaussian process, which seems strange (so a model with rank 100 almost always leads to the proposal being rejected, while a model with rank 10 ends up with some accepted, some rejected).

I can see the logic in your method, but am wondering whether the line

val compensatedTo = from.shapeParameters.parameters + ((to.shapeParameters.parameters - from.shapeParameters.parameters) / stepLength)

could be part of the problem. If I remove this compensation I get much better convergence. But of course, this is likely not mathematically valid.

I could be making an error somewhere else, but it would be great to know if you have tested this change and are confident in it. I see there is still an open issue (https://github.com/unibas-gravis/icp-proposal/issues/2) for this. Hope this question makes sense.

Kind regards,

Paul

[Paul Dickens]

Just to put some numbers to the issue above, I tested the change (in https://github.com/unibas-gravis/icp-proposal/issues/2) on the BfmFitting project in the icp-proposal repo this morning. The logTransitionProbability in the forward direction is always around -75, but in the backwards direction it's around -1050. (These change around a bit as the sampling process is random, but those numbers are typical. The difference (fw - bw) of around 975 is enough to swamp the MH algorithm and virtually ensure that the proposals are never accepted.

[Dennis Madsen]

Hi Paul, 

Thanks a lot for the feedback and your continuous interest in the project. 

I have recently undergone a large number of new experiments using the correct transition probability as proposed in the issue on GitHub. 

With the new experiments, I am confident that the method works and you are able to get a good acceptance rate. But with the change, the hyper parameters does need an update - which might be the problem you are facing. I especially found it difficult for the face experiment due to the likelihood that is used (collective average + Hausdorff). 

The setting where the method works very well is if the likelihood function can be evaluated over a constant number of points, i.e. like the independent point evaluator as used in the femur experiments. 

If the target is a subset of the model, you can then set the evaluation direction to targetToReference. If the model is a subset of the target, you can set it to referenceToTarget. When this is not the case, I would instead stick to the pure “collective average” and remove the Hausdorff evaluation. It is possible to set the hyper parameters for a good acceptance rate, but I found that I had to do this on target basis. 

I am currently working on a way to have the hyper-parameters semi-automatically computed, and build this into a software package that would allow for easier adoption to new domains. So hopefully more to share on that soon. I will also have the current repo updated with usable experiments using the correct transition probability in the coming months. 

Do you have a specific domain you are testing it on? Or purely trying to reproduce the experiments from the paper? 

Best

Dennis

-- 
You received this message because you are subscribed to a topic in the Google Groups "scalismo" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/scalismo/FytrdvlKgk4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to scalismo+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/scalismo/91d6c058-8353-48a0-a0ff-e0a38cc4ebecn%40googlegroups.com.

[Paul Dickens]

Hi Dennis,

Thanks for your detailed reply and advice. I'm testing the code on a head model with face samples, so very similar to the BFM experiment. The target is a subset of the model, so I've changed the evaluator to targetToReference as you suggest. Do you recommend continuing to use  collective average + Hausdorff in this case? Also, I assume that the proposal generator should use the same method (target to reference ICP)?

I'll look out for the update to the repo with experiments, as well as the auto-updating hyperparameters. That sounds like a very worthwhile improvement. Using your updated transition probability calculation, I've found that the acceptance ratio improves if I increase the stepLength from 0.1 to 0.5. I assume stepLength is one of the hyperparameters you mention. Are there others that are worth experimenting with?

Thanks again,

Paul

[Paul Dickens]

I have one other query about the Scalismo code. I am using a MixtureProposal to combine shape updates with translation and rotation updates, as in https://scalismo.org/docs/tutorials/tutorial15, but using the ICP-proposal method for the shape updates. The logTransitionRatio always ends up zero, even though the transition ratio is very large when using the shape update proposal alone. I would appreciate some help understanding the following code (from sampling/proposals/MixtureProposal.scala):

  override def logTransitionProbability(from: A, to: A): Double = {

    val transitions = generators.map(g => g.logTransitionProbability(from, to))

    if (transitions.exists(_.isNaN))

      throw new Exception("NaN transition probability encountered!")

    if (transitions.exists(!_.isInfinite)) {

      val maxExpo = transitions.max

      val sum = mixtureFactors.zip(transitions).map { case (f, t) => f * math.exp(t - maxExpo) }.sum

      val fwd = math.log(sum) + maxExpo

      fwd

    } else

      Double.NegativeInfinity

  }

Using this code, if a shape update is proposed the shape log transition probability is large and negative, while the translation and rotation transition probabilities are very high (since they are unchanged). The result of taking the weighted average after exponentiation is that the rotation and translation terms dominate. This leads to a log transition ratio of zero when the fw - bw calculation is performed.

Would it make more sense to override the logTransitionRatio in Mixture Proposal directly? Maybe something like this:

  /** transition ratio forward probability / backward probability */

  override def logTransitionRatio(from: A, to: A): Double = {

    val sum = generators.map { case g => g.logTransitionRatio(from, to) }.sum

    sum

  }

Using a method like this reduces to the same log transition ratio as we would get using the shape update proposal alone. I may be misunderstood something, but it seems strange for the log transition ratio always to be zero using a MixtureProposal.

Thanks,

Paul

[Dennis Madsen]

Hi Paul, 

I would recommend using the independentPointEvaluator likelihood if possible for your application. I find that the parameters are much easier to adjust and often with a higher acceptance rate. 

As such the proposal and evaluator can be configured separately, so there is not a strict criteria that you e.g. have to use target to model on both. 

For the mixture proposal transitionratio, it does indeed seem like a term is missing in the tutorials if the parameters stay the same. 

I've just had a short talk with some of the other scalismo developers. We'll have to look a bit more into it and then get back to you. Hopefully by next week. 

Best

Dennis

[Paul Dickens]

Hi Dennis,

Thanks for getting back to me. I'd love to hear what you and the other developers think about the logTransitionProbability in the MixtureProposal. At the moment, I'm just using the logTransitionProbability of the last active generator, weighted by the mixture factor for that generator. To me this makes more sense than adding together the unnormalized probabilities for all the different proposal generators. But I'm not sure if this approach is mathematically correct.

Cheers,

Paul

[Dennis Madsen]

Hi Paul, 

We've been looking more into what happens with the proposal mixture. 
For now, a pull-request is pending with some fixes to the proposals, for when they are used in a mixture: https://github.com/unibas-gravis/icp-proposal/pull/3 
One part is returning an appropriate transition probability in case the proposal is unable to transition.

Example for the ICP-proposal:

if (to.copy(shapeParameters = from.shapeParameters).allParameters != from.allParameters) {
Double.NegativeInfinity
}
else {

 // Compute transition probability

}

In case the transition involved something else than shape (rotation, translation), then a negative infinite can be returned as the proposal will never be able to do this transition. 

In the paper, we always only used the same proposal, so this was never an issue. The tutorials will also be updated in the near future with this change. 

Another part is the adjustment for the asymmetry in the ICP-proposal. This can be done with a mixture of ICP-proposals with different step-size.

Alternatively, we've found that even a symmetrical random pose proposal can be helpful. Only using a single ICP-proposal will increase the speed and the mixture with a wide random pose increases the acceptance rate (I get +50%). 

In the pull-request I've also divided the face fitting into 2 files: BfmFittingComplete (uses the independent point evaluator to register a complete face), BfmFittingPartial (uses the Hausdorff and collective likelihood to register a face with a missing nose). 

As for the unnormalized probabilities. As long as proposals operate orthogonal to each other, their unnormalized values can be used. If not, the correct transition probability is necessary.

As mentioned earlier, we are currently working on having many of the parameters automatically adjusted. For now, it is however unsure if it will be in the same repository, or in a new repository. 

Best

Dennis