MATLAB: Knuth’s Algorithm X (Dancing Links; Exact Cover)

dancing linksexact coverknuth's algorithm x

I'm trying to write a Matlab version of Knuth's "Algorithm X." This algorithm makes use of the Dancing Links routine and is a recursive, nondeterministic, depth-first, backtracking algorithm. I am using this algorithm to solve a simple application of the exact cover problem. For more information:
Knuth's Algorithm X: https://en.wikipedia.org/wiki/Knuth%27s_Algorithm_X
Exact cover problem: https://en.wikipedia.org/wiki/Exact_cover
Here is pseudocode for Algorithm X:
 1. If matrix A is empty, the problem is solved; terminate successfully.
 2. Otherwise: choose a column, c, of A with the least number of 1's.
 3. Choose a row, r, that contains a 1 in column c (e.g., A(r,c) = 1).
 4. Add r to the partial solution.
 5. For each column j such that A(r,j) = 1:
         for each row i such that A(i,j) = 1:
              delete row i from matrix A
         delete column j from matrix A
 6. Repeat this algorithm recursively on the reduced matrix A.
Following the Wiki-page, here is a working example:
A = zeros(1,7); B = A; C = B; D = C; E = D; F = E;
A([1 4 7]) = 1; B([1 4]) = 1; C([4 5 7]) = 1; D([3 5 6]) = 1; E([2 3 6 7]) = 1; F([2 7]) = 1; 
M = [A; B; C; D; E; F]; % the binary 'A' matrix
Knuth's algorithm should return the following exact cover of M: {B, D, F}; i.e., rows 2, 4, and 6 of the original (unreduced) matrix collectively contribute a 1 in each column of M (exactly once).
Here is what I have so far:
function sol = dlx(X,sol)
% if isempty(X)
%     return



% end



% Find the first column of X having the lowest number of 1s in any position
c = find(sum(X,1) == min(sum(X,1)), 1);
% If column c is entirely zero, terminate unsuccessfully
% if (sum(X(:,c)) == 0)
%     sol = [];

%     return
% end
% Find the rows that have a 1 in column c
r = find(X(:,c) == 1);
% if isempty(r)

%     sol = [];
%     return
% end
% if isempty(r)
%     sol = sol(1:end-1);
%     return
% end
% Loop over each row that has a 1 in column c
for rr = 1:length(r)
      % include row r in the partial solution
      sol = [sol r(rr)];
      % find the columns in which A(r,j) = 1
      J = find(X(r(rr),:) == 1);
      % initialize the reduced matrix
      Xred = X;
      I = [];
      for jj = 1:length(J)
          I = [I; find(Xred(:,J(jj)) == 1)];
      end
      I = unique(I);
      Xred(I,:) = []; % delete the I rows
      Xred(:,J) = []; % delete the J columns
      % repeat this algorithm recursively on the reduced matrix Xred.
      sol = dlx(Xred, sol);
end
end
For the above example,
sol = dlx(X, [])
sol =
       1     2     1     1
Except that my algorithm should return sol = [2 4 6].
Here is another simple example that also bugs my code:
A = [1 1 0; 0 1 1];
sol = dlx(A, [])
sol =
     1
In this case, the there is no exact cover of the matrix A, and so my algorithm should return an empty solution (e.g., sol = []).
As you can see, I am very close. I believe that I have successfully written the recursion of reducing the matrix, but I am seeking help on the following:
  1. Looking at the pseudocode, should the "Choose a row, r, that contains a 1 in column c" be implemented as a for-loop in my code? That is, should I loop over all such rows?
  2. I'm having trouble formulating the final solution from the partial solution (e.g., so that my function returns sol = [2 4 6]). During the recursion, the partial solution should 'backtrack' when the partial solution turns out not to be valid.
  3. Where should I place the 'if A is empty' condition from the pseudocode? You'll see that I've tried putting this condition at the beginning of the function, and also at the 'isempty(r)' line.
  4. I'm also open to suggestions for making my code more efficient; i.e., should/can I replace the 'find' or 'for-loops' with better code?
Can any one offer me suggestions to my questions? Thanks in advance!

Best Answer

Ok, just read the description of the algorithm on wikipedia. A few comments:
the test for emptiness is the first step according to the algorithm so that should be the first thing you should do. However, if require that the initial matrix is never empty you could move that test to just before calling the recursion since you know that the recursion will only return success. I'm not sure that would make the code any cleaner.
As it is at the moment, your code doesn't check for successful or unsuccessful recursion. You also need to be able to distinguish between success and failure so you need to either add an extra output or find a way to return a different sol in each case.
Furthermore, the algorithm should be able to return multiple solutions if they exist (for example for A = [1 0 0; 0 1 1; 0 1 0; 1 0 1; 0 0 1; 1 1 0] there are 4 solutions). So you can't just return a vector. Since the different solution may have different number of rows, you'd have to return a cell array.
As Thomas says in his comment, you also fails to track the actual row numbers after you've deleted rows.
In terms of loops, you have too many. The deletions can be done without loops (and even without find).
Here is how I'd implement it. I choose to move the recursion into its own function because of the extra argument needed to track the row indices. The distinction between success and failure is done by the sol return value. If it's an empty cell array it's a failure, otherwise it's a success. If the matrix was empty the cell array contains an empty matrix otherwise it contains the possible row combinations so far.
function [sol, success] = dlx(A)
    sol = recursion(A, 1:size(A, 1));
end
function [sol] = recursion(A, rownum)
    sol = {};   %empty cell array indicates failure   
    %step 1
    if isempty(A)
        %success
        sol = {[]};
        return;
    end
    %step 2
    [~, c] = min(sum(A, 1));
    if ~any(A(:, c))
        %failure
        return;
    end
    %step 3
    possiblerows = find(A(:, c));
    for tryrow = possiblerows.'
        %step 4 for later, only if success
        %step 5
        deletecols = A(tryrow, :) == 1;
        deleterows = any(A(:, deletecols) == 1, 2);
        reducedA = A(:,~deletecols);
        reducedA = reducedA(~deleterows, :);
        reducedrows = rownum(~deleterows);
        %step 6
        [subsol] = recursion(reducedA, reducedrows);
        if ~isempty(subsol)
            %step 4 
            %add current row to the all the row combinations returned by the recursion
            sol = [sol; cellfun(@(ss) [rownum(tryrow), ss], subsol, 'UniformOutput', false)]; %#ok<AGROW>
        end
    end
    %if all the tryrow ended in failure then sol is still an empty cell array which indicates failure on return
end
Related Question