Генерация случайного паросочетания

[stiv1995]

Пусть $n > 2$ некоторое натуральное число и дан набор соответствий: $\{1 : \{i_1^{(1)}, i_2^{(1)}, \ldots, i_{k_1}^{(1)}\}, 2:\{i_1^{(2)}, i_2^{(2)}, \ldots, i_{k_2}^{(2)}\}, \ldots, n: \{i_1^{(n)}, i_2^{(n)}, \ldots, i_{k_n}^{(n)}\} \}$, где все $i_j^{(k)} \in \{1, \ldots, n \}$.

Задача состоит в том, чтобы выбрать перестановку $n$ элементов из числа всех допустимых перестановок случайным образом (вероятность каждой допустимой перестановки должна быть одинаковой). Перестановка $\pi$ является допустимой, если для любого $v \in \{1, 2, \ldots, n \}$ верно, что

$$\pi(v) \ne \{i_1^{(v)}, i_2^{(v)}, \ldots, i_{k_v}^{(v)}\}.$$

Иными словами, допустимая перестановка не должна переводить элемент в элемент соответствующего множества. Существование хотя бы одной допустимой перестановки гарантируется.

Идея:
1. Построить двудольный граф с $n$ элементами в каждой из долей
2. Соединить каждую вершину $v$ из левой доли со всеми вершинами из правой доли, которые не принадлежат множеству $\{i_1^{(v)}, i_2^{(v)}, \ldots, i_{k_v}^{(v)}\}$.
3. Присвоить каждому ребру в графе вес, независимо сгенерированный из равномерного распределения на отрезке $[0, 1]$
4. Найти полное паросочетание максимального веса

Верно ли, что соответствующая перестановка удовлетворяет условию задачи, то есть выбрана случайным образом (равномерно) из всех допустимых перестановок? Кажется, что да, так как из соображений симметрии непонятно, почему какая-то допустимая перестановка должна иметь большую вероятность, чем другие. Но формально показать не получается.

[slavav]

Перестановка применяется к числу?

stiv1995 в сообщении #1464049 писал(а):

$v \in \{1, 2, \ldots, n \}$ ...

$$\pi(v) \dots$$

-- 20.05.2020, 10:36 --

Ещё вопрос: вы везде работаете с множествами. Это точно не должны быть кортежи?

[stiv1995]

slavav в сообщении #1464064 писал(а):

Перестановка применяется к числу?

Да, под перестановкой здесь понимается биекция $\pi: \{1, \ldots n\} \to \{1, \ldots n\}$.

slavav в сообщении #1464064 писал(а):

Ещё вопрос: вы везде работаете с множествами. Это точно не должны быть кортежи?

Я имел ввиду, что набор соответствий сопоставляет каждому числу $v \in \{ 1, \ldots n \}$ множество чисел, куда оно не должно переходить под действием допустимой перестановки. Внешние скобки в определении можно заменить на круглые, если вы это имеете ввиду.

Как пример, пусть $n = 3$ и дан набор соответствий $(1: \{1 \}, 2: \{2 \}, 3: \{3 \})$. Тогда по условиям задачи нужно выбрать случайным образом перестановку трех элементов без неподвижных точек. В этом случае мой алгоритм работает верно.

[arseniiv]

stiv1995 в сообщении #1464153 писал(а):

Я имел ввиду, что набор соответствий сопоставляет каждому числу $v \in \{ 1, \ldots n \}$ множество чисел, куда оно не должно переходить под действием допустимой перестановки.

А, то есть в

stiv1995 в сообщении #1464049 писал(а):

$$\pi(v) \ne \{i_1^{(v)}, i_2^{(v)}, \ldots, i_{k_v}^{(v)}\}.$$

на самом деле $\notin$.

Простейшим точно корректным алгоритмом будет генерировать случайные перестановки (Фишером—Йейтсом) и отбрасывать до тех пор, пока они не удовлетворяют ограничениям. Он даже может быть вполне приемлемым, если число неудовлетворительных перестановок мало, для чего можно получить оценки. Может быть алгоритм генерации случайной перестановки можно будет модифицировать, чтобы неподходящие числа не выбирались, но притом оставалось равномерное распределение результатов.

С паросочетаниями в двудольном графе связь была бы очевиднее, если бы в условии было $\in$, но вообще конечно всё едино.

stiv1995 в сообщении #1464049 писал(а):

Кажется, что да, так как из соображений симметрии непонятно, почему какая-то допустимая перестановка должна иметь большую вероятность, чем другие. Но формально показать не получается.

Меня берёт сомнение, что перестановка будет равномерно распределённой с таким алгоритмом. Чуть позже попробую численно проверить например условия $\pi(1)\in\{1,3\}, \pi(2)\in\{1,2\}, \pi(3)\in\{1,2,3\}$ (или какие-нибудь более несимметричные).

[arseniiv]

Интуиция кажется оказалась верна, вот результаты проверки:

10_000 samples:
mean: 0.33333   stdev: 0.046388
( 0.3548  0.3651  0.2801 )
100_000 samples:
mean: 0.33333   stdev: 0.037151
( 0.35213  0.35733  0.29054 )
1_000_000 samples:
mean: 0.33333   stdev: 0.037085
( 0.35473  0.35476  0.29051 )
10_000_000 samples:
mean: 0.33333   stdev: 0.03676
( 0.35447  0.35464  0.29089 )

Видно, что сходиться к $1/3$ они не собираются. В скобках тут частоты выпадения каждой из трёх перестановок.

-- Ср май 20, 2020 21:09:51 --

Код, которым проверялось, на всякий случай:

Используется синтаксис Python

from __future__ import annotations
from typing import Tuple, FrozenSet, Iterable, Iterator
from random import random
import statistics as s
from itertools import permutations

MatchingRestrictions = FrozenSet[Tuple[int, int]]
Permutation = Tuple[Tuple[int, int], ...]

def restricted_permutations(restrictions: MatchingRestrictions) -> Iterator[Permutation]:
    keys_tuple = tuple(set(x for x, _ in restrictions))
    for perm in permutations(keys_tuple):
        candidate = tuple(zip(keys_tuple, perm))
        if restrictions.issuperset(candidate):
            yield candidate

def matchings_v1(restricted_perms: Iterable[Permutation],
                 restrictions: MatchingRestrictions) -> Permutation:
    weights = {pair: random() for pair in restrictions}
    cur_weight = -1.0
    cur_perm = None
    for perm in restricted_perms:
        weight = sum(weights[pair] for pair in perm)
        if weight > cur_weight:
            cur_weight, cur_perm = weight, perm
    assert cur_perm is not None
    return cur_perm

def divided_ranges(first_stop: int, mult: int, max_n: int) -> Iterator[range]:
    current_start = 0
    current_stop = first_stop
    while current_stop < max_n:
        yield range(current_start, current_stop)
        current_start = current_stop
        current_stop *= mult
    yield range(current_start, max_n)

def main() -> None:
    restrictions_dict = {1: {1, 3}, 2: {1, 2}, 3: {1, 2, 3}}
    restrictions = frozenset((x, y) for x, ys in restrictions_dict.items()
                                    for y in ys)
    perms = tuple(restricted_permutations(restrictions))
    perm_counts = {perm: 0 for perm in perms}

    def stats(total: int) -> None:
        total_ = sum(perm_counts.values())
        assert total == total_
        perm_freqs = tuple(count / total for count in perm_counts.values())
        mean = s.fmean(perm_freqs)
        stdev = s.stdev(perm_freqs, mean)
        print(f'{total:_} samples:')
        print(f'mean: {mean:2.5}   stdev: {stdev:2.5}')
        perm_freqs_fmt = (f'{x:2.5}' for x in perm_freqs)
        print('(', '  '.join(perm_freqs_fmt), ')', sep=' ')

    for r in divided_ranges(10_000, 10, 10_000_000):
        for _ in r:
            perm_counts[matchings_v1(perms, restrictions)] += 1
        stats(r.stop)

if __name__ == '__main__':
    main()

[mihaild]

И даже понятно, почему получается по-разному. Перестановку $123$ мы получим, если максимальной оказалась сумма $x_{11}+x_{22}+x_{33}$, $312$ - $x_{13}+x_{21}+x_{32}$, $321$ - $x_{13} + x_{22} + x_{31}$.
Последние слагаемые можно выкинуть - если у нас есть три случайных величины, из которых какая-то оказывается максимальной чаще остальных, то после добавления к каждой из них своей (все добавляемые независимы от всего остального и распределены одинаково), то та, которая была максимальной чаще, всё еще будет распределена чаще.
После выкидывания, получаем что для вывода в лидеры $123$ нужны условия $x_{11} > x_{13}$ и $x_{11} + x_{22} > x_{13} + x_{21}$, а вот для вывода в лидеры $321$ - $x_{13} > x_{11}$ и $x_{22} > x_{21}$. Первые условия симметричные, а вот второе при условии первого для $123$ более вероятно: если у нас уже оказалось, что одно слагаемое в левой части больше слагаемого в правой, то нам достаточно, но не необходимо, чтобы и второе оказалось больше.

[arseniiv]

Ага, я хотел сначала проанализировать этот случай и прям на месте, но не осилил и стал писать численный тест. Как хорошо, что доказали точно!

[stiv1995]

Да, спасибо большое, мое предположение действительно не верное :)

Я нашел, что вообще эта задача известна в приложении к оценкам перманента матриц. Вот здесь обсуждают приближенный и нетривиальный алгоритм на основе цепей Маркова https://www.cc.gatech.edu/~vigoda/Permanent.pdf?fbclid=IwAR1kYWBh1oUpqHOjfF-GgVbegvCTZULut788TeGOtX3qgRaxQDpw7gG9PJw.

Спасибо!!