Would an "importance Gibbs" sampling method work?

Question

I suspect this is a fairly unusual and exploratory question, so please bear with me.

I am wondering if one could apply the idea of importance sampling to Gibbs sampling. Here's what I mean: in Gibbs sampling, we change the value of one variable (or block of variables) at a time, sampling from the conditional probability given the remaining variables.

However, it may not be possible or easy to sample from the exact conditional probability. So instead we sample from a proposal distribution $q$ and use, for example, Metropolis-Hastings (MH).

So far, so good. But here is a divergent path: what happens if we, instead of using MH, use the same idea used in importance sampling, namely we sample from $q$ and keep an importance weight $p/q$ of the current sample?

In more detail: assume we have variables $x_1,\dots,x_n$ and a factored distribution $\phi_1,\dots,\phi_m$ so that $p \propto \prod_{i=1}^m \phi_i$. We keep the proposal probability $q_i$ used to sample the current value of each variable $x_i$. At each step we change a subset of the variables and update $p(x)/q(x)$ (only the factors of $p$ and $q$ that are affected). We take the samples and their importance weight to compute whatever statistic we are interested in.

Would this algorithm be correct? If not, any clear reasons why not? Intuitively it makes sense to me as it seems to be doing the same thing importance sampling does, but with dependent samples instead.

I did implement this for a Gaussian random walk model and observed that the weights become smaller and smaller (but not monotonically), so the initial samples end up having too much importance and dominate the statistic. I'm pretty certain the implementation is not buggy, because at each step I compare the updated weight to an explicit brute-force computation of it. Note that the weights do not go down indefinitely to zero, because they are $p/q$ where both $p$ and $q$ are products of a finite number of densities, and each sample is obtained from a Normal distribution that only rarely will be zero.

So I am trying to understand why the weights go down like that, and whether this is a consequence of this method being actually not correct.

---

Here's a more precise definition of the algorithm, applied to a Gaussian random walk on variables $X_1,\dots,X_n$. The code follows below.

The model is simply $X_i \sim \mathcal N(X_{i-1}, \sigma^2), i = 1,\dots,n$, with $X_0$ fixed to $0$.

The weight of the current sample is $\frac{\prod_i p(x_i)}{\prod_i q(x_i)}$, where $p$ are the Gaussian densities and $q$ are the distributions from which the current values have been sampled. Initially, we simply sample the values in a forward manner, so $q = p$ and the initial weight is $1$.

Then at each step I choose $j \in \{1,\dots,n\}$ to change. I sample a new value $x'_j$ for $X_j$ from $\mathcal N(X_{j-1},\sigma^2)$, so this density becomes the new used proposal distribution for $X_j$.

To update the weight, I divide it by the densities $p(x_j | x_{j-1})$ and $p(x_{j+1} | x_j)$ of old value $x_j$ according to $x_{j-1}$ and $x_{j+1}$, and multiply by the densities $p(x'_j | x_{j-1})$ and $p(x_{j+1} | x'_j)$ of new value $x'_j$ according to $x_{j-1}$ and $x_{j+1}$. This updates the numerator $p$ of the weight.

