matrix - Algorithm for equiprobable random square binary matrices with two non-adjacent non-zeros in each row and column

Question

It would be great if someone could point me towards an algorithm that would allow me to :

1. create a random square matrix, with entries 0 and 1, such that
2. every row and every column contain exactly two non-zero entries,
3. two non-zero entries cannot be adjacent,
4. all possible matrices are equiprobable.

Right now I manage to achieve points 1 and 2 doing the following : such a matrix can be transformed, using suitable permutations of rows and columns, into a diagonal block matrix with blocks of the form

1 1 0 0 ... 0
0 1 1 0 ... 0
0 0 1 1 ... 0
.............
1 0 0 0 ... 1

So I start from such a matrix using a partition of [0, ..., n-1] and scramble it by permuting rows and columns randomly. Unfortunately, I can't find a way to integrate the adjacency condition, and I am quite sure that my algorithm won't treat all the matrices equally.

Update

I have managed to achieve point 3. The answer was actually straight under my nose : the block matrix I am creating contains all the information needed to take into account the adjacency condition. First some properties and definitions:

• a suitable matrix defines permutations of [1, ..., n] that can be build like so: select a 1 in row 1. The column containing this entry contains exactly one other entry equal to 1 on a row a different from 1. Again, row a contains another entry 1 in a column which contains a second entry 1 on a row b, and so on. This starts a permutation 1 -> a -> b ....

For instance, with the following matrix, starting with the marked entry

v
1 0 1 0 0 0 | 1
0 1 0 0 0 1 | 2
1 0 0 1 0 0 | 3
0 0 1 0 1 0 | 4
0 0 0 1 0 1 | 5
0 1 0 0 1 0 | 6
------------+--
1 2 3 4 5 6 |

we get permutation 1 -> 3 -> 5 -> 2 -> 6 -> 4 -> 1.

• the cycles of such a permutation lead to the block matrix I mentioned earlier. I also mentioned scrambling the block matrix using arbitrary permutations on the rows and columns to rebuild a matrix compatible with the requirements.

But I was using any permutation, which led to some adjacent non-zero entries. To avoid that, I have to choose permutations that separate rows (and columns) that are adjacent in the block matrix. Actually, to be more precise, if two rows belong to a same block and are cyclically consecutive (the first and last rows of a block are considered consecutive too), then the permutation I want to apply has to move these rows into non-consecutive rows of the final matrix (I will call two rows incompatible in that case).

So the question becomes : How to build all such permutations ?

The simplest idea is to build a permutation progressively by randomly adding rows that are compatible with the previous one. As an example, consider the case n = 6 using partition 6 = 3 + 3 and the corresponding block matrix

1 1 0 0 0 0 | 1
0 1 1 0 0 0 | 2
1 0 1 0 0 0 | 3
0 0 0 1 1 0 | 4
0 0 0 0 1 1 | 5
0 0 0 1 0 1 | 6
------------+--
1 2 3 4 5 6 |

Here rows 1, 2 and 3 are mutually incompatible, as are 4, 5 and 6. Choose a random row, say 3.

We will write a permutation as an array: [2, 5, 6, 4, 3, 1] meaning 1 -> 2, 2 -> 5, 3 -> 6, ... This means that row 2 of the block matrix will become the first row of the final matrix, row 5 will become the second row, and so on.

Now let's build a suitable permutation by choosing randomly a row, say 3:

• p = [3, ...]

The next row will then be chosen randomly among the remaining rows that are compatible with 3 : 4, 5and 6. Say we choose 4:

• p = [3, 4, ...]

Next choice has to be made among 1 and 2, for instance 1:

• p = [3, 4, 1, ...]

And so on: p = [3, 4, 1, 5, 2, 6].

Applying this permutation to the block matrix, we get:

1 0 1 0 0 0 | 3
0 0 0 1 1 0 | 4
1 1 0 0 0 0 | 1
0 0 0 0 1 1 | 5
0 1 1 0 0 0 | 2
0 0 0 1 0 1 | 6
------------+--
1 2 3 4 5 6 |

Doing so, we manage to vertically isolate all non-zero entries. Same has to be done with the columns, for instance by using permutation p' = [6, 3, 5, 1, 4, 2] to finally get

0 1 0 1 0 0 | 3
0 0 1 0 1 0 | 4
0 0 0 1 0 1 | 1
1 0 1 0 0 0 | 5
0 1 0 0 0 1 | 2
1 0 0 0 1 0 | 6
------------+--
6 3 5 1 4 2 | 

