週刊 spaghetti_source

Dinkelbach algorithm

09:36

関数 f(x) と g(x) の比 f(x)/g(x) の最小化問題のことを分数計画（fractional programming），比の最小化などとよびます．発散を防ぐため g(x) > 0 を仮定します．今回紹介する Dinkelbach algorithm は，連続・離散どちらについても，広い範囲で適用できる比の最小化に対するアルゴリズムです．

比の最小化に対する基本的な方針は，[Jagannathan 1966] と [Dinkelbach 1967] が確立したパラメトリックアプローチです．これは t* = min f(x)/g(x) と思って分母を払う気持ちで h(t) = min { f(x) - t g(x) | x \in X } というパラメタ t をもつ関数を考えると

• h(t) > 0 iff t は最適値より大
• h(t) = 0 iff t は最適値
• h(t) < 0 iff t は最適値より小

という関係がなりたつ，というものです．この関係から，比の最小化は h(t) の零点計算に帰着するのですが，h(t) は狭義単調減少する凹関数なので（∵ h は「x でパラメタ付けられた傾きが負の線型関数」の min），零点計算がわりと簡単に行えます．h(t) の関数値評価では x に関する最小化問題を解く必要があり，一般にはもとの最小化とどちらが簡単かは言えませんが，よく扱われる問題では h(t) の関数値評価のほうが簡単です（eg., f, g が両方共線型なら f - t g も線型，f 凸, g 凹, t ≧ 0 なら f - t g 凸など）．

さて，h(t) = 0 を満たす t の計算法で主流のものが 3 つあります．

• 二分探索 [see: Radzik 1992]
• Dinkelbach algrithm [Dinkelbach 1967]
• Parametric search [Meggido 1979]

二分探索は単純に h(t) = 0 を二分探索するもので普通ですが，特に早くはありません．Parametric Search は二分探索を t を不定元だと考えて実行して t の分岐木を作るタイプのアルゴリズムで，理論計算量では最良のものができることが多いのですが，実装が複雑化して定数が大きい [Cole 1987] と言われています．最近は少し実用的になったみたい [van Oostrum-Veltkamp 2004] ですが，自分の印象では，あまり実装したくないタイプです（※ この論文はタイトルが「Parametric search made practical」．つまり，実用的になったよ！で論文が書ける程度には複雑）．

Dinkelbach algorithm は，実装容易かつ実用上早い手法です．手続きを述べると以下のようになります．

初期値 x[0] を与える
for k = 1, 2, ...
  t[k] := f(x[k]) / g(x[k])
  x[k+1] := argmin { f(x) - t[k] g(x) | x \in X }

収束性として次の結果が知られています．

• コンパクト集合上の連続関数に対して，このアルゴリズムは収束
• 線型などの多くのケースで超一次収束（cf. 二分探索は一次収束）

• 有限集合上の離散関数に対し，このアルゴリズムは収束
• {0,1}線型や劣モ/優モなどの多くの例で反復回数が強多項式（cf. 二分探索は最大関数値がかかる）

歴史

連続：特殊な場合については，Dinkelbach よりも前から同種の手続きが知られていました．例えば線型の場合は [Saaty-Gass 1954] や [Isbell-Marlow 1956]，二次の場合は [Ritter 1962] があります．[Jagannathan 1966] と [Dinkelbach 1967] は一般の連続関数について議論を行い，凸/凹 の場合について Dinkelbach が上のアルゴリズムを定式化しました．一般の収束証明は [Jeflea 2003] にあります（簡単なので，もっと前から知られていそうですが）．

離散：[Radzik 1992] が {0,1} 整数線型計画に対してこのアルゴリズムの計算量が強多項式で抑えられることを証明しました（それまでは parametric search が主流）．その後は具体的な問題に対するケーススタディが多いように見受けられます．また，ごく最近劣モジュラ / 優モジュラに対する同種の結果 [Kawahara-Nagano-Okamoto 2011] が出ていて，個人的には興味深いです．

