2025年1月25日土曜日

Pythonでソルバから双対変数を取り出す方法とその落とし穴

以下の最適化問題が与えられたとする。

maximizef(x)subject togi(x)0for i=1,...,kxRn \begin{array}{lll} \text{maximize} & f(x)\\ \text{subject to} & g_i(x) \leq 0 & \text{for } i = 1, ..., k\\ & x \in \mathbb R^n\\ \end{array}

このうちのひとつの制約gj(x)0g_j(x)\le 0を選んで、他の制約を固定しつつこの制約の右辺の定数を微小に増減させた時に目的関数値がどう変わるか知りたいことがある。これを制約gj(x)0g_j(x) \le 0のシャドウプライスと呼ぶ。つまり、

Lj(b):=maximizef(x)subject togi(x)0for i=1,...,j1,j+1,...,kgj(x)bxRn \begin{array}{llll} L_j(b) :=& \text{maximize} & f(x)\\ & \text{subject to} & g_i(x) \leq 0 & \text{for } i = 1, ..., j-1, j+1, ..., k\\ && g_j(x) \le b \\ && x \in \mathbb R^n\\ \end{array}

としたときのLjL_jb=0b=0における傾きがシャドウプライスである。[1] 微分可能な場合は微分だし、そうでない場合も劣微分を考えればよい。ちなみにLPの場合にはLjL_jは区分線形凹関数になることが知られている。[2] これは不等式制約の例だが、等式制約についても同様に定義できる。

LPや一部の凸最適化問題など、ラグランジュ双対問題との強双対性が成り立つものについては、シャドウプライスは制約に対応する双対変数の値に一致する。[3] たいていのソルバは求解の過程で双対問題も同時に解いているので問題が解ければこの値も簡単に得られる。例えばSoPlexをCLIで起動する場合はsoplex <input file> -yなどで参照できる。

Pythonのソルバインターフェースやモデリングツールも普通は双対変数を取り出す機能を備えている。例えばPuLPだったら制約オブジェクトのpiというアトリビュートにシャドウプライスが格納される[4]

import pulp as pl

problem = pl.LpProblem("problem", pl.LpMaximize)
x = pl.LpVariable("x", lowBound=0)
y = pl.LpVariable("y", lowBound=0)
constr = x + 2 * y <= 1
problem += constr
problem += 2 * x + y
problem.solve()
print(constr.pi)

上の例を実行すると、x+2y1x + 2y \le 1のシャドウプライスが2.0であると正しく表示される。

しかし、「符号の正しいシャドウプライス」を得るためには注意が必要である。一般にシャドウプライスは制約の右辺の増分に対して定義される。これはx+2y1x+2y \le 11x+2y1 \ge x + 2yではシャドウプライスの符号が逆になると解釈できる。もっとも、シャドウプライスは普通はf(x)0f(x) \le 0とかAxbAx \ge bのように制約のある種の「標準形」に対して定義されるので、非標準的な立式に対して符号がどちらになるべきかは曖昧と言えるかもしれない。ただ、少なくとも使っているツールのシャドウプライスの符号に一貫的でわかりやすいルールがあることは期待したくなる。しかし、現実にはそれはかなわないようだ。

import pulp as pl

x = pl.LpVariable("x", lowBound=0)
y = pl.LpVariable("y", lowBound=0)


def solve(constr: pl.LpConstraint) -> float:
    problem = pl.LpProblem("problem", pl.LpMaximize)
    problem += constr
    problem += 2 * x + y

    problem.solve()
    return constr.pi


print(solve(x + 2 * y <= 1))
print(solve(-x -2 * y >= -1))
print(solve(1 >= x + 2 * y))
print(solve(1 - x >= 2 * y))

上はx+2y1x+2y \le 1を様々な形で表した場合にどうなるかを調べていて、結果は下のようになる:

2.0  :  x + 2y <= 1
-2.0 : -x - 2y >= -1
2.0  :       1 >= x + 2y
-2.0 :   1 - x >= 2y

