The Random Skipping Sequential (RSS) Monte Carlo algorithm

Question

• Research paper: The Random Skipping Sequential (RSS) Monte Carlo algorithm

The main algorithm described in the given link is a Random Skipping Sequential Monte Carlo Algorithm (RSS Algorithm) designed to accelerate Markov chain Monte Carlo simulations by relaxing the strict detailed balance condition traditionally required for convergence.

Key Elements of the Algorithm are as follows:

1. Sequential Updating:

   • Each component (e.g., particle/spin) is updated sequentially using Metropolis-type acceptance rules.
   • This approach improves mobility in phase space compared to random updates.

2. Random Skipping Mechanism:

   • Alternates between standard sequential updates and updates that skip a randomly selected site.
   • This reduces oscillatory behavior and ensures irreducibility and aperiodicity, which are necessary for convergence.

3. Balance Without Strict Detailed Balance:

   • Satisfies a weaker "balance" condition instead of strict detailed balance, allowing for faster convergence while maintaining the correct equilibrium distribution.

Please verify the closeness of match of the implementation with the linked paper:

import random
import math
def metropolis_acceptance(delta_energy, temperature):
 """Calculate the Metropolis acceptance probability."""
 if delta_energy <= 0:
 return 1.0
 return math.exp(-delta_energy / temperature)
def calculate_energy_difference(grid, x, y, trial_state):
 """Calculate the energy difference for a trial move."""
 n = len(grid)
 m = len(grid[0])
 
 # Neighbors with periodic boundary conditions
 left = grid[x][y - 1 if y > 0 else m - 1]
 right = grid[x][(y + 1) % m]
 up = grid[x - 1 if x > 0 else n - 1][y]
 down = grid[(x + 1) % n][y]
 
 current_energy = grid[x][y] * (left + right + up + down)
 trial_energy = trial_state * (left + right + up + down)
 
 return trial_energy - current_energy
def generate_trial_state(current_state):
 """Generate a trial state by flipping the spin."""
 return -current_state
def update_component(grid, x, y, temperature):
 """Perform a Metropolis update on a single grid component."""
 current_state = grid[x][y]
 trial_state = generate_trial_state(current_state)
 delta_energy = calculate_energy_difference(grid, x, y, trial_state)
 if random.random() < metropolis_acceptance(delta_energy, temperature):
 grid[x][y] = trial_state # Accept the trial state
def random_skipping_sequential_update(grid, temperature, skip_probability=0.5):
 """Perform one sweep of the Random Skipping Sequential (RSS) update."""
 n = len(grid)
 m = len(grid[0])
 
 # Randomly select a site to skip with the given probability
 skip_x, skip_y = None, None
 if random.random() < skip_probability:
 skip_x = random.randint(0, n - 1)
 skip_y = random.randint(0, m - 1)
 
 # Sequentially update each site, skipping the selected one
 for x in range(n):
 for y in range(m):
 if (x, y) != (skip_x, skip_y):
 update_component(grid, x, y, temperature)
def initialize_grid(size):
 """Initialize the grid with random ±1 spins."""
 return [[random.choice([-1, 1]) for _ in range(size)] for _ in range(size)]
def calculate_total_energy(grid):
 """Calculate the total energy of the system."""
 n = len(grid)
 total_energy = 0
 for x in range(n):
 for y in range(n):
 # Neighbors with periodic boundary conditions
 left = grid[x][y - 1 if y > 0 else n - 1]
 up = grid[x - 1 if x > 0 else n - 1][y]
 total_energy -= grid[x][y] * (left + up)
 return total_energy
def print_grid(grid):
 """Print the grid state."""
 for row in grid:
 print(' '.join(f"{spin:2}" for spin in row))
 print()
def rss_simulation(grid_size, temperature, num_sweeps):
 """Simulate the 2D system using the RSS algorithm."""
 grid = initialize_grid(grid_size)
 print("Initial grid:")
 print_grid(grid)
 print("Initial energy:", calculate_total_energy(grid))
 
 for sweep in range(num_sweeps):
 random_skipping_sequential_update(grid, temperature)
 if sweep % 10 == 0:
 print(f"Sweep {sweep}: Energy = {calculate_total_energy(grid)}")
 
 print("Final grid:")
 print_grid(grid)
 print("Final energy:", calculate_total_energy(grid))
 return grid
if __name__ == "__main__":
 # Parameters
 grid_size = 10 # Size of the grid (10x10)
 temperature = 2.0 # Simulation temperature
 num_sweeps = 100 # Number of Monte Carlo sweeps
 # Run the simulation
 final_grid = rss_simulation(grid_size, temperature, num_sweeps)

• Toby Speight: You should be showing your tests that verify the correctness - Code Review may identify missing tests, but you're expected to have confirmed general correctness yourself before asking for a review.

Answer 1

Please verify the closeness of match of the implementation with the linked paper

That's your obligation.

As a first-pass review of your implementation:

