# TemperingConductance.py

### From Werner KRAUTH

This is the algorithm TemperingConductance.py, which computes the conductance of the Hildebrand problem.It was presented in the 8th lecture of the 2022 course "Advanced topics in MCMC"

# Computing the conductance of the Metropolis algorithm for the V-shaped # stationary distribution. This, for the moment, is through complete enumeration # of the power set of the set of all samples. A rigorous argument can be made # also. # import random import numpy as np import scipy.linalg as la from itertools import chain, combinations def powerset(iterable): #powerset, but without the empty set. s = list(iterable) return chain.from_iterable(combinations(s, r) for r in range(1, len(s) + 1)) def OutsideFlow(S): Flow = 0.0 Weight = 0.0 for i in S: Weight += PiStat[Table[i]] for j in range(2 * n): if j not in S: Flow += PiStat[Table[i]] * P[i, j] return Weight, Flow / Weight n = 8 for ReplicaChange in [k/100 for k in range(50 + 1)]: const = 4.0 / n ** 2 PiStat = {} Table = [] for x in range(1, n + 1): Table.append((x, 0)) Table.append((x, 1)) # # factor of 1/2 because the total must be normalized # PiStat[(x, 0)] = 1.0 / float(n) / 2.0 PiStat[(x, 1)] = const * abs( (n + 1) / 2 - x) / 2.0 PiStat[(0, 0)] = 0.0 PiStat[(0, 1)] = 0.0 PiStat[(n + 1, 0)] = 0.0 PiStat[(n + 1, 1)] = 0.0 Position = (1, 0) PTrans = np.eye(2 * n) for x in range(1, n + 1): for Rep in [0, 1]: i = Table.index((x, Rep)) for Dir in [-1, 1]: if PiStat[(x + Dir, Rep)] > 0.0: j = Table.index((x + Dir, Rep)) PTrans[i, j] = min(1.0, PiStat[(x + Dir, Rep)] / PiStat[(x, Rep)]) / 2.0 PTrans[i, i] -= PTrans[i, j] PReplica = np.zeros((2 * n,2 * n)) for x in range(1, n + 1): i = Table.index((x,0)) j = Table.index((x,1)) PReplica[i, j] = ReplicaChange * min(1.0, PiStat[(x, 1)] / PiStat[(x, 0)]) PReplica[i, i] = 1.0 - PReplica[i, j] PReplica[j, i] = ReplicaChange * min(1.0, PiStat[(x, 0)] / PiStat[(x, 1)]) PReplica[j, j] = 1.0 - PReplica[j, i] P = PTrans @ PReplica # # compute the conductance. # x = powerset([i for i in range(2 * n)]) MinFlow = 100.00 for i in x: S = set(list(i)) Weight, Flow = OutsideFlow(S) if Weight <= 0.5: if Flow < MinFlow: MinFlow = Flow MinS = S print(MinS, MinFlow, '2 / 2n + 1 / n^2 = ', 1.0 / (2.0 * n) + 1.0 / (n ** 2))