シャドウプライスを右辺の増分に対して定義されると考える場合、上から2.0, -2.0, -2.0, -2.0になるはずなので、3つ目だけ符号が逆ということになる。これだけが異なるというのは自然な説明が難しく、驚きを小さくするという原則に反しているようにも思える。PuLPのLpAffineExpressionの実装を見るとこうなっている理由がわかる:

class LpAffineExpression(_DICT_TYPE):
    ...
    def __le__(self, other):
        return LpConstraint(self - other, const.LpConstraintLE)

    def __ge__(self, other):
        return LpConstraint(self - other, const.LpConstraintGE)

__le__<=の実装、__ge__>=の実装である。a <= ba - b <= 0と等価であると読める。これはシャドウプライスの観点でも等価である。しかし、aがfloatの場合、Pythonは最初にfloatの__le__を調べ、floatにはLpAffineExpressionとの比較が実装されていないので次に対称律に従ってb__ge__を試すというステップを踏む。その結果b - a >= 0と解釈されてしまう。上の例でも、1 >= x + 2yの代わりにpl.LpAffineExpression(1) >= x + 2yとあらかじめ変換しておけばシャドウプライスが意図通りの符号になる。

つまりこの問題は、Pythonが比較演算子を対称律の成り立つ述語として想定しているところを、述語でもない上にa >= bb <= aで意味が異なりうるような使い方をしたことが原因で生じていると言え、解決が難しい。実際、OR-ToolsやPyomoなど他のモデリングツールも同じ問題を抱えているようだ。

OR-Tools:

from collections.abc import Callable
from ortools.linear_solver.pywraplp import Solver, Constraint, Variable


def solve(constr_fn: Callable[[Variable, Variable], Constraint]) -> float:
    solver = Solver.CreateSolver("GLOP")
    x = solver.NumVar(0, Solver.Infinity(), "x")
    y = solver.NumVar(0, Solver.Infinity(), "y")
    solver.Maximize(2 * x + y)
    constr = constr_fn(x, y)
    solver.Add(constr, "constr")
    solver.Solve()
    solver.ComputeConstraintActivities()
    return next(c.dual_value() for c in solver.constraints() if c.name() == "constr")


print(solve(lambda x, y: x + 2 * y <= 1))
print(solve(lambda x, y: -x - 2 * y >= -1))
print(solve(lambda x, y: 1 >= x + 2 * y))
print(solve(lambda x, y: 1 - x >= 2 * y))

# 2.0
# -2.0
# 2.0
# -2.0

Pyomo:

# Pyomo
from collections.abc import Callable
import pyomo.environ as pyo


def solve(constr_fn: Callable[[pyo.ScalarVar, pyo.ScalarVar], pyo.Constraint]) -> float:
    model = pyo.ConcreteModel()
    model.x = pyo.ScalarVar(bounds=(0, None))
    model.y = pyo.Var(bounds=(0, None))
    model.obj = pyo.Objective(expr=2 * model.x + model.y, sense=pyo.maximize)
    model.constr = pyo.Constraint(expr=constr_fn(model.x, model.y), name="constr")
    model.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)
    pyo.SolverFactory("glpk").solve(model)
    return model.dual[model.constr]


print(solve(lambda x, y: x + 2 * y <= 1))
print(solve(lambda x, y: -x - 2 * y >= -1))
print(solve(lambda x, y: 1 >= x + 2 * y))
print(solve(lambda x, y: 1 - x >= 2 * y))

# 2.0
# -2.0
# 2.0
# 2.0

OR-ToolsはPuLPと同じだが、Pyomoはさらに1x2y1-x \ge 2yの場合も予想と違っている。

結局、意図通りの符号のシャドウプライスが欲しい場合は使っているツールの実装の詳細に気を配る必要があるということになる。上で見たようにPythonの基本的な挙動に関わる問題なので、完全に「自然」に動くツールがあるとは思えない。

ちなみにAMPLでも試すと、2, -2, -2, -2となって意図通りになる:

option solver 'gurobi';

var x >= 0;
var y >= 0;

maximize obj: 2*x + y;

subject to c1: x + 2*y <= 1;

solve;
display c1.dual;
drop c1;

subject to c2: -x - 2*y >= -1;
solve;
display c2.dual;
drop c2;