So this seems to work quite efficiently, but building these permutations needs to be done with caution, because one can easily be stuck: for instance, with n=6 and partition 6 = 2 + 2 + 2, following the construction rules set up earlier can lead to p = [1, 3, 2, 4, ...]. Unfortunately, 5 and 6 are incompatible, so choosing one or the other makes the last choice impossible. I think I've found all situations that lead to a dead end. I will denote by r the set of remaining choices:

1. p = [..., x, ?], r = {y} with x and y incompatible
2. p = [..., x, ?, ?], r = {y, z} with y and z being both incompatible with x (no choice can be made)
3. p = [..., ?, ?], r = {x, y} with x and y incompatible (any choice would lead to situation 1)
4. p = [..., ?, ?, ?], r = {x, y, z} with x, y and z being cyclically consecutive (choosing x or z would lead to situation 2, choosing y to situation 3)
5. p = [..., w, ?, ?, ?], r = {x, y, z} with xwy being a 3-cycle (neither x nor y can be chosen, choosing z would lead to situation 3)
6. p = [..., ?, ?, ?, ?], r = {w, x, y, z} with wxyz being a 4-cycle (any choice would lead to situation 4)
7. p = [..., ?, ?, ?, ?], r = {w, x, y, z} with xyz being a 3-cycle (choosing w would lead to situation 4, choosing any other would lead to situation 4)

Now it seems that the following algorithm gives all suitable permutations:

• As long as there are strictly more than 5 numbers to choose, choose randomly among the compatible ones.
• If there are 5 numbers left to choose: if the remaining numbers contain a 3-cycle or a 4-cycle, break that cycle (i.e. choose a number belonging to that cycle).
• If there are 4 numbers left to choose: if the remaining numbers contain three cyclically consecutive numbers, choose one of them.
• If there are 3 numbers left to choose: if the remaining numbers contain two cyclically consecutive numbers, choose one of them.

I am quite sure that this allows me to generate all suitable permutations and, hence, all suitable matrices.

Unfortunately, every matrix will be obtained several times, depending on the partition that was chosen.

Answer 1

Intro

Here is some prototype-approach, trying to solve the more general task of uniform combinatorial sampling, which for our approach here means: we can use this approach for everything which we can formulate as SAT-problem.

It's not exploiting your problem directly and takes a heavy detour. This detour to the SAT-problem can help in regards to theory (more powerful general theoretical results) and efficiency (SAT-solvers).

That being said, it's not an approach if you want to sample within seconds or less (in my experiments), at least while being concerned about uniformity.

Theory

The approach, based on results from complexity-theory, follows this work:

GOMES, Carla P.; SABHARWAL, Ashish; SELMAN, Bart. Near-uniform sampling of combinatorial spaces using XOR constraints. In: Advances In Neural Information Processing Systems. 2007. S. 481-488.

The basic idea:

• formulate the problem as SAT-problem
• add randomly generated xors to the problem (acting on the decision-variables only! that's important in practice)
  • this will reduce the number of solutions (some solutions will get impossible)
  • do that in a loop (with tuned parameters) until only one solution is left!
    • search for some solution is being done by SAT-solvers or #SAT-solvers (=model-counting)
    • if there is more than one solution: no xors will be added but a complete restart will be done: add random-xors to the start-problem!

The guarantees:

• when tuning the parameters right, this approach achieves near-uniform sampling
  • this tuning can be costly, as it's based on approximating the number of possible solutions
  • empirically this can also be costly!
  • Ante's answer, mentioning the number sequence A001499 actually gives a nice upper bound on the solution-space (as it's just ignoring adjacency-constraints!)

The drawbacks:

• inefficient for large problems (in general; not necessarily compared to the alternatives like MCMC and co.)
  • need to change / reduce parameters to produce samples
  • those reduced parameters lose the theoretical guarantees
  • but empirically: good results are still possible!

Parameters:

In practice, the parameters are:

• N: number of xors added
• L: minimum number of variables part of one xor-constraint
• U: maximum number of variables part of one xor-constraint

N is important to reduce the number of possible solutions. Given N constant, the other variables of course also have some effect on that.

Theory says (if i interpret correctly), that we should use L = R = 0.5 * #dec-vars.

This is impossible in practice here, as xor-constraints hurt SAT-solvers a lot!

