リリース時刻付きジョブスケジューリングの定式化

問題: $n$個のジョブを1つの機械で任意順に処理する。同時に取り組めるジョブは1つであり、1つのジョブに着手したら途中で他のジョブに切り替えることはできない。ジョブ$i$を始めてから完了するまでの時間は$p_i \ (> 0)$であり、ジョブ$i$を始めるのは時刻$r_i \ (\ge 0)$以降でなければならない。また、各ジョブには重み $w_i > 0$が与えられており、ジョブの処理を完了する時刻の重み付き総和を最小化したい。つまり、ジョブ$i$を時刻$t_i$から始めるとした時、$\sum_i w_i(t_i + p_i)$を最小化したい。最適な$t_1, ..., t_n$を求めよ。

ジョブスケジューリング問題の一種であり、分類としては$1|r_j|\sum w_jC_j$という表記で表される。NP-hardである。

これをMIPとして定式化したい。自然な考え方としては、ジョブ$i$にジョブ$j$以前に取り組むことに対応するバイナリ変数$x_{i, j}$を使った定式化が考えられる。目的関数は定数項を無視すると$\sum_i w_i t_i$であり、制約は以下のようになる:

• 任意のジョブ$j$について $t_j \ge r_j$
• 任意の異なるジョブ$i, j$について $t_j \ge t_i + p_i - Mx_{j, i}$
  • $M$は十分に大きな定数
• 任意の異なるジョブ$i, j$について $x_{i, j} + x_{j, i} = 1$
• 任意の異なるジョブ$i, j, k$について $x_{i, j} + x_{j, k} \le 1 + x_{i, k}$
• 任意の異なるジョブ$i, j$について $x_{i, j} \in \{0, 1\}$
• 任意のジョブ$j$について $t_j \in \mathbb R$

これで定式化になっているが、NemhauserとSavelsberghによるとbig Mを使わない定式化が可能らしい。一般性を失わずに$r_1 \le r_2 \le ... \le r_n$と仮定するとき、$t_j \ge r_j$と$t_j \ge t_i + p_i - Mx_{j, i}$の代わりに

$$\tag{*} t_j \ge r_ix_{i, j} + \sum_{k < i, k \neq j}p_k(x_{i, k} + x_{k, j}-1) + \sum_{k \ge i, k \neq j}p_kx_{k, j}$$

とbig Mを使わずに書けるというのが論文の内容である。ただし、便宜上$x_{i,i}=1$とする。上のbig Mを使った定式化を$\operatorname{FormM}$、$(*)$を使った定式化を$\operatorname{Form*}$とする。$\operatorname{Form*}$の正当性を確かめたい。

ジョブを処理する順番を固定して考えることにする。具体的には、次のことを示す:

定理: $1, ..., n$の順列$\pi = (\pi(1), ..., \pi(n))$が与えられたとする。冒頭の問題に「ジョブを$\pi$の順で処理する」という条件を追加したものの最適解の集合は、$\operatorname{Form*}$のすべての$x_{i, j}$を順列に従った値を取る定数としたLPの最適解の集合と一致する。

これが示されれば、

冒頭の問題の最適解の集合
= $\bigcup_{\pi: 順列}$(ジョブを$\pi$の順で処理するという条件の下での冒頭の問題の最適解の集合) の要素で目的関数を最小化するものの集合
= $\bigcup_{\pi: 順列}$($\operatorname{Form*}$のすべての$x_{i, j}$を順列に従った値を取る定数としたLPの最適解の集合) の要素で目的関数を最小化するものの集合
= $\operatorname{Form*}$の最適解の集合

となって$\operatorname{Form*}$は正しい定式化であると言える。

定理の証明

以降、ジョブ$j$を$\operatorname{pos}(j)$番目に処理するとする。$\operatorname{pos}(\pi(k)) = k$である。

まず、ジョブを$\pi(1), ..., \pi(n)$の順で処理するという条件下でのスケジューリング問題の最適解$(t_1, ..., t_n) = (T_1, ..., T_n)$の性質を考える。明らかに、可能な最も早い時刻にジョブを順次処理していくのが最適だから、$T_{\pi(k+1)} = \max(r_{\pi(k+1)}, T_{\pi(k)} + p_{\pi(k)})$である。これを利用して$T_j$を$r$と$p$のみで表すことを考える。

