# TemperingConductance.py

(Difference between revisions)
 Revision as of 16:11, 16 March 2022Werner (Talk | contribs)← Previous diff Current revisionWerner (Talk | contribs) Line 1: Line 1: + 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 # Computing the conductance of the Metropolis algorithm for the V-shaped # stationary distribution. This, for the moment, is through complete enumeration # stationary distribution. This, for the moment, is through complete enumeration Line 8: Line 11: import scipy.linalg as la import scipy.linalg as la from itertools import chain, combinations from itertools import chain, combinations - + def powerset(iterable): #powerset, but without the empty set. def powerset(iterable): #powerset, but without the empty set. s = list(iterable) s = list(iterable) return chain.from_iterable(combinations(s, r) for r in range(1, len(s) + 1)) return chain.from_iterable(combinations(s, r) for r in range(1, len(s) + 1)) - + - + def OutsideFlow(S): def OutsideFlow(S): Flow = 0.0 Flow = 0.0 Line 23: Line 26: Flow += PiStat[Table[i]] * P[i, j] Flow += PiStat[Table[i]] * P[i, j] return Weight, Flow / Weight return Weight, Flow / Weight - + n = 8 n = 8 for ReplicaChange in [k/100 for k in range(50 + 1)]: for ReplicaChange in [k/100 for k in range(50 + 1)]: Line 29: Line 32: PiStat = {} PiStat = {} Table = [] Table = [] - + for x in range(1, n + 1): for x in range(1, n + 1): Table.append((x, 0)) Table.append((x, 0))

## Current revision

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))
```