The indexing is a little inside-out-and-backwards. For "left" you're modifying y, and for "up" you're modifying x; those seem swapped. You also index on [x, y] when you should index on [y, x] for row-major order. At the least this is deeply-confusing naming; at most it's a (collection of) transposition errors.

In the loop of random_skipping_sequential_update, almost all calls to update_component rely on the new values from this same update step rather than the old values. Are you absolutely sure this is what you need to do? If so, fine (though the algorithm will be stuck as slow); if not, then vectorise and update an entire matrix at a time.

You need to write unit tests, and support setting of a known random seed to make those unit tests deterministic.

Add PEP484 type hints to your signatures.

You should do at least partial vectorisation with Numpy. You probably can't vectorise the outermost "sweep" loop, though.

You mutate grid. That isn't the end of the world, but you need to call it out in comments, and especially you need to make it clear in signatures like that of random_skipping_sequential_update.

rss_simulation does too much printing. The start and end prints should be moved to a main method.

Your __main__ guard is good but insufficient. Those symbols are still in the global namespace.

[[random.choice([-1, 1]) for _ in range(size)] for _ in range(size)]

should be written as

1 - 2*rand.integers(low=0, high=2, size=(size, size), dtype=np.int32)

with rand as an np.random.Generator.

A light refactor looks like

import random
import math
def metropolis_acceptance(delta_energy: int, temperature: float) -> float:
 """Calculate the Metropolis acceptance probability."""
 if delta_energy <= 0:
 return 1.
 return math.exp(-delta_energy / temperature)
def calculate_energy_difference(
 grid: list[list[int]],
 x: int, y: int, trial_state: int,
) -> int:
 """Calculate the energy difference for a trial move."""
 n = len(grid)
 m = len(grid[0])
 # Neighbors with periodic boundary conditions
 left = grid[x][y - 1 if y > 0 else m - 1]
 right = grid[x][(y + 1) % m]
 up = grid[x - 1 if x > 0 else n - 1][y]
 down = grid[(x + 1) % n][y]
 current_energy = grid[x][y]*(left + right + up + down)
 trial_energy = trial_state*(left + right + up + down)
 return trial_energy - current_energy
def generate_trial_state(current_state: int) -> int:
 """Generate a trial state by flipping the spin."""
 return -current_state
def update_component(
 grid: list[list[int]],
 x: int, y: int, temperature: float,
) -> None:
 """Perform a Metropolis update on a single grid component."""
 current_state = grid[x][y]
 trial_state = generate_trial_state(current_state)
 delta_energy = calculate_energy_difference(grid, x, y, trial_state)
 if random.random() < metropolis_acceptance(delta_energy, temperature):
 grid[x][y] = trial_state # Accept the trial state
def random_skipping_sequential_update(
 grid: list[list[int]],
 temperature: float,
 skip_probability: float = 0.5,
) -> None:
 """Perform one sweep of the Random Skipping Sequential (RSS) update."""
 n = len(grid)
 m = len(grid[0])
 # Randomly select a site to skip with the given probability
 skip_x, skip_y = None, None
 if random.random() < skip_probability:
 skip_x = random.randint(0, n - 1)
 skip_y = random.randint(0, m - 1)
 # Sequentially update each site, skipping the selected one
 for x in range(n):
 for y in range(m):
 if (x, y) != (skip_x, skip_y):
 update_component(grid, x, y, temperature)
def initialize_grid(size: int = 10) -> list[list[int]]:
 """Initialize the grid with random ±1 spins."""
 return [[random.choice([-1, 1]) for _ in range(size)] for _ in range(size)]
def calculate_total_energy(grid: list[list[int]]) -> int:
 """Calculate the total energy of the system."""
 n = len(grid)
 total_energy = 0
 for x in range(n):
 for y in range(n):
 # Neighbors with periodic boundary conditions
 left = grid[x][y - 1 if y > 0 else n - 1]
 up = grid[x - 1 if x > 0 else n - 1][y]
 total_energy -= grid[x][y] * (left + up)
 return total_energy
def print_grid(grid: list[list[int]]) -> None:
 """Print the grid state."""
 for row in grid:
 print(' '.join(f"{spin:2}" for spin in row))
 print()
def rss_simulation(
 grid: list[list[int]], # mutated
 temperature: float, num_sweeps: int,
) -> None:
 """Simulate the 2D system using the RSS algorithm."""
 for sweep in range(num_sweeps):
 random_skipping_sequential_update(grid, temperature)
 if sweep % 10 == 0:
 print(f'Sweep {sweep}: energy={calculate_total_energy(grid)}')
def main() -> None:
 grid_size = 10 # of both x and y dimension
 temperature = 2. # of the simulation
 num_sweeps = 100 # in the Monte Carlo process
 grid = initialize_grid(grid_size)
 print('Initial grid:')
 print_grid(grid)
 print('Initial energy:', calculate_total_energy(grid))
 rss_simulation(grid, temperature, num_sweeps)
 print('Final grid:')
 print_grid(grid)
 print('Final energy:', calculate_total_energy(grid))