Here some more scientific slides about the impact of L and U.

They call xors of size 8-20 short-XORS, while we will need to use even shorter ones later!

Implementation

Final version

Here is a pretty hacky implementation in python, using the XorSample scripts from here.

The underlying SAT-solver in use is Cryptominisat.

The code basically boils down to:

• Transform the problem to conjunctive normal-form
  • as DIMACS-CNF
• Implement the sampling-approach:
  • Calls XorSample (pipe-based + file-based)
  • Call SAT-solver (file-based)
• Add samples to some file for later analysis

Code: (i hope i did warn you already about the code-quality)

from itertools import count
from time import time
import subprocess
import numpy as np
import os
import shelve
import uuid
import pickle
from random import SystemRandom
cryptogen = SystemRandom()

""" Helper functions """
# K-ARY CONSTRAINT GENERATION
# ###########################
# SINZ, Carsten. Towards an optimal CNF encoding of boolean cardinality constraints.
# CP, 2005, 3709. Jg., S. 827-831.

def next_var_index(start):
    next_var = start
    while(True):
        yield next_var
        next_var += 1

class s_index():
    def __init__(self, start_index):
        self.firstEnvVar = start_index

    def next(self,i,j,k):
        return self.firstEnvVar + i*k +j

def gen_seq_circuit(k, input_indices, next_var_index_gen):
    cnf_string = ''
    s_index_gen = s_index(next_var_index_gen.next())

    # write clauses of first partial sum (i.e. i=0)
    cnf_string += (str(-input_indices[0]) + ' ' + str(s_index_gen.next(0,0,k)) + ' 0
')
    for i in range(1, k):
        cnf_string += (str(-s_index_gen.next(0, i, k)) + ' 0
')

    # write clauses for general case (i.e. 0 < i < n-1)
    for i in range(1, len(input_indices)-1):
        cnf_string += (str(-input_indices[i]) + ' ' + str(s_index_gen.next(i, 0, k)) + ' 0
')
        cnf_string += (str(-s_index_gen.next(i-1, 0, k)) + ' ' + str(s_index_gen.next(i, 0, k)) + ' 0
')
        for u in range(1, k):
            cnf_string += (str(-input_indices[i]) + ' ' + str(-s_index_gen.next(i-1, u-1, k)) + ' ' + str(s_index_gen.next(i, u, k)) + ' 0
')
            cnf_string += (str(-s_index_gen.next(i-1, u, k)) + ' ' + str(s_index_gen.next(i, u, k)) + ' 0
')
        cnf_string += (str(-input_indices[i]) + ' ' + str(-s_index_gen.next(i-1, k-1, k)) + ' 0
')

    # last clause for last variable
    cnf_string += (str(-input_indices[-1]) + ' ' + str(-s_index_gen.next(len(input_indices)-2, k-1, k)) + ' 0
')

    return (cnf_string, (len(input_indices)-1)*k, 2*len(input_indices)*k + len(input_indices) - 3*k - 1)

# K=2 clause GENERATION
# #####################
def gen_at_most_2_constraints(vars, start_var):
    constraint_string = ''
    used_clauses = 0
    used_vars = 0
    index_gen = next_var_index(start_var)
    circuit = gen_seq_circuit(2, vars, index_gen)
    constraint_string += circuit[0]
    used_clauses += circuit[2]
    used_vars += circuit[1]
    start_var += circuit[1]

    return [constraint_string, used_clauses, used_vars, start_var]

def gen_at_least_2_constraints(vars, start_var):
    k = len(vars) - 2
    vars = [-var for var in vars]

    constraint_string = ''
    used_clauses = 0
    used_vars = 0
    index_gen = next_var_index(start_var)
    circuit = gen_seq_circuit(k, vars, index_gen)
    constraint_string += circuit[0]
    used_clauses += circuit[2]
    used_vars += circuit[1]
    start_var += circuit[1]

    return [constraint_string, used_clauses, used_vars, start_var]

# Adjacency conflicts
# ###################
def get_all_adjacency_conflicts_4_neighborhood(N, X):
    conflicts = set()
    for x in range(N):
        for y in range(N):
            if x < (N-1):
                conflicts.add(((x,y),(x+1,y)))
            if y < (N-1):
                conflicts.add(((x,y),(x,y+1)))

    cnf = ''  # slow string appends
    for (var_a, var_b) in conflicts:
        var_a_ = X[var_a]
        var_b_ = X[var_b]
        cnf += '-' + var_a_ + ' ' + '-' + var_b_ + ' 0 
'

    return cnf, len(conflicts)

# Build SAT-CNF
  #############
def build_cnf(N, verbose=False):
    var_counter = count(1)
    N_CLAUSES = 0
    X = np.zeros((N, N), dtype=object)
    for a in range(N):
        for b in range(N):
            X[a,b] = str(next(var_counter))

    # Adjacency constraints
    CNF, N_CLAUSES = get_all_adjacency_conflicts_4_neighborhood(N, X)

    # k=2 constraints
    NEXT_VAR = N*N+1

    for row in range(N):
        constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_most_2_constraints(X[row, :].astype(int).tolist(), NEXT_VAR)
        N_CLAUSES += used_clauses
        CNF += constraint_string

        constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_least_2_constraints(X[row, :].astype(int).tolist(), NEXT_VAR)
        N_CLAUSES += used_clauses
        CNF += constraint_string

    for col in range(N):
        constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_most_2_constraints(X[:, col].astype(int).tolist(), NEXT_VAR)
        N_CLAUSES += used_clauses
        CNF += constraint_string

        constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_least_2_constraints(X[:, col].astype(int).tolist(), NEXT_VAR)
        N_CLAUSES += used_clauses
        CNF += constraint_string

    # build final cnf
    CNF = 'p cnf ' + str(NEXT_VAR-1) + ' ' + str(N_CLAUSES) + '
' + CNF

    return X, CNF, NEXT_VAR-1


# External tools
# ##############
def get_random_xor_problem(CNF_IN_fp, N_DEC_VARS, N_ALL_VARS, s, min_l, max_l):
    # .cnf not part of arg!
    p = subprocess.Popen(['./gen-wff', CNF_IN_fp,
                          str(N_DEC_VARS), str(N_ALL_VARS),
                          str(s), str(min_l), str(max_l), 'xored'],
                          stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    result = p.communicate()

    os.remove(CNF_IN_fp + '-str-xored.xor')  # file not needed
    return CNF_IN_fp + '-str-xored.cnf'

def solve(CNF_IN_fp, N_DEC_VARS):
    seed = cryptogen.randint(0, 2147483647)  # actually no reason to do it; but can't hurt either
    p = subprocess.Popen(["./cryptominisat5", '-t', '4', '-r', str(seed), CNF_IN_fp], stdin=subprocess.PIPE, stdout=subprocess.PIPE)
    result = p.communicate()[0]

    sat_line = result.find('s SATISFIABLE')

    if sat_line != -1:
        # solution found!
        vars = parse_solution(result)[:N_DEC_VARS]

        # forbid solution (DeMorgan)
        negated_vars = list(map(lambda x: x*(-1), vars))
        with open(CNF_IN_fp, 'a') as f:
            f.write( (str(negated_vars)[1:-1] + ' 0
').replace(',', ''))

        # assume solve is treating last constraint despite not changing header!
        # solve again

        seed = cryptogen.randint(0, 2147483647)
        p = subprocess.Popen(["./cryptominisat5", '-t', '4', '-r', str(seed), CNF_IN_fp], stdin=subprocess.PIPE, stdout=subprocess.PIPE)
        result = p.communicate()[0]
        sat_line = result.find('s SATISFIABLE')
        if sat_line != -1:
            os.remove(CNF_IN_fp)  # not needed anymore
            return True, False, None
        else:
            return True, True, vars
    else:
        return False, False, None

def parse_solution(output):
    # assumes there is one
    vars = []
    for line in output.split("
"):
        if line:
            if line[0] == 'v':
                line_vars = list(map(lambda x: int(x), line.split()[1:]))
                vars.extend(line_vars)
    return vars

# Core-algorithm
# ##############
def xorsample(X, CNF_IN_fp, N_DEC_VARS, N_VARS, s, min_l, max_l):
    start_time = time()
    while True:
        # add s random XOR constraints to F
        xored_cnf_fp = get_random_xor_problem(CNF_IN_fp, N_DEC_VARS, N_VARS, s, min_l, max_l)
        state_lvl1, state_lvl2, var_sol = solve(xored_cnf_fp, N_DEC_VARS)

        print('------------')

        if state_lvl1 and state_lvl2:
            print('FOUND')

            d = shelve.op