まず、$T_{\pi(1)}=r_{\pi(1)}$である。また、$\pi(1)$以外の任意のジョブ$j$については、$j$より前に行うジョブでリリース時刻ちょうどに処理を始めるものが必ず存在する。その中で最後に処理するジョブを$j'$とするとき

$$T_j=\max \bigl(r_j, r_{j'} + \sum_{k : \operatorname{pos}(j') \le \operatorname{pos} (k) < \operatorname{pos}(j)}p_k \bigr)$$

である。したがって、次の3つの命題が示せればよい。以下、$(*)$の右辺を $R_{i, j} := r_ix_{i, j} + \sum_{k < i, k \neq j}p_k(x_{i, k} + x_{k, j}-1) + \sum_{k \ge i, k \neq j}p_kx_{k, j}$ とする。ただし、ここでは$x_{i, j}$は与えられた順列$\pi$に応じて値が定まっている定数である。

命題1: ジョブ$j \neq \pi(1)$が与えられたとする。$j$より前に行いかつリリース時刻ちょうどに処理を始めるジョブで最後に処理するものを$j'$とするとき、$R_{j', j} = r_{j'} + \sum_{k: \operatorname{pos}(j') \le \operatorname{pos} (k) < \operatorname{pos}(j)}p_k$である。

命題2: 任意のジョブ$j$について、$T_j = r_j$ならば$R_{j, j} = r_j$である。

命題3: 任意のジョブ$i, j$について$R_{i, j} \le T_j$である。

命題1と2は$\operatorname{Form*}$が$t_j$の下限に関するタイトな制約をすべて含むと主張している。命題3はそれ以外の制約が冗長であると主張している。

命題1の証明

$R_{j', j} = r_{j'} + \sum_{k < j', k \neq j}p_k(x_{j', k} + x_{k, j}-1) + \sum_{k \ge j', k \neq j}p_kx_{k, j}$である。任意のジョブ$\lambda \neq j$について、$p_\lambda$が$R_{j', j}$にどう現れるか調べる:

• $\operatorname{pos}(\lambda) < \operatorname{pos}(j')$の場合
  • $\lambda < j'$なら$p_\lambda$の係数は$x_{j', k} + x_{k, j}-1=0$
  • $\lambda \ge j'$であることはない。なぜならジョブ$j'$は$r_{j'}$ちょうどに開始するため、$\operatorname{pos}(\lambda) < \operatorname{pos}(j')$なら$r_\lambda < r_{j'}$でなければならず、つまり$\lambda < j'$だから。
• $\operatorname{pos}(j') \le \operatorname{pos}(\lambda) < \operatorname{pos}(j)$の場合
  • $\lambda < j'$なら$p_\lambda$の係数は$x_{j', \lambda} + x_{\lambda, j} -1=1$
  • $\lambda \ge j'$なら$p_\lambda$の係数は$x_{\lambda, j}=1$
• $\operatorname{pos}(j) < \operatorname{pos}(\lambda)$の場合
  • $\lambda < j'$なら$p_\lambda$の係数は$x_{j', \lambda} + x_{\lambda, j} -1=0$
  • $\lambda \ge j'$なら$p_\lambda$の係数は$x_{\lambda, j}=0$