if __name__ == '__main__':
 main()

Initial grid:
-1 1 -1 1 1 1 1 1 -1 -1
-1 1 1 -1 1 -1 -1 1 1 1
-1 -1 -1 1 -1 -1 1 1 -1 -1
 1 -1 -1 1 -1 1 -1 1 1 -1
 1 -1 -1 -1 1 -1 -1 -1 -1 -1
 1 1 1 -1 1 -1 1 1 -1 1
 1 -1 -1 1 -1 1 1 1 1 -1
-1 1 -1 1 -1 1 1 -1 1 -1
 1 -1 -1 1 -1 1 1 1 -1 1
-1 -1 1 1 -1 1 1 -1 1 1
Initial energy: 20
Sweep 0: energy=76
Sweep 10: energy=132
Sweep 20: energy=180
Sweep 30: energy=184
Sweep 40: energy=184
Sweep 50: energy=168
Sweep 60: energy=176
Sweep 70: energy=200
Sweep 80: energy=160
Sweep 90: energy=164
Final grid:
 1 -1 1 -1 1 -1 1 -1 1 -1
-1 1 -1 1 -1 1 -1 1 -1 1
 1 -1 1 -1 1 -1 1 -1 1 -1
-1 1 -1 1 -1 1 -1 1 -1 1
 1 -1 1 -1 1 -1 1 -1 1 -1
-1 -1 -1 1 -1 1 -1 1 -1 1
 1 1 1 -1 1 -1 1 -1 1 -1
-1 1 -1 1 -1 1 -1 1 -1 1
 1 -1 1 -1 1 -1 1 -1 1 -1
-1 1 -1 -1 -1 -1 -1 1 -1 1
Final energy: 172

• STerliakov: Do you really think that code of this size (small one-file CLI script) would benefit from type hints? I'd warn against doing that, TBH, to keep it as simple as possible. Given that even a basic test suite would make a 100% coverage, net benefit from mypy would be nearly zero. This snippet is so short and simple that type hints won't improve readability either (YMMV, of course, I'm really asking for opinion)
• J_H: @STerliakov, in python and in other languages I choose to write in a restricted subset, following rules which improve readability. For example, when defining a predicate I name it is_foo() instead of foo(). Here, I feel the optional type annotations definitely improve human-to-human communication, e.g. it would never have occurred to me that delta_energy was restricted to less than float, nor that we would have a grid of strictly integer values. A docstring like """Print the grid state.""" is good, but a signature that mentions what kind of grid we accept is sooo much better.
• STerliakov: @J_H Yep, consistent naming is certainly a way to go. Not to the extent of hungarian prefix notation, but some common patterns improve readability. But regarding typing... All functions here work on the same grid, I'd really prefer to see a module-level docstring saying "grid is a 2d integer array, here's what we do with that" instead of carrying that list[list[int]] around. Should we make it a TypeAlias at least?)
• J_H: @STerliakov. Agree!
• Reinderien: Type aliases are fine. But I'm more aligned with @J_H's explanation than that of STerliakov's first comment. Script length (one way or the other) is no excuse to weaken static type analysis and self-descriptive code.

Answer 2

This is great, thank you for sharing this piece of code. I'd welcome a contribution to some project written like this.

I love how this code is organized, helper functions are small, testable and descriptive.

I agree with most of the points made in the answer by @Reinderien except two (type hints aren't mandatory and should not be treated as such; mutate/update verbs are sufficient indication that grid is mutated, I don't see any strong reason to insist on commenting that fact).

Here are a few style remarks:

• Don't grid[x][y - 1 if y > 0 else m - 1]. Python supports negative indexes, for y==0 that would become grid[x][-1], same as grid[x][m-1], so the ternary is redundant. Just grid[x][y - 1].

• You may want to combine skip_x and skip_y together - now you rebuild the tuple every time for no benefit:

  def random_skipping_sequential_update(grid, temperature, skip_probability=0.5):
   """Perform one sweep of the Random Skipping Sequential (RSS) update."""
   n = len(grid)
   m = len(grid[0])
   # Randomly select a site to skip with the given probability
   skip_coord = None
   if random.random() < skip_probability:
   skip_coord = (random.randint(0, n - 1), random.randint(0, m - 1))
   # Sequentially update each site, skipping the selected one
   for x in range(n):
   for y in range(m):
   if (x, y) != skip_coord:
   update_component(grid, x, y, temperature)
  

• Docs: please add a module-level docstring explaining in a sentence or two what this code is doing. Put a link to the article there along with its DOI (so that we can find it later if sci-hub goes extinct, shall we survive after such a loss).

• Naming: I understand the temptation to use short names. But please pick conventional ones: for a grid there's no "standard" n and m, there are height and width. x and y are fine, but may be better as row and column.