To update the denominator $q$, I multiply the weight by the old proposal $q(x_j)$ (thus removing it from the denominator) and divide it by $q(x'_j)$.

(Because I sample $x'_j$ from the normal centered on $x_{j-1}$, $q(x'_j)$ is always equal to $p(x'_j | x_{j-1})$ so they cancel out and the implementation does not actually use them).

Like I mentioned before, in the code I compare this incremental weight computation to the actual explicit computation just to be sure.

---

Here's the code for reference.

println("Original sample: " + currentSample);
int flippedVariablesIndex = 1 + getRandom().nextInt(getVariables().size() - 1);
println("Flipping: " + flippedVariablesIndex);
double oldValue = getValue(currentSample, flippedVariablesIndex);
NormalDistribution normalFromBack = getNormalDistribution(getValue(currentSample, flippedVariablesIndex - 1));
double previousP = normalFromBack.density(oldValue);
double newValue = normalFromBack.sample();
currentSample.set(getVariable(flippedVariablesIndex), newValue);
double previousQ = fromVariableToQ.get(getVariable(flippedVariablesIndex));
fromVariableToQ.put(getVariable(flippedVariablesIndex), normalFromBack.density(newValue));
if (flippedVariablesIndex < length - 1) {
    NormalDistribution normal = getNormalDistribution(getValue(currentSample, flippedVariablesIndex + 1));
    double oldForwardPotential = normal.density(oldValue);
    double newForwardPotential = normal.density(newValue);
    // println("Removing old forward potential " + oldForwardPotential);
    currentSample.removePotential(new DoublePotential(oldForwardPotential));
    // println("Multiplying new forward potential " + newForwardPotential);
    currentSample.updatePotential(new DoublePotential(newForwardPotential));
}

// println("Removing old backward potential " + previousP);
currentSample.removePotential(new DoublePotential(previousP));
// println("Multiplying (removing from divisor) old q " + previousQ);
currentSample.updatePotential(new DoublePotential(previousQ));

println("Final sample: " + currentSample);
println();

// check by comparison to brute force calculation of weight:
double productOfPs = 1.0;
for (int i = 1; i != length; i++) {
    productOfPs *= getNormalDistribution(getValue(currentSample, i - 1)).density(getValue(currentSample, i));
}
double productOfQs = Util.fold(fromVariableToQ.values(), (p1, p2) -> p1*p2, 1.0);
double weight = productOfPs/productOfQs;
if (Math.abs(weight - currentSample.getPotential().doubleValue()) > 0.0000001) {
    println("Error in weight calculation");
    System.exit(0);
}

Answer 1

This is an interesting idea, but I see several difficulties with it:

1. contrary to standard importance sampling, or even Metropolised importance sampling the proposal is not acting in the same space as the target distribution, but in a space of smaller dimension so validation is unclear [and may impose to keep weights across iterations, hence facing degeneracy]
2. the missing normalising constants in the full conditionals change at each iteration but are not accounted for [see below]
3. the weights are not bounded, in that along iterations, there will eventually be simulations with a very large weight, unless one keeps track of the last occurrence of an update for the same index $j$, which may clash with the Markovian validation of the Gibbs sampler. Running a modest experiment with $n=2$ and $T=10^3$ iterations shows a range of weights from 7.656397e-07 to 3.699364e+04.

To get into more details, consider a two-dimensional target $p(\cdot,\cdot)$, including the proper normalising constant, and implement the importance Gibbs sampler with proposals $q_X(\cdot|y)$ and $q_Y(\cdot|x)$. Correct importance weights [in the sense of producing the correct expectation, i.e., an unbiased estimator, for an arbitrary function of $(X,Y)$] for successive simulations are either $$\dfrac{p(x_t,y_{t-1})}{q_X(x_t|y_{t-1})m_Y(y_{t-1})}\qquad\text{or}\qquad\dfrac{p(x_{t-1},y_{t})}{q_Y(y_t|x_{t-1})m_X(x_{t-1})}$$ where $m_X(\cdots)$ and $m_Y(\cdot)$ are the marginals of $p(\cdot,\cdot)$. Or equivalently $$\dfrac{p_X(x_t|y_{t-1})}{q_X(x_t|y_{t-1})}\qquad\text{or}\qquad\dfrac{p_Y(y_{t}|x_{t-1})}{q_Y(y_t|x_{t-1})}$$ In either case, this requires the [intractable] marginal densities of $X$ and $Y$ under the target $p(\cdot,\cdot)$.

It is worthwhile to compare what happens here with the parallel importance weighted Metropolis algorithm. (See for instance Schuster und Klebanov, 2018.) If the target is again $p(\cdot,\cdot)$ and the proposal is $q(\cdot,\cdot|x,y)$, the importance weight $$\dfrac{p(x',y')}{q(x',y'|x,y)}$$is correct [towards producing an unbiased estimate] and does not update the earlier weight but starts from scratch at each iteration.

(C.) A correction to the original importance Gibbs proposal is to propose a new value for the entire vector, e.g., $(x,y)$, from the Gibbs proposal $q_X(x_t|y_{t-1})q_Y(y_t|x_{t})$, because then the importance weight $$\dfrac{p(x_t,y_t)}{q_X(x_t|y_{t-1})q_Y(y_t|x_{t})}$$ is correct [missing a possible normalising constant that is now truly constant and does not carry from previous Gibbs iterations].

A final note: for the random walk target considered in the code, direct simulation is feasible by cascading: simulate $X_1$, then $X_2$ given $X_1$, &tc.