Couverture exacte, algorithme X et Dancing Links

Le problème de la couverture exacte

Le problème de la couverture exacte, énoncé par Karp dans l'article Reducibility Among Combinatorial Problems, peut être formulé de deux manières équivalentes:

Mais peu importe la formulation, il s'agit d'un problème difficile à résoudre et depuis longtemps connu comme étant NP complet.

L'algorithme X

L'algorithme X, décrit par Knuth dans l'article Dancing Links, permet de solutionner le problème de la couverture exacte. Il s'agit d'un algorithme de recherche récursif utilisant du retour sur trace. Tout comme le problème de la couverture exacte, celui-ci peut être énoncé de deux manières équivalentes.

Une première formulation de l'algorithme, correspondant à la première formulation du problème de la couverture exacte, opère sur un ensemble d'éléments U et une famille F de sous-ensembles de U:

  • si l'ensemble U est vide, le problème est résolu et retourner avec succès;
  • choisir arbitrairement un élément x dans l'ensemble U;
  • pour chaque sous-ensemble S dans F contenant x:
    • former une famille d'ensembles F' comprenant tous les ensembles dans F qui sont disjoints de S;
    • résoudre récursivement la couverture exacte pour U/S, F';
    • si une solution a été trouvée, ajouter S à la solution et retourner.

Une seconde formulation de l'algorithme, correspondant à la seconde formulation du problème de la couverture exacte, opère sur des matrices binaires:

  • si la matrice A est vide, le problème est résolu et retourner avec succès;
  • choisir une colonne c dans A;
  • choisir une ligne r dans la colonne c;
  • inclure la ligne r dans la solution partielle;
  • pour chaque indice j tel que A[r, j] = 1:
    • supprimer la colonne d'indice j de la matrice A;
    • pour chaque indice i tel que A[i, j] = 1:
      • supprimer la ligne d'indice i de la matrice A;
  • répéter l'algorithme récursivement sur la matrice réduite A.

Une analyse très intéressante de l'algorithme X peut être consultée ici.

Dancing Links

La technique des Dancing Links permet d'implanter la formulation de l'algorithme X opérant sur des matrices binaires d'une manière très ingénieuse. Au cœur de la technique se trouve une structure de données quadruplement chaînée:

Matrice binaire encodant le problème de couverture exacte.

Structure chaînée représentant la matrice binaire encodant le problème de couverture exacte.

Cette structure de données particulière permet essentiellement deux choses:

Au demeurant, la lecture de la première partie de l'article de Knuth est certainement recommandée pour une bonne compréhension de la suite de ce texte.

Implantation

Matrices clairsemées

Comme les matrices binaires encodant les instances du problème de la couverture exacte sont la plupart du temps clairsemées, les représenter avec des tableaux de vecteurs serait inutilement coûteux en espace mémoire. Pour diminuer ce coût, nous utiliserons plutôt un vecteur de std::map. Aussi, seul les éléments non nuls d'une ligne y seront stockés.

L'implantation de cette structure de données se divise en deux parties: le conteneur et un proxy. L'utilisation d'un proxy est nécessaire car nous devons être en mesure de retourner à l'utilisateur de la classe des références non constantes vers des éléments qui n'existent pas.

template <typename T>
class MappedMatrix {
public:
    MappedMatrix(size_t row, size_t col) :
        size(make_subscript(row, col)), zero() {
        data.resize(row);
    }

    const MappedMatrix& operator=(MappedMatrix rhs) {
        std::swap(data, rhs.data);
        size = rhs.size;
        return *this;
    }

    const T& operator()(uint32 row, uint32 col) const {
        if (row >= size.row || col >= size.col)
            throw std::invalid_argument("Invalid subscripts.");
        return value_at(make_subscript(row, col));
    }

    Proxy<T, MappedMatrix> operator()(uint32 row, uint32 col) {
        if (row >= size.row || col >= size.col)
            throw std::invalid_argument("Invalid subscripts.");
        return Proxy<T, MappedMatrix>(*this, make_subscript(row, col));
    }

    uint32 rows() const { return size.row; }
    uint32 cols() const { return size.col; }

private:
    template <typename, template <typename> class>
    friend class Proxy;