subject to c3: 1 >= x + 2*y;
solve;
display c3.dual;
drop c3;

subject to c4: 1 - x >= 2*y;
solve;
display c4.dual;
drop c4;

  1. シャドウプライスはWikipediaでは経済学の言葉として説明されている。感度 (sensitivity) という用語もあるが、これは言葉の指し示す範囲が広く、制約や目的関数の係数を増減させた場合の効果もそう呼ばれていると認識している。最適化系の文献ではシャドウプライスは感度の一種として特別な名前がついていないこともある。また、後述するように、ラグランジュ双対問題との強双対性がある場合はシャドウプライスは双対変数の値に一致するので、dual multiplierとかLagrange multiplierといった概念と同一視されていることも多いと思う。 ↩︎

  2. 例えばMartin, R. Kipp. “Large scale linear and integer optimization.” のCorollary 2.49 ↩︎

  3. Boyd, Stephen P. and Lieven Vandenberghe. “Convex Optimization.” の 5.6. Perturbation and sensitivity analysis などを参照 ↩︎

  4. ただし、ドキュメントには説明がない。関連するPR (#202, #781) やイシューを眺める限り、半公式の機能とみなされているとは思う。 ↩︎

2024年2月23日金曜日

SCIPにLPの基底解を出させるのは難しい

SCIPに整数性のあるLPを解かせてみる。公式のPySCIPOptを使う。

from itertools import combinations
from pyscipopt import Model, quicksum, SCIP_PARAMSETTING

n = 4
model = Model("graph_matching")
xs = {
    (i, j): model.addVar(lb=0.0, name=f"x_{i}_{j}")
    for i, j in combinations(range(n), 2)
}
for v in range(n):
    model.addCons(quicksum(x for (i, j), x in xs.items() if i == v or j == v) == 1)
model.setObjective(quicksum(xs.values()), "minimize")
model.optimize()
sol = {(i, j): model.getVal(x) for (i, j), x in xs.items()}
print(f"objective: {model.getObjVal()}")
print(f"solution: {sol}")

制約はn=4n=4の完全グラフに関する完全マッチング多面体を表していて、変数が6つ、等号制約が4つなので、基底解はゼロになる変数を少なくとも2つ含む。例えば(1,0,0,0,0,1)(1, 0, 0, 0, 0, 1)など。しかし、SCIPは基底解を返さない:

objective: 2.0
solution: {(0, 1): 0.3333333333333333, (0, 2): 0.3333333333333333, (0, 3): 0.3333333333333333, (1, 2): 0.33333333333333337, (1, 3): 0.33333333333333337, (2, 3): 0.33333333333333337}

SCIPはLPソルバではないので、分枝限定法の1回目のLPを行う前に問題を変形したりカットを追加したりはしているはずで、この結果はおかしくはない。ただ、READMEを眺めてpresolveをオフにしてみたり、他の設定を弄ったりしても特に結果に変わりなかった:

model.setPresolve(SCIP_PARAMSETTING.OFF)
model.setHeuristics(SCIP_PARAMSETTING.OFF)
model.disablePropagation()

model.writeProblem("graph_matching.lp")でLPを書き出して直接SCIPに解かせてみることにする。LPファイルの中身は以下のようになる:

\ SCIP STATISTICS
\   Problem name     : graph_matching
\   Variables        : 6 (0 binary, 0 integer, 0 implicit integer, 6 continuous)
\   Constraints      : 4
Minimize
 Obj: +1 x01 +1 x02 +1 x03 +1 x12 +1 x13 +1 x23
Subject to
 c1: +1 x01 +1 x02 +1 x03 = +1
 c2: +1 x01 +1 x12 +1 x13 = +1
 c3: +1 x02 +1 x12 +1 x23 = +1
 c4: +1 x03 +1 x13 +1 x23 = +1
Bounds
End

これをSCIPのシェルを直接使って解いてもやはり結果は変わらない:

SCIP> read graph_matching.lp
SCIP> set presolving emphasis off
SCIP> set heuristics emphasis off
SCIP> set separating emphasis off
SCIP> optimize
SCIP> display solution

objective value:                                    2
x01                                 0.333333333333333   (obj:1)
x02                                 0.333333333333333   (obj:1)
x03                                 0.333333333333333   (obj:1)
x12                                 0.333333333333333   (obj:1)
x13                                 0.333333333333333   (obj:1)
x23                                 0.333333333333333   (obj:1)

write transproblemでpresolve後のLPを出力してみる:

\ SCIP STATISTICS
\   Problem name     : t_graph_matching.lp
\   Variables        : 6 (0 binary, 0 integer, 0 implicit integer, 6 continuous)
\   Constraints      : 12
Minimize
 Obj: +1 t_x01 +1 t_x02 +1 t_x03 +1 t_x12 +1 t_x13 +1 t_x23
Subject to
 c1: +1 t_x01 +1 t_x02 +1 t_x03 = +1
 c2: +1 t_x01 +1 t_x12 +1 t_x13 = +1
 c3: +1 t_x02 +1 t_x12 +1 t_x23 = +1
 c4: +1 t_x03 +1 t_x13 +1 t_x23 = +1
 SSTcut_0_1: -1 t_x01 +1 t_x02 <= +0
 SSTcut_0_2: -1 t_x01 +1 t_x03 <= +0
 SSTcut_0_3: -1 t_x01 +1 t_x12 <= +0
 SSTcut_0_4: -1 t_x01 +1 t_x13 <= +0
 SSTcut_0_5: -1 t_x01 +1 t_x23 <= +0
 SSTcut_1_2: -1 t_x02 +1 t_x03 <= +0
 SSTcut_1_3: -1 t_x02 +1 t_x12 <= +0
 SSTcut_1_4: -1 t_x02 +1 t_x13 <= +0
Bounds
 0 <= t_x01 <= 1
 0 <= t_x02 <= 1
 0 <= t_x03 <= 1
 0 <= t_x12 <= 1
 0 <= t_x13 <= 1
 0 <= t_x23 <= 1
End

SSTcutが犯人で、x01x02x_{01} \ge x_{02}のような制約を等号にすると確かに上のような「基底解」が出ることになる。[1] 問題はset presolve emphasis offでもこのカットが無効にならないことだが、とりあえずset misc usesymmetry 0とすることで解決した:

SCIP> display solution

objective value:                                    2
x03                                                 1   (obj:1)
x12                                                 1   (obj:1)

後から気づいたが、PySCIPOptからも同等の設定がmodel.setIntParam("misc/usesymmetry", 0)でできるようだ。

結論として、set presolve emphasis offやそれに対応するPySCIPOptのsetPresolve(SCIP_PARAMSETTING.OFF)はpresolveのフェーズで起こることを完全に無しにしてくれるわけではない。おそらくmisc/usesymmetry = 0を加えても完全ではなく、stackoverflowでも別の人が別の問題で困っているようだった。

そもそもLPソルバが使いたいならSCIPの背後にあるSoPlexを直接使うべきなんだろうけど、SoPlexには公式のPythonバインディングがないという問題がある。[2] Pythonで単体法のLPソルバが使いたい時はSCIPではなく、CBC/CLPやHiGHSを検討したほうがいいかもしれない。特にCBC/CLPのPythonバインディングであるCyLPは、READMEを見た感じだとLPソルバを使って自前でILPを解くような状況を意識しているようだ。(試してない。)


  1. 問題の対称性がある時に等価な枝を調べるのを避けるためのカット ↩︎

  2. 試みはある: https://github.com/SBRG/soplex_cython ↩︎

2024年1月27日土曜日

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

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

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

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

  • 任意のジョブjjについて tjrjt_j \ge r_j
  • 任意の異なるジョブi,ji, jについて tjti+piMxj,it_j \ge t_i + p_i - Mx_{j, i}
    • MMは十分に大きな定数
  • 任意の異なるジョブi,ji, jについて xi,j+xj,i=1x_{i, j} + x_{j, i} = 1
  • 任意の異なるジョブi,j,ki, j, kについて xi,j+xj,k1+xi,kx_{i, j} + x_{j, k} \le 1 + x_{i, k}
  • 任意の異なるジョブi,ji, jについて xi,j{0,1}x_{i, j} \in \{0, 1\}
  • 任意のジョブjjについて tjRt_j \in \mathbb R

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

tjrixi,j+k<i,kjpk(xi,k+xk,j1)+ki,kjpkxk,j(*) \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を使わずに書けるというのが論文の内容である。ただし、便宜上xi,i=1x_{i,i}=1とする。上のbig Mを使った定式化をFormM\operatorname{FormM}()(*)を使った定式化をForm*\operatorname{Form*}とする。Form*\operatorname{Form*}の正当性を確かめたい。

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

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

これが示されれば、

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

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

定理の証明

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

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

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

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

である。したがって、次の3つの命題が示せればよい。以下、()(*)の右辺を Ri,j:=rixi,j+k<i,kjpk(xi,k+xk,j1)+ki,kjpkxk,jR_{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} とする。ただし、ここではxi,jx_{i, j}は与えられた順列π\piに応じて値が定まっている定数である。

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

命題2: 任意のジョブjjについて、Tj=rjT_j = r_jならばRj,j=rjR_{j, j} = r_jである。

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

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

命題1の証明

Rj,j=rj+k<j,kjpk(xj,k+xk,j1)+kj,kjpkxk,jR_{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}である。任意のジョブλj\lambda \neq jについて、pλp_\lambdaRj,jR_{j', j}にどう現れるか調べる:

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

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

命題2の証明

Rj,j=rjxj,j+k<jpk(xj,k+xk,j1)+k>jpkxk,jR_{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}である。

kjk \neq jなら常に xj,k+xk,j1=0x_{j, k} + x_{k, j} - 1 = 0 なので k<jpk(xj,k+xk,j1)=0\sum_{k < j}p_k(x_{j, k} + x_{k, j}-1) = 0。次にk>jpkxk,j\sum_{k > j}p_kx_{k, j}について考える。k>jk > jならばrkrjr_k \ge r_jだから、jjrjr_jちょうどに開始するという仮定の下ではkkは必ずjjの後に行うことになる。つまり、xk,j=0x_{k, j} = 0。従って総和もゼロであり、Rj,j=rjxj,j=rjR_{j, j} = r_j x_{j, j} = r_j

命題3の証明

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

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

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

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

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

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

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

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

Rλ,j=rλ+k<λ,pos(λ)<pos(k)<pos(j)pk+kλ,pos(k)<pos(j)pk 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 Ri,j=ri+k<i,pos(i)<pos(k)<pos(j)pk+ki,pos(k)<pos(j)pk 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

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

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

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

次に、pos(j)<pos(i)\operatorname{pos}(j) < \operatorname{pos}(i)の場合を検討する。Ri,j=k<i,kjpk(xi,k+xk,j1)+ki,kjpkxk,jR_{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}である。やはり任意のジョブλj\lambda \neq jについて、pλp_\lambdaRi,jR_{i, j}にどう現れるか調べる。

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

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

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

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

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

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

Rλ,j=k<λ,pos(λ)<pos(k)<pos(j)pk+kλ,pos(k)<pos(j)pk 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 Ri,j=k<i,pos(j)<pos(k)<pos(i)pk+ki,pos(k)<pos(j)pkki,pos(k)<pos(j)pk 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

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

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

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

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

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


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

実験

big Mを避けることで実際に速くなるかどうか知りたい。FormM\operatorname{FormM}Form*\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(2nn)O(2^nn)の自明なDPで扱える規模(n<30n < 30)を超えてSCIPで新たに解けるようになる範囲は+10~20程度のようだ。ただし、上の定式化は最低限の制約しか使っておらず、冗長な制約(妥当不等式)を追加していくなどによってさらに速くなる可能性はある。妥当不等式に関しては、NemhauserとSavelsberghの論文にも議論がある。

その他

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

2024年1月18日木曜日

functional graphの閉路の情報を得る

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

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

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

アルゴリズム

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

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

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

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

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

どの頂点についても高々2回しかffを適用しないのでO(N)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 == []

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

2021年11月16日火曜日

ABC 226 F - Score of Permutations

F(M)=F(M) =長さNNの順列で、すべてのサイクルの長さがMMの約数であるようなものの数 f(m)=f(m)=長さNNの順列で、サイクルの長さのLCMがmmであるようなものの数

とするとき、求める数はNNの分割のLCMとして得られるようなすべてのmmについてf(m)mKf(m)m^Kの総和を取ったものに等しい。

F(M)=mMf(m)F(M) = \sum_{m | M} f(m)であり、ffFFのメビウス反転として計算できる。また、NNの分割のLCMの最大値はランダウ関数g(N)g(N)として知られていて、OEISによるとg(50)=180180g(50) = 180180なので、個々のF(M)F(M)を高速に計算できればよさそう。

F(M)F(M)は単純なDPで求まる。dp[x]=1,2,...,N\operatorname{dp}[x] = 1, 2, ..., Nからxx個使って任意個数の長さがMMの約数であるようなサイクルを作る場合の数、とするとき、MMの各約数m(N)m ( \le N)について

dp[x+m]+=(m1)!(Nx1m1)dp[x] \operatorname{dp}[x+m] \mathrm{{+}{=}} (m-1)! \binom{N-x-1}{m-1}\operatorname{dp}[x]

と遷移すればよい。(m1)!(m-1)!は長さmmのサイクルの総数であり、(Nx1m1)\binom{N-x-1}{m-1}はまだ使っていないNxN-x個の数から(完全な数え上げをするために)必ず最小の数は選ぶとして残りのm1m-1個を自由に選ぶことに対応している。計算量はO((N2+logK)g(N)){O}((N^2 +\log K)g(N))[1]

また、FPSで考えることもできる。MMの約数だけに限って各長さのサイクルの総数を数えたEGFをDM(x)\mathcal D_M(x)とするとき、DM(x)=mM(m1)!xm/m!=mMxm/m\mathcal D_M(x)=\sum_{m | M}(m-1)!x^m/m! = \sum_{m|M}x^m/mである。このとき、F(M)=[xN/N!]exp(DM(x))F(M) = [x^N/N!]\exp(\mathcal D_M(x))である。[2] FFTを使うと計算量はO((NlogN+logK)g(N)){O}((N\log N + \log K) g(N))


解説を見ると、LCMごとにまとめなくてもすべての分割を試せるらしい。確かに。


  1. 実際には1,2,...,g(N)1, 2, ..., g(N)の約数でNN以下であるものの総数が問題になる。これがo(N2){o}(N^2)であればもっとよい評価ができることになるけれど、よくわからなかった。 ↩︎

  2. 仮に「箱」が区別できる設定だったら単にEGFの掛け算をすればいいので[xn/n!]n1DM(x)n[x^n/n!]\sum_{n \ge 1}\mathcal D_M(x)^nでよくて、箱の区別を無くすために各項をn!n!で割るとn1DM(x)n/n!=exp(DM(x))\sum_{n \ge 1}\mathcal D_M(x)^n/n! = \exp(\mathcal D_M(x))になる。 ↩︎

2021年11月9日火曜日

ACL Beginner Contest F - Heights and Pairs

hih_iが互いに異なる場合は(2N1)!!=(2N)!/2NN!(2N-1)!! = (2N)! / 2^N N!通りという有名事実がある(後述)。以下、pi=(2i1)!!p_i = (2i-1)!!とする。

身長hhの人がaha_h人いるとする。ii人の同じ身長を持つ人からjjペアを作る場合の数は(i2j)pj\binom{i}{2j}p_jである。F(x)=hj(ah2j)pjxjF(x) = \prod_h \sum_j \binom{a_h}{2j}p_j x^{j}とすると、F(x)F(x)xkx^kの係数[xk]F(x)[x^k]F(x)2N2N人から身長が同じであるようなkkペアを作る場合の数に等しい。このとき、残った人、つまり2(Nk)2(N-k)人でペアを作る場合の数はpNkp_{N-k}通りであり、あとは包除でk(1)kpNk([xk]F(x))\sum_k(-1)^kp_{N-k}([x^k]F(x))通りと求まる。

(2N1)!!(2N-1)!!について

2N2N人に1,2,...,2N1, 2, ..., 2Nという番号を振ってこの順に並べ、二人ずつ取り除いていくとする。この際、ペアのうちの一人目は必ずその時点で番号が最小の人を選ぶ、としておくと完全な数え上げができる。一人目には2N12N-1通りの人がマッチできて、二人目には2N32N-3通りの人がマッチできて、……となるので二重階乗の記号を使うと(2N1)!!(2N-1)!!と表せる。

別の考え方もできる。2N2N人を自由な順番で一列に並べて、隣り合う二人をくっつけてNNペアを作ることにする。並べ方は(2N)!(2N)!通りあるが、二人組の内の順番は考慮しないので2N2^Nで割る必要がある。さらに、NNペアの順番も考慮しないのでN!N!で割って(2N)!/(2NN!)(2N)!/(2^N N!)と求まる。これらは二重階乗の公式も与えている:

(2N1)!!=(2N)!2NN! (2N-1)!! = \frac{(2N)!}{2^N N!}

2021年11月7日日曜日

エイシング プログラミングコンテスト 2020 F - Two Snuke

2(s1+n1+u1+k1+e1)+Δs+Δn+Δu+Δk+ΔeN2(s_1+n_1+u_1+k_1+e_1)+\Delta s + \Delta n + \Delta u + \Delta k + \Delta e \le Nで各変数は非負という条件下でΔsΔnΔuΔkΔe\Delta s\Delta n\Delta u\Delta k\Delta eの総和を取る問題。

d=N(Δs+Δn+Δu+Δk+Δe)d = N - (\Delta s + \Delta n + \Delta u + \Delta k + \Delta e )とするとき、ΔsΔnΔuΔkΔe\Delta s\Delta n\Delta u\Delta k\Delta es1+n1+u1+k1+e1d/2s_1+n_1+u_1+k_1+e_1 \le d/2の非負整数解の数だけ足される。これは(d/2)+55)\binom{\lfloor d/2 \rfloor)+5}{5}通りだから、求める数は

