Local search in Julia

Julia's array syntax and broadcasting fit physics problems like a glove. Local search converges fast and reads cleanly.

1. 01
   Represent the graph as a weighted adjacency matrix. Symmetric, zero diagonal — and Julia's BitArray is the natural fit for the partition assignment.

   using Random
   
   function build_graph(n::Int, p::Float64; seed=42)
       rng = MersenneTwister(seed)
       A = zeros(Float64, n, n)
       for i in 1:n, j in i+1:n
           if rand(rng) < p
               A[i, j] = A[j, i] = 1.0
           end
       end
       A
   end

2. 02
   Cut value: count edge weights where the two endpoints are on opposite sides. Vectorize over the upper triangle to keep it tight.

   function cut_value(A::Matrix{Float64}, side::BitVector)
       s = 0.0
       n = size(A, 1)
       for i in 1:n, j in i+1:n
           if side[i] != side[j]
               s += A[i, j]
           end
       end
       s
   end

3. 03
   Greedy local search: starting from a random partition, flip any vertex whose flip increases the cut, until no such flip exists. Polynomial-time, lands at a local optimum within a factor of 2 of the global maximum.

   function local_search(A::Matrix{Float64}; restarts=20)
       n = size(A, 1)
       best_side = falses(n); best_cut = 0.0
       for _ in 1:restarts
           side = bitrand(n)
           improved = true
           while improved
               improved = false
               for v in 1:n
                   δ = sum(j -> (side[v] == side[j] ? 1 : -1) * A[v, j], 1:n)
                   if δ > 0
                       side[v] = !side[v]
                       improved = true
                   end
               end
           end
           c = cut_value(A, side)
           if c > best_cut
               best_cut = c; best_side = copy(side)
           end
       end
       (cut=best_cut, side=best_side)
   end

4. 04
   Run it. For exact solutions on small graphs, branch-and-bound; for large ones, Goemans–Williamson SDP via Convex.jl; for cutting-edge research, IBM's Qiskit ports this exact problem to quantum hardware via QAOA. Julia's strength is that all three approaches share the same data types.

   julia> A = build_graph(20, 0.3);
   julia> r = local_search(A)
   (cut = 21.0, side = Bool[1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0])