したがって、$R_{j', j} = r_{j'} + \sum_{\lambda: \operatorname{pos}(j') \le \operatorname{pos} (\lambda) < \operatorname{pos}(j)}p_\lambda$。

命題2の証明

$R_{j, j} = r_jx_{j, j} + \sum_{k < j}p_k(x_{j, k} + x_{k, j}-1) + \sum_{k > j}p_kx_{k, j}$である。

$k \neq j$なら常に $x_{j, k} + x_{k, j} - 1 = 0$ なので $\sum_{k < j}p_k(x_{j, k} + x_{k, j}-1) = 0$。次に$\sum_{k > j}p_kx_{k, j}$について考える。$k > j$ならば$r_k \ge r_j$だから、$j$を$r_j$ちょうどに開始するという仮定の下では$k$は必ず$j$の後に行うことになる。つまり、$x_{k, j} = 0$。従って総和もゼロであり、$R_{j, j} = r_j x_{j, j} = r_j$。

命題3の証明

$\operatorname{pos}(i) \le \operatorname{pos}(j)$の場合と$\operatorname{pos}(j) < \operatorname{pos}(i)$の場合にわけて検討する。

まず、$\operatorname{pos}(i) \le \operatorname{pos}(j)$の場合を扱う。任意のジョブ$\lambda \neq j$について、命題1の証明と同じように$p_\lambda$が$R_{i, j}$にどう現れるか調べる。

• $\operatorname{pos}(\lambda) < \operatorname{pos}(i)$の場合
  • $\lambda < i$なら$p_\lambda$の係数は $x_{i, \lambda}+x_{\lambda, j} - 1 = 0$
  • $\lambda \ge i$なら$p_\lambda$の係数は $x_{\lambda, j} = 1$ ――(A)
• $\operatorname{pos}(i) \le \operatorname{pos}(\lambda) < \operatorname{pos}(j)$の場合、命題1の証明と同様に$p_\lambda$の係数は1である
• $\operatorname{pos}(j) < \operatorname{pos}(\lambda)$の場合、命題1の証明と同様に$p_\lambda$の係数は0である

(A)に該当する$\lambda$が存在しなければ $R_{i, j} = r_i + \sum_{\lambda: \operatorname{pos}(i) \le \operatorname{pos} (\lambda) < \operatorname{pos}(j)}p_\lambda$ だから、次のことが成り立つ。

補題1: $\operatorname{pos}(i) \le \operatorname{pos}(j)$とする。(A)に該当するジョブが存在しなければ$R_{i, j} \le T_j$である。

また、(A)に該当する$\lambda$が存在する場合は実は次のことが成り立つ:

補題2: $\operatorname{pos}(i) \le \operatorname{pos}(j)$とし、(A)に該当するジョブが存在するとする。そのようなジョブの中でもっとも早く行うもの$\lambda$に注目する。言い換えると、$\lambda$は$\lambda \ge i$であるジョブの中でもっとも早く行うとする。このとき、$R_{i, j} \le R_{\lambda, j}$である。

証明: $\operatorname{pos}(i) \le \operatorname{pos}(j)$に注意すると、$R_{\lambda, j}, R_{i, j}$は以下のように表せる:

$$R_{\lambda, j} = r_\lambda + \sum_{k < \lambda, \operatorname{pos}(\lambda) < \operatorname{pos}(k) < \operatorname{pos}(j)}p_k + \sum_{k \ge \lambda, \operatorname{pos}(k) < \operatorname{pos}(j)}p_k$$ $$R_{i, j} = r_i + \sum_{k < i, \operatorname{pos}(i) < \operatorname{pos}(k) < \operatorname{pos}(j)}p_k + \sum_{k \ge i, \operatorname{pos}(k) < \operatorname{pos}(j)}p_k$$

$R_{i, j}$にどのような$p_k$が現れるか調べるために、総和に登場する添字$k$の集合を分割する。次の1～4のうち、1が一つ目の総和の$k$の集合で、2～4が二つ目の総和の$k$の集合の分割である。

1. $\{k \in [n]: k < i, \operatorname{pos}(i) < \operatorname{pos}(k)<\operatorname{pos}(j)\}$
2. $\{k \in [n]: i \le k < \lambda, \operatorname{pos}(k) < \operatorname{pos}(\lambda)\}$
3. $\{k \in [n] : i \le k < \lambda, \operatorname{pos}(\lambda) < \operatorname{pos}(k)<\operatorname{pos}(j)\}$
4. $\{k \in [n]: k \ge \lambda, \operatorname{pos}(k) < \operatorname{pos}(j)\}$

$i \le \lambda$である$\lambda$の中でもっとも早く行うものを選んでいるので2は空である。また、1と3の合併は$\{k \in [n]: k < \lambda, \operatorname{pos}(\lambda) < \operatorname{pos}(k) <\operatorname{pos}(j)\}$に包含され、後者は$R_{\lambda, j}$の一つ目の総和の添字集合に等しい。また、4は$R_{\lambda, j}$の二つ目の総和の添字集合に等しい。従って、$R_{i, j} \le R_{\lambda, j}$である。$\square$

次に、$\operatorname{pos}(j) < \operatorname{pos}(i)$の場合を検討する。$R_{i, j} = \sum_{k < i, k \neq j}p_k(x_{i, k} + x_{k, j}-1) + \sum_{k \ge i, k \neq j}p_kx_{k, j}$である。やはり任意のジョブ$\lambda \neq j$について、$p_\lambda$が$R_{i, j}$にどう現れるか調べる。

• $\operatorname{pos}(\lambda) < \operatorname{pos}(j)$の場合
  • $\lambda < i$ なら$p_\lambda$の係数は$x_{i, \lambda}+x_{\lambda, j} -1 = 0$
  • $\lambda \ge i$ なら$p_\lambda$の係数は$x_{\lambda, j} = 1$ ――(B)
• $\operatorname{pos}(j) < \operatorname{pos}(\lambda) < \operatorname{pos}(i)$の場合
  • $\lambda < i$なら$p_\lambda$の係数は$x_{i, \lambda}+x_{\lambda, j} -1 = -1$
  • $\lambda \ge i$なら$p_\lambda$の係数は$x_{\lambda, j} = 0$
• $\operatorname{pos}(i) \le \operatorname{pos}(\lambda)$の場合
  • $\lambda < i$なら$p_\lambda$の係数は$x_{i, \lambda} + x_{\lambda, j} -1=0$
  • $\lambda \ge i$なら$p_\lambda$の係数は$x_{\lambda, j}=0$

(B)に該当する$\lambda$が存在しなければ、$R_{i, j} \le 0 \le T_j$だから、次のことが成り立つ:

補題3: $\operatorname{pos}(j) < \operatorname{pos}(i)$とする。(B)に該当する$\lambda$が存在しなければ$R_{i, j} \le T_j$である。

また、(B)に該当する$\lambda$が存在する場合は次のことが成り立つ:

補題4: $\operatorname{pos}(j) < \operatorname{pos}(i)$とする。(B)に該当するジョブが存在するとする。そのようなジョブの中でもっとも早く行うもの$\lambda$に注目する。言い換えると、$\lambda$は$\lambda \ge i$であるジョブの中でもっとも早く行うとする。このとき、$R_{i, j} \le R_{\lambda, j}$である。

証明: $\operatorname{pos}(\lambda) < \operatorname{pos}(j) < \operatorname{pos}(i)$に注意すると、$R_{\lambda, j}, R_{i, j}$は以下のように評価できる:

$$R_{\lambda, j} = \sum_{k < \lambda, \operatorname{pos}(\lambda) < \operatorname{pos}(k) < \operatorname{pos}(j)}p_k + \sum_{k \ge \lambda, \operatorname{pos}(k) < \operatorname{pos}(j)}p_k$$ $$R_{i, j} = -\sum_{k < i, \operatorname{pos}(j) < \operatorname{pos}(k) < \operatorname{pos}(i)}p_k + \sum_{k \ge i, \operatorname{pos}(k) < \operatorname{pos}(j)}p_k \le \sum_{k \ge i, \operatorname{pos}(k) < \operatorname{pos}(j)}p_k$$

$R_{i, j}$の上界$\sum_{k \ge i, \operatorname{pos}(k) < \operatorname{pos}(j)}p_k$にどのような$p_k$が現れるか調べるために、添字$k$の集合を3つに分割する:

1. $\{k \in [n]: i \le k < \lambda, \operatorname{pos}(k) < \operatorname{pos}(\lambda)\}$
2. $\{k \in [n] : i \le k < \lambda, \operatorname{pos}(\lambda) < \operatorname{pos}(k)<\operatorname{pos}(j)\}$
3. $\{k \in [n]: k \ge \lambda, \operatorname{pos}(k) < \operatorname{pos}(j)\}$

$i \le \lambda$であるジョブ$\lambda$の中でもっとも早く行うものを選んでいるので1は空である。また、2は$R_{\lambda, j}$の一つ目の総和の添字集合に包含され、3は二つ目の総和の添字集合に等しい。従って、$R_{i, j} \le R_{\lambda, j}$である。$\square$

補題1～4を適用することで、命題3が示せる: 補題1と3より、(A)や(B)に該当するような$\lambda$が存在しない$(i, j)$に対しては$R_{i, j} \le T_j$である。また、(A)や(B)に該当するような$\lambda$が存在する場合、どちらの場合も$i < \lambda$だから、補題2または4を繰り返し適用することで、$i < \lambda_1 < \lambda_2 < ...$, $R_{i, j} \le R_{\lambda_1, j} \le R_{\lambda_2, j} \le ...$ という列が得られる。高々有限回の適用で(A)や(B)に該当する項がないものが見つかり、補題1または3によって$R_{i, j} \le T_j$が従う。

これで$\operatorname{Form*}$の正当性が示された。

---

上のNemhauserとSavelsberghの論文ではLP双対を経由した非常に短い「証明」が書かれている……が、読んでも双対を取ることで問題が自明になっていると思えなかった。双対に注目しても、ギャップを埋めようとすれば$R_{i, j}$のうちどれがどれより大きいのかという細かい議論を避けられないように思える。

実験

big Mを避けることで実際に速くなるかどうか知りたい。$\operatorname{FormM}$と$\operatorname{Form*}$をSCIPに与えて、ランダムインスタンスで比較した。

import math
import random
import statistics
from dataclasses import dataclass
from itertools import pairwise, product
from typing import Any

import pandas as pd
import pulp


@dataclass
class SchedulingInstance:
    n: int
    release_times: list[float]
    weights: list[float]
    processing_times: list[float]


@dataclass
class SchedulingSolution:
    n: int
    start_times: list[float]
    model: pulp.LpProblem  # artifact


def validate_solution(
    instance: SchedulingInstance, solution: SchedulingSolution
) -> None:
    n = instance.n
    release_times = instance.release_times
    processing_times = instance.processing_times
    start_times = solution.start_times
    eps = 1e-8

    for i in range(n):
        assert release_times[i] <= start_times[i] + eps

    for i1, i2 in pairwise(sorted(range(n), key=lambda i: start_times[i])):
        assert start_times[i1] + processing_times[i1] <= start_times[i2] + eps
        # optimality condition
        assert math.isclose(
            max(start_times[i1] + processing_times[i1], release_times[i2]),
            start_times[i2],
            rel_tol=eps,
            abs_tol=eps,
        )


def make_random_instance(n: int, mag: float, seed: int) -> SchedulingInstance:
    rng = random.Random(seed)
    release_times = sorted((rng.randrange(0, mag) for _ in range(n)))
    weights = [rng.randrange(1, mag) for _ in range(n)]
    processing_times = [rng.randrange(1, mag) for _ in range(n)]
    return SchedulingInstance(n, release_times, weights, processing_times)


def solve(
    instance: SchedulingInstance, is_big_m: bool, time_out_sec: float = 600.0
) -> SchedulingSolution | None:
    n = instance.n
    release_times = instance.release_times
    weights = instance.weights
    processing_times = instance.processing_times

    model = pulp.LpProblem("single_machine_scheduling", pulp.LpMinimize)

    start_times = [pulp.LpVariable(f"start_time_{i}") for i in range(n)]
    is_priors = [
        [pulp.LpVariable(f"is_prior_{i}_{j}", cat=pulp.LpBinary) for i in range(n)]
        for j in range(n)
    ]

    model += pulp.lpSum(start_times[i] * weights[i] for i in range(n))

    def is_prior(i: int, j: int) -> pulp.LpAffineExpression:
        assert i != j
        if i > j:
            return pulp.LpAffineExpression(1 - is_priors[j][i])
        else:
            return pulp.LpAffineExpression(is_priors[i][j])

    for i, j, k in product(range(n), range(n), range(n)):
        if i == j or j == k or k == i:
            continue
        model += is_prior(i, j) + is_prior(j, k) <= 1 + is_prior(i, k)
    for i in range(n):
        model += start_times[i] >= release_times[i]

    if is_big_m:
        big_m = 2 * (release_times[-1] + sum(processing_times)) + 1
        for i, j in product(range(n), range(n)):
            if i == j:
                continue
            model += start_times[j] >= start_times[i] + processing_times[
                i
            ] - big_m * is_prior(j, i)
    else:
        for i, j in product(range(n), range(n)):
            if i == j:
                continue
            model += start_times[j] >= release_times[i] * is_prior(i, j) + pulp.lpSum(
                processing_times[k] * (is_prior(i, k) + is_prior(k, j) - 1)
                for k in range(i)
                if k != j
            ) + pulp.lpSum(
                processing_times[k] * is_prior(k, j) for k in range(i, n) if k != j
            )

    solver = pulp.SCIP_CMD(timeLimit=time_out_sec)
    model.solve(solver)
    print(f"{model.solutionTime=}")
    if model.status == pulp.LpStatusNotSolved:
        return None
    assert model.status == pulp.LpStatusOptimal, pulp.LpStatus[model.status]
    solution_start_times = [pulp.value(start_times[i]) for i in range(n)]
    sol = SchedulingSolution(n, solution_start_times, model)
    validate_solution(instance, sol)
    return sol


def time_with_param(n: int, mag: float, is_big_m: bool, n_iter: int) -> dict[str, Any]:
    times: list[float] = []
    objs: list[float] = []
    for seed in range(n_iter):
        instance = make_random_instance(n, mag, seed)
        sol = solve(instance, is_big_m, time_out_sec=600.0)
        times.append(sol.model.solutionTime if sol else math.inf)
        objs.append(pulp.value(sol.model.objective) if sol else math.nan)
    print(objs)
    return (
        {"n": n, "method": "big M" if is_big_m else "no M"}
        | {f"seed {i}": times[i] for i in range(n_iter)}
        | {"median": statistics.median(times)}
    )


df = pd.DataFrame(
    [
        time_with_param(10, 10, True, 5),
        time_with_param(15, 15, True, 5),
        time_with_param(10, 10, False, 5),
        time_with_param(15, 15, False, 5),
        time_with_param(20, 20, False, 5),
        time_with_param(30, 30, False, 5),
        time_with_param(50, 50, False, 5),
    ],
)
print(df.to_markdown(index=False, floatfmt=".2f"))

結果は次のようになった。infは600秒以内に解けなかったことを表す。

n	method	seed 0	seed 1	seed 2	seed 3	seed 4	median
10	big M	2.26	0.94	3.95	1.30	9.83	2.26
15	big M	inf	298.36	inf	233.59	inf	inf
10	no M	0.21	3.47	2.25	1.20	0.07	1.20
15	no M	4.04	4.07	2.56	0.70	4.65	4.04
20	no M	3.29	0.32	3.60	0.36	3.41	3.29
30	no M	12.25	4.03	5.63	5.40	14.96	5.63
50	no M	20.19	65.13	inf	inf	37.51	65.13

big Mを避けた定式化のほうが速く解けると言える。しかし、この定式化であっても、$O(2^nn)$の自明なDPで扱える規模($n < 30$)を超えてSCIPで新たに解けるようになる範囲は+10～20程度のようだ。ただし、上の定式化は最低限の制約しか使っておらず、冗長な制約（妥当不等式）を追加していくなどによってさらに速くなる可能性はある。妥当不等式に関しては、NemhauserとSavelsberghの論文にも議論がある。

その他

NemhauserとSavelsberghのテクニックは、機械が1つであれば冒頭の問題設定以外にも様々なリリース時刻付き問題に使える。例えば、ジョブに依存関係がある$1|\textrm{prec};r_j|\sum w_jC_j$や、納期に対する遅れ時間$\max_i(T_i - d_i)$を最小化する$1|r_j|L_{\max}$など。

Khowalaが単一マシンの様々なジョブスケジューリング問題と定式化のベンチマークを行っていて、ここで扱ったような処理順序を表すバイナリ変数を導入するやり方以外の定式化も知ることができる。

この問題に関しては、ジョブの順番に制約がなくて順番を決めれば最適解が自明に求まるという性質上、局所探索が明らかに簡単かつ有効ではある。実用的にはサイズがかなり大きくなっても困らないタイプの問題と言えるかもしれない。（実際はだいたい面倒な制約が付くだろうけど……）

文献

• Brucker, Peter, Jan Karel Lenstra and Alexander H. G. Rinnooy Kan. “Complexity of machine scheduling problems.” (1975).
• Khowala, Ketan. “Single Machine Scheduling: Comparison of MIP Formulations and Heuristics for Interfering Job Sets.” (2012).
• Nemhauser, George L. and Mwp Martin Savelsbergh. “A Cutting Plane Algorithm for the Single Machine Scheduling Problem with Release Times.” (1992).

functional graphの閉路の情報を得る

問題: 関数$f: \{0, 1, ..., N-1\} \rightarrow \{0, 1, ..., N-1\}$が与えられる。$f$によって定まるfunctional graph上の閉路の情報を$O(N)$で得たい。具体的には、次の配列を作りたい:

• $\operatorname {cycle}[v] :=$ 閉路上の頂点で$v$に最も近いもの。より形式的には、$f^k(v)$が閉路に含まれるような最小の非負整数$k$に関する$f^k(v)$。
• $\operatorname{length}[v] := v$が閉路に含まれるならその閉路の長さ。そうでないなら$v$から$\operatorname{cycle}[v]$までの距離。

以下はグラフと求めたい値の例。頂点番号 (cycleの値, lengthの値)という表記になっている。

アルゴリズム

長さ$N$の配列$\mathrm{tmp}, \mathrm{length}, \mathrm{cycle}$を初期化する。$\mathrm{tmp}$は$-1$で埋め、$\mathrm{length}, \mathrm{cycle}$は$0$で埋める。

頂点$r = 0, 1, ..., N-1$について、それぞれ未訪問($\mathrm{tmp}[r] = -1$)であれば$r$に関するイテレーションを開始する。

1回のイテレーションにおいては、$v \leftarrow r$を始点にして$v \leftarrow f(v)$を繰り返してグラフを辿っていく。この際、$r$から$v$までの距離$d$を保持して、$\mathrm{tmp}[v] \leftarrow d$と記録していく。以前に訪れた頂点($\mathrm{tmp}[v] \neq -1$)を再訪したところでストップするが、このとき、$v$がこのイテレーションで初めて訪問した頂点である($\mathrm{length}[v] = 0$)か、以前のイテレーションで訪問した頂点かで必要な処理が異なる。

$v$がこのイテレーションで初めて訪問した頂点であれば新しい閉路を発見したことになるので、閉路の長さ($d - \mathrm{tmp}[v]$)に従って、閉路をもう一周しながら$\mathrm{length}$と$\mathrm{cycle}$を埋める。さらに、始点$r$から閉路に達するまでの頂点についても$\mathrm{length}$と$\mathrm{cycle}$を埋める。

$v$が以前のイテレーションで訪問した頂点の場合、始点$r$から辿って$\mathrm{cycle}$を$\mathrm{cycle}[v]$で埋めていく。$\mathrm{length}$については、$v$が閉路中の頂点($v = \mathrm{cycle}[v]$)であれば$\mathrm{tmp}[v]$を採用し、$v$が閉路外の頂点であれば$\mathrm{tmp}[v]+\mathrm{length}[v]$を採用して始点から埋めていく。

どの頂点についても高々2回しか$f$を適用しないので$O(N)$。

以下はPythonの実装例:

def get_cycle_info(f: list[int]) -> tuple[list[int], list[int]]:
    """{0, 1, ..., N-1}上の写像fを入力として2つのリストcycle, lengthを返す

    cycle[v] := 閉路に含まれるような最初のf(f(...f(v)...))
    length[v] := vが閉路上の頂点の場合は閉路の長さ、そうでない場合はvから閉路までの距離
    """
    n = len(f)
    tmp = [-1] * n
    cycle = [0] * n
    length = [0] * n
    for start in range(n):
        v = start
        dist = 0
        while tmp[v] == -1:
            tmp[v] = dist
            v = f[v]
            dist += 1
        if length[v] == 0:
            u = start
            path_length = tmp[v]
            for i in range(path_length):
                cycle[u] = v
                length[u] = path_length - i
                u = f[u]
            cycle_length = dist - tmp[v]
            for _ in range(cycle_length):
                cycle[v] = v
                length[v] = cycle_length
                v = f[v]
        else:
            u = start
            path_length = dist if v == cycle[v] else dist + length[v]
            for i in range(dist):
                cycle[u] = cycle[v]
                length[u] = path_length - i
                u = f[u]
    return cycle, length

cycle, length = get_cycle_info([1, 2, 3, 4, 5, 6, 3, 8, 1, 10, 11, 10, 6, 12, 14])
assert cycle == [3, 3, 3, 3, 4, 5, 6, 3, 3, 10, 10, 11, 6, 6, 14]
assert length == [3, 2, 1, 4, 4, 4, 4, 4, 3, 1, 2, 2, 1, 2, 1]

cycle, length = get_cycle_info([1, 2, 2])
assert cycle == [2, 2, 2]
assert length == [2, 1, 1]

cycle, length = get_cycle_info([])
assert cycle == []
assert length == []

たぶん基本的な問題だけど書いてある文献が見つからなかったのでメモ。