Estimate the number of self-avoiding walks of length $n$

An answer to this question on the Scientific Computing Stack Exchange.

Question

For the past couple of days, I have been thinking a lot and searching online for an algorithm capable of estimating the number of Self-Avoiding Walks (SAW) of length $n$ in $\mathbb{Z}^2$. There is a very simple way to do so by simply generating all random walks of length $n$ and then checking whether each one is self-avoiding or not. This works fine for very small values of $n$ but as soon as $n\ge 10$, this method becomes completely obsolete because of the exponential complexity. I was capable of creating a function in Python capable of generating a self-avoiding path of length $n$ using the following algorithm:

$1)$ start from an arbitrary SAW of length $n$ (it's easy to just pick the walk of length $n$ that simply goes up the $y$-axis for example).

$2)$ select a random point in the SAW other than its extremities and let this point be the pivot, and then pick an angle of rotation from $\pi$, $\frac\pi2$ and $-\frac\pi2$.

$3)$ for every point placed strictly after the pivot in the SAW, rotate it around the pivot with the selected angle of rotation.

$4)$ if the newly obtained random walk is self-avoiding then we keep it, otherwise repeat steps $2)$ and $3)$ until you find a pair that works. You can repeat these steps as many times as you wish to generate a seemingly unrelated SAW to the one we started with.

The problem is this method gives one SAW of length $n$ and doesn't exactly help to estimate the number of SAWs of length $n$. What algorithm can one thus use to estimate the number of SAWs of length $n$ and how could it be implemented in Python?

Answer

Background

The number of self-avoiding walks of length N on a square lattice is:

0   1 
1   4 
2   12 
3   36 
4   100 
5   284 
6   780 
7   2172 
8   5916 
9   16268 
10   44100 
11   120292 
12   324932 
13   881500 
14   2374444 
15   6416596 
16   17245332 
17   46466676 
18   124658732 
19   335116620 
20   897697164 
21   2408806028 
22   6444560484 
23   17266613812 
24   46146397316 
25   123481354908 
26   329712786220 
27   881317491628 
28   2351378582244 
29   6279396229332 
30   16741957935348 
31   44673816630956 
32   119034997913020 
33   317406598267076 
34   845279074648708 
35   2252534077759844 
36   5995740499124412 
37   15968852281708724 
38   42486750758210044 
39   113101676587853932 
40   300798249248474268 
41   800381032599158340 
42   2127870238872271828 
43   5659667057165209612 
44   15041631638016155884 
45   39992704986620915140 
46   106255762193816523332 
47   282417882500511560972 
48   750139547395987948108 
49   1993185460468062845836 
50   5292794668724837206644 
51   14059415980606050644844 
52   37325046962536847970116 
53   99121668912462180162908 
54   263090298246050489804708 
55   698501700277581954674604 
56   1853589151789474253830500 
57   4920146075313000860596140 
58   13053884641516572778155044 
59   34642792634590824499672196 
60   91895836025056214634047716 
61   243828023293849420839513468 
62   646684752476890688940276172 
63   1715538780705298093042635884 
64   4549252727304405545665901684 
65   12066271136346725726547810652 
66   31992427160420423715150496804 
67   84841788997462209800131419244 
68   224916973773967421352838735684 
69   596373847126147985434982575724 
70   1580784678250571882017480243636 
71   4190893020903935054619120005916

This is drawn from Iwan Jensen's 2013 paper A new transfer-matrix algorithm for exact enumerations: self-avoiding walks on the square lattice.

The number of SAW walks is believed to be approximated by the equation $$c_n \sim Aµ^n n^{γ−1}$$ where $\gamma=43/32, µ\sim2.638, A\sim1.17704242$. These values are estimated in Jensen (2013). The value of γ has been known since at least Nienhuis (1984, "Critical behavior of two-dimensional spin models and charge asymmetry in the Coulomb gas").

The integer sequence above is also known as A001411 in the Online Encyclopedia of Integer Sequences.

The Jensen (2013) algorithm used ~16,500 CPU hours to enumerate self-avoiding walks up to length 79. Zbarsky (2019) suggests an asymptotically faster algorithm, but is not sure if it would be faster than the Jensen (2013) algorithm for such small values of N.

Nathan Clisby's 2021 paper Enumerative Combinatorics of Lattice Polymers doesn't deal with 2D walks specifically, but seems like a good overview of the broader uses of self-avoiding walks and adjacent areas.

Algorithms

Your algorithm generates paths randomly. This makes it poorly suited to enumeration because:

  1. It may generate the same path multiple times.
  2. Due to this, it will take it a long time to explore the entirety of the space.

A good algorithm should exhaustively explore the space by visiting each path only once.

The easiest way to achieve this is with recursive backtracking. I've implemented an algorithm to do so below:

from typing import Final, List, Set, Tuple
def walk(n_max: int, n: int, x: int, y: int, used: Set[Tuple[int, int]]) -> int:
  """
  Recursive function to calculate the number of self-avoiding walks
  Our strategy here is to build a function that generates self-avoiding walks in
  a depth-first manner using backtracking to avoid duplication and a hash-set to
  avoid overlaps.
  Args:
    n_max   - Length of the walk
    n       - How far we are into the walk
    x       - Current x location
    y       - Current y location
  Returns:
    The number of self-avoiding walks of length `n_max`
  """
  # Possible displacement values on a square lattice
  dr_vals: Final[List[Tuple[int, int]]] = [(-1, 0), (1, 0), (0, 1), (0, -1)]
  # End recursion. We have discovered 1 walk!
  if n==n_max:
    return 1
  # Make a note that we're now at location (x,y) so it is unavailable
  used.add((x,y))
  # Body of recursion: we will count the walks leading away from (x,y)
  walks = 0
  # Look at all the possible directions we can go from here
  for dr in dr_vals:
    # Determine next coordinates given the displacement
    nx = x + dr[0]
    ny = y + dr[1]
    # Don't visit the neighbouring square twice!
    if (nx, ny) in used:
      continue
    # Visit the neighbouring square and count how many walks it has
    walks += walk(n_max, n + 1, nx, ny, used)
  # We've finished the body of the recursion, we now go back up the stack
  # Since we're going back up the stack we no longer occupy (x,y)
  used.remove((x,y))
  return walks
def main():
  # Make a table of how many walks of length n there are
  # Given the exponential increase the cost of computation and simplicity of this
  # algorithm, this will probably run forever
  for n in range(30):
    used: Set[Tuple[int, int]] = set([(1000, 1000)])
    print(n, walk(n, 0, 0, 0, used))
main()