d1+d2+d3+d4+d5+d=N(d/2+55)d1d2d3d4d5 \sum_{d_1+d_2+d_3+d_4+d_5 + d = N}\binom{\lfloor d/2 \rfloor+5}{5}d_1d_2d_3d_4d_5

に等しい。FPSで考えるなら、nnxn=x/(1x)2\sum_n nx^n = x/(1-x)^2を5乗してi(i/2+55)xi\sum_{i} \binom{\lfloor i/2 \rfloor +5}{5}x^{i}という謎の級数を掛けると、xNx^Nの係数が求める数である。後者にBerlekamp-Massey法を適用してみると21次の多項式を分母とする有理式が出てくるので、xNx^Nの係数はBostan-Mori法で求まる。


そもそも一般にi/2\lfloor i/2 \rfloorの多項式を係数にするFPSが有理べき級数になるのかどうかとか、基本的なことがよくわかっていなかった。

与えられたFPSf(x)f(x)の偶数項、奇数項がそれぞれ有理べき級数g(x),h(x)g(x), h(x)の偶数項、奇数項と一致しているなら、ggの(0-basedで)偶数項だけ抜き出したもの(g(x)+g(x))/2(g(x)+g(-x))/2hhの奇数項だけ抜き出したもの(h(x)h(x))/2(h(x)-h(-x))/2を足せばf(x)f(x)になるので、確かに有理べき級数になる。一般にi/k\lfloor i/k \rfloorの多項式を係数とするFPSについても、11kk乗根を考えれば同じことが言える。

一般化すると、次の性質を満たすf(x)=nanxnf(x) = \sum_n a_n x^nは有理べき級数と言える:

MM個の有理べき級数fk(x)=nbk,nxn (k=0,1,...,M1)f_k(x) = \sum_n b_{k, n}x^n \ (k=0, 1, ..., M-1)NNN \in \mathbb Nが存在して、nNn \ge Nならばan=bnmodM,na_n = b_{n \bmod M, n}