他にも周辺はいろいろありますが，ありすぎて，あんまり簡潔にまとまっていないのが現状だと思います．「Recent Development と Fractional Programming」で検索して出てくるサーベイや，Handbook of Global Optimization を見るのがいいんでしょうか．

なお，こういった経緯があるのでアルゴリズム名も様々な名前で呼ばれています．自分の印象では連続系の人は Dinkelbach algorithm と呼ぶことが多く，離散系の人は Radzik の呼び名にしたがって Newton method や Discrete Newton method と呼ぶことが多いように思います．

実装例

ここでは離散分数計画における典型例である「最小比全域木問題」を実装します．これは無向グラフ G に2種類の非負重み a, b が与えられたとき w[a](T)/w[b](T) を最小化する全域木を求める問題です．

解法は上に述べたとおりで，パラメタ化された最小化問題 h(t) = min { w[a](T) - t w[b](T) | T \in Spanning Trees } を繰り返しとけば良いのですが，この問題は {0,1} 整数線型計画，すなわち w[a] - t w[b] = w[a - t b] という性質があるので，h(t) の評価は普通の最小全域木のルーチンで重みを a - t b に設定したものが使えます．

C++11 の勉強がてら，それっぽい実装ですが，まあ読めるとおもいます．コンパイルは g++ 4.7 で確認．

struct Graph {
  int n;
  struct Edge { int u, v, a, b; };
  vector<Edge> E;
  Graph(int n) : n(n) { }
  void addEdge(int u, int v, int a, int b) { E.push_back({u, v, a, b}); }

  double fcost, gcost;
  double minimumSpanningTree(double t) {
    fcost = gcost = 0;
    UnionFind unionFind(n);
    sort(E.begin(), E.end(), [t](Edge e, Edge f) { 
      return e.a - t*e.b < f.a - t*f.b; });
    for (Edge e : E) 
      if (unionFind.unite(e.u, e.v)) fcost += e.a, gcost += e.b;
    return fcost - t * gcost;
  }
  double fractionalMin() {
    for (double t = 0; ; t = fcost / gcost)
      if (fabs( minimumSpanningTree(t) ) < 1e-8) return t;
  }
};

参考文献

• T. Saaty and S. Gass (1954): "Parametric objective functions (Part 1)", Operations Research, vol. 2, pp. 316-319.
• J. R. Isbell and W. H. Marlow (1956): "Attrition games", Naval Research Logistic Quarterly, vol. 3, pp. 71-93.
• K. Ritter (1962): "Ein Verfahren zur Losung parameterabhangiger, nichtlinearer Maximum-Probleme", Unternehmensforschung, vol. 6, pp. 149-166.
• R. Jagannathan (1966): "On some propeties of progamming problems in parametric form pertaining to fractional programming", Management Science, vol. 12, pp. 609-615.
• W. Dinkelbach (1967): "On nonlinear fractional programming", Management Science, vol. 13, pp. 492-498.
• R. Cole (1987): "Slowing down sorting networks to obtain faster sorting algorithms", Journal of the ACM, vol. 37, no. 1, pp. 200-208.
• T. Radzik (1992): "Newton's Method for Fractional Combinatorial Optimization", Proceedings, 33rd Annual Symposium on Proceedings of Foundation of Computer Science, pp. 659-669.
• A. Jeflea (2003): "A parametric study for solving nonlinear fractional problems", Analele Stiintifice ale Universitatii Ovidus Constanta Seria: Mathematica, vol. 11, no. 2, pp. 72-92.
• R. van Oostrum and R. C. Veltkamp (2004): "Parametric search made practical", Computational Geometry: Theory and Applications - Special issue on the 18th annual symposium on computational geometry, vol. 28, no. 2-3, pp. 75-88.
• Y. Kawahara, K. Nagano, and Y. Okamoto (2011): "Submodular fractional programming for balanced clustering", Pattern Recognition Letters vol. 32, no. 2, pp. 235-243.