    //! \brief Set the value of an element at a given subscript.
    void set_value(const Subscript<uint32>& subscript, const T& value) {
        std::map<uint32, T>& row = data[subscript.row];

        typedef typename std::map<uint32, T>::iterator iterator;
        iterator none = row.end();
        iterator elem = row.find(subscript.col);

        // Inserting a non-zero value over a zero value.
        // Inserting a non-zero value over a non-zero value.
        // Inserting a zero value over a non-zero value.
        // Inserting a zero value over a zero value: NOP.
        if (value != zero && elem == none)
            row.insert(std::make_pair(subscript.col, value));
        else if (value != zero && elem != none)
            elem->second = value;
        else if (value == zero && elem != none)
            row.erase(elem);
    }

    const T& value_at(const Subscript<uint32>& subscript) const {
        typedef typename std::map<uint32, T>::const_iterator iterator;
        const std::map<uint32, T>& row = data[subscript.row];
        iterator iter = row.find(subscript.col);
        return iter != row.end() ? iter->second : zero;
    }

    Subscript<uint32> size;
    std::vector<std::map<uint32, T> > data;
    T zero;
};

Le proxy prend en argument une référence sur la matrice qui l'a instancié ainsi que l'indice de l'élément dont il joue le rôle. Afin de demeurer transparent, le proxy surcharge tous les opérateurs pouvant être appliqués au type sous-jacent. Ces surcharges ne font, en pratique, que faire suivre l'appel.

template <typename M, template <typename M> class Matrix>
class Proxy {
public:
    Proxy(Matrix<M>& matrix, const Subscript<uint32>& subscript) :
        matrix(matrix), subscript(subscript) {
    }

    Proxy& operator=(const Proxy& rhs) {
        matrix.set_value(subscript, M(rhs));
        return *this;
    }

    template <typename A>
    Proxy& operator=(const A& value) {
        matrix.set_value(subscript, value);
        return *this;
    }

    template <typename A>
    bool operator==(const A& rhs) const {
        return matrix.value_at(subscript) == M(rhs);
    }

    template <typename A>
    bool operator!=(const A& rhs) const {
        return matrix.value_at(subscript) != M(rhs);
    }

    /* ... */

    operator const M&() const {
        return matrix.value_at(subscript);
    }

private:
    Matrix<M>& matrix;
    Subscript<uint32> subscript;
};

Algorithme X

Avant d'exécuter l'algorithme X, il importe d'allouer la structure de données chaînée sur laquelle ce dernier opérera. Les nœuds de cette structure sont représentés par un struct contenant quatre pointeurs vers ses voisins immédiats, un pointeur vers l'entête de sa colonne et un champ data. Selon que le nœud est utilisé en tant qu'en-tête de colonne ou comme élément de la structure chaînée, le champ data stockera soit le nombre d'éléments appartenant à la colonne ou bien le numéro de la ligne d'origine du nœud dans la matrice binaire encodant l'instance du problème de couverture exacte (une information nécessaire afin de reconstruire la solution).

struct Node {
    Node* header;
    Node* up, *down;
    Node* left, *right;
    uint32 data;
};

La méthode build_cover_matrix prend en entrée une matrice binaire et construit, ligne par ligne, la structure chaînée de la matrice de couverture correspondante. Évidemment, puisqu'il s'agit d'un template, l'argument n'a pas a véritablement être une matrice, mais n'a qu'à implanter l'interface idoine. Comme nous le verrons dans un texte subséquent, cette caractéristique ouvre la porte à des optimisations intéressantes.

template <typename Matrix>
Node* build_cover_matrix(const Matrix& matrix)  {
    Node* root = build_root();

    // These vectors will be used to get fast access to row and
    // column headers during the initialization phase of the cover
    // matrix. Row headers will immediately be deleted at the
    // end of the initialization process.
    std::vector<Node*> row_headers = build_row_headers(root, matrix.rows());
    std::vector<Node*> column_headers = build_column_headers(root, matrix.cols());

    for (size_t row = 0; row < matrix.rows(); row++) {

        Node* row_first = row_headers[row];
        Node* row_last = row_headers[row];

        for (size_t col = 0; col < matrix.cols(); col++) {
            if (matrix(row, col)) {
                Node* header = column_headers[col];

                // Initialize new node.
                Node* node = new Node;
                node->header = header;
                node->up = header->up;
                node->down = header;
                node->left = row_last;
                node->right = row_first;
                node->data = row;

                // Update neighbors.
                node->up->down = node;
                node->left->right = node;

                header->up = node;
                header->data++;

                row_first->left = node;
                row_last = node;
            }
        }
    }

    delete_row_headers(row_headers);

    return root;
}

