libqueso/queso — pull request #262 by dmcdougall (2014-09-04 18:07:40+00:00)
Title: Speed improvement: Use `logit` transformed transition kernel

Introduction
--

This is a pretty exciting change. I'll go slow so this is understandable.

To illustrate the change, we'll solve a simple Bayesian inverse problem in three dimensions, the seed of the chain is near the boundary of the domain and the proposal covariance is about three times the size of the domain. It's a little convoluted but this is just an illustration. For the more detail-oriented readers, I've included the example code (and input file) at the bottom of this post.

The posterior we will sample is Gaussian with mean at the mid-point of the domain.

Result
--

Here's the main timing result before this change:

```
Total run time = 236 seconds
```

And here are the samples from the first component of the chain (before the change):

![samples-before](https://cloud.githubusercontent.com/assets/328518/4153533/dbc7b038-3459-11e4-931b-14fe909192b2.png)

For comparison, here's the timing result after this change:

```
Total run time = 57 seconds
```

And here are the samples from the first component of the chain (after the change):

![samples-after](https://cloud.githubusercontent.com/assets/328518/4153648/edcfc3c8-345a-11e4-9e21-73ab44cb6191.png)

Explanation
--

So, in summary, we're getting a nice speedup with this change. Why? Well, `queso` uses a Gaussian proposal distribution, so if the problem we're solving is posed on a bounded domain, it can be the case that `queso` will propose a state outside of the domain of the problem, which is a probability zero event. In fact, as the dimension of the problem increases, the likelihood of generating a proposal outside of the domain grows (and grows very quickly).

Instead, this PR applies a coordinate transformation to the Gaussian transition kernel so that all proposals generated are within the domain of the problem. This way, we can still draw a state X from a Gaussian transition kernel, which is easy to do, and simply propose f(X), where

f(x) = (b exp(x) + a) / (1 + exp(x))

This transformation ensures that f(X) is contained in [a, b].

What is supported?
--

Of course, for some problems, the domain is not bounded. Or there might be directions in the domain that are bounded and some that are not. Or there might be directions of the domain that are only bounded on one side. This code handles all of those cases.

For directions of the domain that are unbounded, the transformation above is just the identity transformation and the transition kernel falls back to the usual Gaussian distribution. For directions that are only bounded on one side, a half-transform is done.

Code design
--

The PR is designed such that the user should not have to change their code.

The only thing I'm not certain on is where, exactly, should all the transform logic sit. Currently, there's transform logic in the PDF of the random variable f(X), which is unavoidable. There's also transform logic in the transition kernel, and this is needed because when the Markov chain wants to update the mean of the transition kernel, the transition kernel and map between coordinate spaces and set the mean of the *underlying* Gaussian distribution, not the mean of the transformed distribution.

It would be nice to just have all the transform logic sitting in place, rather than in several, but I can't think of how to do that. Maybe it is not possible.

What needs to be done?
--

I've had to update a couple of the regression test reference solutions, because different samples are drawn after this change due to the change in transition kernel.

I need to write some tests for the `InvLogit*` classes, but that's not hard to do.

Appendix
--

Example code:

```cpp
#include
#include
#include
#include
#include
#include
#include

template
class Likelihood : public QUESO::BaseScalarFunction
{
public:

Likelihood(const char * prefix, const QUESO::VectorSet & domain)
: QUESO::BaseScalarFunction(prefix, domain)
{
}

virtual ~Likelihood()
{
}

virtual double lnValue(const V & domainVector, const V * domainDirection,
V * gradVector, M * hessianMatrix, V * hessianEffect) const
{
// Just evaluate a Gaussian likelihood
double sum = 0.0;

for (unsigned int i = 0; i < domainVector.sizeLocal(); i++) {
sum += (domainVector[i] - 4.5) * (domainVector[i] - 4.5);
}

return -50.0 * sum;
}

virtual double actualValue(const V & domainVector, const V * domainDirection,
V * gradVector, M * hessianMatrix, V * hessianEffect) const
{
return std::exp(this->lnValue(domainVector, domainDirection, gradVector,
hessianMatrix, hessianEffect));
}
};

int main(int argc, char ** argv) {
MPI_Init(&argc, &argv);

QUESO::FullEnvironment env(MPI_COMM_WORLD, argv[1], "", NULL);

QUESO::VectorSpace paramSpace(env,
"param_", 3, NULL);

double min_val = 4.0;
double max_val = 5.0;

QUESO::GslVector paramMins(paramSpace.zeroVector());
paramMins.cwSet(min_val);
QUESO::GslVector paramMaxs(paramSpace.zeroVector());
paramMaxs.cwSet(max_val);

QUESO::BoxSubset paramDomain("param_",
paramSpace, paramMins, paramMaxs);

QUESO::UniformVectorRV priorRv("prior_",
paramDomain);

Likelihood lhood("llhd_", paramDomain);

QUESO::GenericVectorRV
postRv("post_", paramSpace);

QUESO::StatisticalInverseProblem
ip("", NULL, priorRv, lhood, postRv);

QUESO::GslVector paramInitials(paramSpace.zeroVector());

paramInitials[0] = 4.99;
paramInitials[1] = 4.99;
paramInitials[2] = 4.99;

QUESO::GslMatrix proposalCovMatrix(paramSpace.zeroVector());

for (unsigned int i = 0; i < 3; i++) {
proposalCovMatrix(i, i) = 3.0;
}

ip.solveWithBayesMetropolisHastings(NULL, paramInitials, &proposalCovMatrix);

MPI_Finalize();

return 0;
}
```

Example input file:

```
###############################################
# UQ Environment
###############################################
#env_help = anything
env_numSubEnvironments = 1
env_subDisplayFileName = outputData/display
env_subDisplayAllowAll = 0
env_subDisplayAllowedSet = 0
env_displayVerbosity = 2
env_syncVerbosity = 0
env_seed = 0

###############################################
# Statistical inverse problem (ip)
###############################################
#ip_help = anything
ip_computeSolution = 1
ip_dataOutputFileName = outputData/sipOutput
ip_dataOutputAllowedSet = 0

###############################################
# 'ip_': information for Metropolis-Hastings algorithm
###############################################
#ip_mh_help = anything
ip_mh_dataOutputFileName = outputData/sipOutput
ip_mh_dataOutputAllowedSet = 0 1

ip_mh_rawChain_dataInputFileName = .
ip_mh_rawChain_size = 1000000
ip_mh_rawChain_generateExtra = 0
ip_mh_rawChain_displayPeriod = 50000
ip_mh_rawChain_measureRunTimes = 1
ip_mh_rawChain_dataOutputFileName = outputData/ip_raw_chain
ip_mh_rawChain_dataOutputAllowedSet = 0 1
ip_mh_rawChain_computeStats = 1

ip_mh_displayCandidates = 0
ip_mh_putOutOfBoundsInChain = 0
ip_mh_tk_useLocalHessian = 0
ip_mh_tk_useNewtonComponent = 1
ip_mh_dr_maxNumExtraStages = 1
ip_mh_dr_listOfScalesForExtraStages = 5.
ip_mh_am_initialNonAdaptInterval = 0
ip_mh_am_adaptInterval = 10
ip_mh_am_eta = 1.92
ip_mh_am_epsilon = 1.e-5

ip_mh_filteredChain_generate = 0
```