La méthode associée, delete_cover_matrix permet de désallouer la matrice de couverture une fois le problème solutionné.

void delete_cover_matrix(Node* root) {
    Node* header = root->right;
    Node* old, *node;

    // Go over each header of the cover matrix.
    while (header != root) {
        node = header->down;

        // Deallocate each element in the
        // column associated to this header.
        while (node != header) {
            old = node;
            node = node->down;
            delete old;
        }

        // Deallocate the header.
        old = header;
        header = header->right;
        delete old;
    }

    delete root;
}

La méthode solve implante l'algorithme X à proprement parler. Elle prend en argument un pointeur sur la racine de la matrice de couverture ainsi qu'une référence sur un vecteur dans lequel l'indice des lignes (de la matrice binaire encodant le problème de couverture exacte) formant la solution seront stockés.

bool solve(Node* root, std::vector<uint32>& cover) {
    bool solved = false;

    Node* header = choose_next_column(root);

    if (header == root)
        return true;

    cover_column(header);

    Node* column_element = header->down;
    while (column_element != header) {
        Node* row_element = column_element->right;
        while (row_element != column_element) {
            cover_column(row_element->header);
            row_element = row_element->right;
        }

        solved = solve(root, cover);

        row_element = column_element->left;
        while (row_element != column_element) {
            uncover_column(row_element->header);
            row_element = row_element->left;
        }

        // If we've solved the exact cover problem,
        // set the cell value accordingly.
        if (solved) {
            cover.push_back(column_element->data);
            break;
        }

        column_element = column_element->down;
    }

    uncover_column(header);
    return solved;
}

La méthode choose_next_column implante une heuristique qui permet d'accélérer substantiellement l'exécution de l'algorithme. Tout simplement, il s'agit de toujours choisir les prochains pivots dans la colonne qui contient le moins d'éléments. Ce faisant, le degré de ramification (branching factor) demeure minimal.

Node* choose_next_column(Node* root)  {
    uint32 lower = std::numeric_limits<uint32>::max();

    Node* current = root->right;
    Node* next = current;

    while (current != root) {
        if (current->data < lower) {
            lower = current->data;
            next = current;
        }
        current = current->right;
    }

    return next;
}

Finalement, les méthodes cover_column et uncover_column permettent de supprimer et de réinsérer des colonnes au sein de la matrice de couverture; en ne touchant toujours qu'à quelques pointeurs. C'est ici que les pointeurs font leur danse.

void cover_column(Node* header) {
    header->left->right = header->right;
    header->right->left = header->left;

    Node* col = header->down;
    while (col != header) {
        Node* row = col->right;
        while (row != col) {
            row->up->down = row->down;
            row->down->up = row->up;
            row->header->data--;
            row = row->right;
        }
        col = col->down;
    }
}

void uncover_column(Node* header) {
    Node* col = header->up;
    while (col != header) {
        Node* row = col->left;
        while (row != col) {
            row->up->down = row;
            row->down->up = row;
            row->header->data++;
            row = row->left;
        }
        col = col->up;
    }

    header->left->right = header;
    header->right->left = header;
}

Conclusion

Plusieurs améliorations pourraient être apportées à cette implantation de l'algorithme X, mais la plus intéressante serait sans doute d'utiliser un pool de mémoire pour allouer les nœuds de la matrice de couverture. Cela permettrait très certainement d'accélérer l'allocation des nœuds en plus d'améliorer la localité des données.

Le code source est disponible sur github: exact-cover.

Knuth a donné une lecture sur la technique des Dancing Links. Celle-ci peut être visionnée en ligne. La version originale de son article peut être téléchargée à partir de son site web sous la forme d'un document PostScript compressé.