algorithm

Implementation of various algorithms

View the Project on GitHub drken1215/algorithm

:warning: Alien Trick --- min f(x) s.t. g(x) = K (by ラグランジュ緩和)
(Optimization/alien_trick.hpp)

概要

問題設定

変数 $x$ について、なんらかの「個数」を表す関数を $g(x)$ とおく(たとえば $x$ がグラフのパスを表し、$g(x)$ がパス中の辺の本数を表すなど)。 $g(x)$ は $L$ 以上 $R$ 以下の整数値をすべてとる関数とする。

$K$ を $L \le K \le R$ を満たす正の整数として、次の最適化問題を考える ($f(x)$ は整数値をとる関数とする)。

$\min_{x} f(x)$ $\mathrm{s.t.}$ $g(x) = K$

$g(x) = k$ となる $x$ についての $f(x)$ の最小値を $h(k)$ とする。 $h(k)$ の定義域は $k = L, L+1, \dots, R$ であり、求めたい値は $h(K)$ である。 ここで、次の仮定をおく。


【仮定】

$h(k)$ が定義域内で離散凸性を満たすとする。すなわち、

$h(k-1) + h(k+1) \ge 2h(k)$

を満たすとする。


この仮定を満たす具体的な状況については、Alien Trick で、

という 4 つの場面が紹介されている。

ラグランジュ緩和

上記の仮定の下で、ラグランジュ強双対性が成り立つことが知られている。

を与えている。 具体的には、ラグランジュ双対変数 $\lambda$ を整数として、ラグランジュ緩和して得られる最適化問題

$\max_{\lambda} \min_{x} (f(x) + \lambda(g(x) - K))$

の最適値は、もとの最適化問題の最適値に一致する。補足として、ラグランジュ強双対性が成り立たなくても、以下のことが成り立つ。

ここでは、三分探索法よりも高速な二分探索法を考える。

ラグランジュ緩和問題の解き方

$\lambda$ を固定して、最適化問題

$P(\lambda): \min_{x} (f(x) + \lambda g(x))$

を考える(最適値を $L(\lambda)$ とおく)。 この最適化問題 $P(\lambda)$ は、$\min_{x} (f(x) + \lambda(g(x) - K))$ と等価であることに注意する。 また、最適化問題 $P(\lambda)$ には、Fenchel 双共役性に基づく幾何学的解釈も存在する(Alien Trick を参照)。

さて、$g(x) = K$ となるように $\lambda$ をコントロールすることを考える。 最適化問題 $P(\lambda)$ の最適解 $x^{\mathrm{opt}}$ を 1 つとり、$k(\lambda) = g(x^{\mathrm{opt}})$ とおくと、 $k(\lambda)$ は $\lambda$ について広義単調減少となる(Alien Trick を参照)。

よって、二分探索法によって、

$k(\lambda_{1}) \le K \le k(\lambda_{2})$ ($\lambda_{1} = \lambda_{2} + 1$)

を満たすような $\lambda_{1}, \lambda_{2}$ を求めることができる。 このとき、$\lambda = \lambda_{1}, \lambda_{2}$ のいずれかについて、 $g(x) = K$ であるような最適化問題 $P(\lambda)$ の最適解 $x$ が存在することが示せる。 すなわち、

$h(K) = \max(L(\lambda_{1}) - \lambda_{1} K, L(\lambda_{2}) - \lambda_{2} K)$

が成立する。

復元について

$k(\lambda_{1}) \le K \le k(\lambda_{2})$ ($\lambda_{1} = \lambda_{2} + 1$)

を満たすような $\lambda_{1}, \lambda_{2}$ を求めることができても、実際に

$\min_{x} f(x)$ $\mathrm{s.t.}$ $g(x) = K$

の最適解 $x$ を復元することは、一般には難しいと考えられる。ただし、以下に挙げるような、最適解に交換性があるような個別の場合については、最適解を復元できる場合がある。

これらの場合には、$g(x)$ が $K$ よりも大きい最適解と小さい最適解とから、最適性を保ちながら $g(x) = K$ であるような解を構成できる。

仕様

template<class VAL> struct AlienTrick;

制約条件が g(x) = K の場合

template<class FUNC> pair<VAL, VAL> solve(const FUNC &lag, VAL K)

制約条件が L <= g(x) <= R の場合

template<class FUNC> pair<VAL, VAL> solve(const FUNC &lag, VAL L, VAL R) 

インターフェースについては同上。

使用例

Alien Trick は、とくに Monge 単一始点 d-辺最短路問題への応用 (Alien DP) がよく問われる。以下のページに、Alien DP における使用例を載せている。

問題例

参考文献

Required by

Verified with

Code

#pragma once
#include <bits/stdc++.h>
#include <cassert>
using namespace std;

// Alien Trick (by ラグランジュ緩和)
// lag(λ): x✳︎ = argmin_{x} {f(x) + λ g(x)} としたときのペア値 {f(x✳︎), g(x✳︎)} を返す関数 (最大化の方では argmax)
// 返り値: {最適値, そのときの λ}
template<class VAL> struct AlienTrick {
    // solver
    template<class FUNC> pair<VAL, VAL> solve(const FUNC &lag, VAL K) {
        // zero check
        auto [f0, g0] = lag(VAL(0));
        if (g0 == K) return {f0, VAL(0)};

        // find an initial value of binary search (by doubling)
        VAL lo, hi, flo, fhi, glo, ghi;
        if (g0 > K) {
            lo = 0, flo = f0, glo = g0;
            tie(hi, fhi, ghi) = find_upper_lambda(lag, K);
        } else {
            hi = 0, fhi = f0, ghi = g0;
            tie(lo, flo, glo) = find_lower_lambda(lag, K);
        }
        if (glo == K) return {flo, lo};
        if (ghi == K) return {fhi, hi};
        assert(glo > K && ghi < K);

        // binary search
        // invariant: lo < hi, g(lo) > K, g(hi) < K
        while (hi - lo > VAL(1)) {
            VAL mid = (lo + hi) / 2;
            auto [fmid, gmid] = lag(mid);
            if (gmid == K) return {fmid, mid};
            if (gmid > K) lo = mid, flo = fmid, glo = gmid;
            else hi = mid, fhi = fmid, ghi = gmid;
        }

        // When the minimizer skips g(x) = K due to ties, recover the optimum from the Lagrangian dual value
        VAL reslo = flo + lo * (glo - K), reshi = fhi + hi * (ghi - K);
        if (reslo >= reshi) return {reslo, lo};
        else return {reshi, hi};
    }
    template<class FUNC> pair<VAL, VAL> minimize(const FUNC &lag, VAL K) {
        return solve(lag, K);
    }
    template<class FUNC> pair<VAL, VAL> maximize(const FUNC &lag, VAL K) {
        auto reverse_lag = [&](VAL lambda) -> pair<VAL, VAL> {
            auto [f, g] = lag(-lambda);
            return {-f, g};
        };
        auto [f, g] = solve(reverse_lag, K);
        return {-f, -g};
    }

    // in case: L <= g(x) <= R
    template<class FUNC> pair<VAL, VAL> solve(const FUNC &lag, VAL L, VAL R) {
        assert(L <= R);
        auto [f0, g0] = lag(VAL(0));
        if (L <= g0 && g0 <= R) return {f0, VAL(0)};
        if (g0 > R) return solve(lag, R);
        else return solve(lag, L);
    }
    template<class FUNC> pair<VAL, VAL> minimize(const FUNC &lag, VAL L, VAL R) {
        return solve(lag, L, R);
    }
    template<class FUNC> pair<VAL, VAL> maximize(const FUNC &lag, VAL L, VAL R) {
        auto reverse_lag = [&](VAL lambda) -> pair<VAL, VAL> {
            auto [f, g] = lag(-lambda);
            return {-f, g};
        };
        auto [f, g] = solve(reverse_lag, L, R);
        return {-f, -g};
    }

    // find lo such that g(lo) >= R
    template<class FUNC> tuple<VAL, VAL, VAL> find_lower_lambda(const FUNC &lag, VAL K) {
        auto [f0, g0] = lag(VAL(0));
        VAL lo = 0, hi = 0, flo = f0, glo = g0;
        if (glo >= K) return tie(lo, flo, glo);
        lo = VAL(-1);
        tie(flo, glo) = lag(lo);
        while (glo < K) {
            assert(lo < hi);
            VAL dif = hi - lo;
            hi = lo, lo -= VAL(2) * dif;
            tie(flo, glo) = lag(lo);
        }
        return tie(lo, flo, glo);
    }

    // find hi such that g(hi) <= L
    template<class FUNC> tuple<VAL, VAL, VAL> find_upper_lambda(const FUNC &lag, VAL K) {
        auto [f0, g0] = lag(VAL(0));
        VAL lo = 0, hi = 0, fhi = f0, ghi = g0;
        if (ghi <= K) return tie(hi, fhi, ghi);
        hi = VAL(1);
        tie(fhi, ghi) = lag(hi);
        while (ghi > K) {
            assert(lo < hi);
            VAL dif = hi - lo;
            lo = hi, hi += VAL(2) * dif;
            tie(fhi, ghi) = lag(hi);
        }
        return tie(hi, fhi, ghi);
    }
};
#line 2 "Optimization/alien_trick.hpp"
#include <bits/stdc++.h>
#line 4 "Optimization/alien_trick.hpp"
using namespace std;

// Alien Trick (by ラグランジュ緩和)
// lag(λ): x✳︎ = argmin_{x} {f(x) + λ g(x)} としたときのペア値 {f(x✳︎), g(x✳︎)} を返す関数 (最大化の方では argmax)
// 返り値: {最適値, そのときの λ}
template<class VAL> struct AlienTrick {
    // solver
    template<class FUNC> pair<VAL, VAL> solve(const FUNC &lag, VAL K) {
        // zero check
        auto [f0, g0] = lag(VAL(0));
        if (g0 == K) return {f0, VAL(0)};

        // find an initial value of binary search (by doubling)
        VAL lo, hi, flo, fhi, glo, ghi;
        if (g0 > K) {
            lo = 0, flo = f0, glo = g0;
            tie(hi, fhi, ghi) = find_upper_lambda(lag, K);
        } else {
            hi = 0, fhi = f0, ghi = g0;
            tie(lo, flo, glo) = find_lower_lambda(lag, K);
        }
        if (glo == K) return {flo, lo};
        if (ghi == K) return {fhi, hi};
        assert(glo > K && ghi < K);

        // binary search
        // invariant: lo < hi, g(lo) > K, g(hi) < K
        while (hi - lo > VAL(1)) {
            VAL mid = (lo + hi) / 2;
            auto [fmid, gmid] = lag(mid);
            if (gmid == K) return {fmid, mid};
            if (gmid > K) lo = mid, flo = fmid, glo = gmid;
            else hi = mid, fhi = fmid, ghi = gmid;
        }

        // When the minimizer skips g(x) = K due to ties, recover the optimum from the Lagrangian dual value
        VAL reslo = flo + lo * (glo - K), reshi = fhi + hi * (ghi - K);
        if (reslo >= reshi) return {reslo, lo};
        else return {reshi, hi};
    }
    template<class FUNC> pair<VAL, VAL> minimize(const FUNC &lag, VAL K) {
        return solve(lag, K);
    }
    template<class FUNC> pair<VAL, VAL> maximize(const FUNC &lag, VAL K) {
        auto reverse_lag = [&](VAL lambda) -> pair<VAL, VAL> {
            auto [f, g] = lag(-lambda);
            return {-f, g};
        };
        auto [f, g] = solve(reverse_lag, K);
        return {-f, -g};
    }

    // in case: L <= g(x) <= R
    template<class FUNC> pair<VAL, VAL> solve(const FUNC &lag, VAL L, VAL R) {
        assert(L <= R);
        auto [f0, g0] = lag(VAL(0));
        if (L <= g0 && g0 <= R) return {f0, VAL(0)};
        if (g0 > R) return solve(lag, R);
        else return solve(lag, L);
    }
    template<class FUNC> pair<VAL, VAL> minimize(const FUNC &lag, VAL L, VAL R) {
        return solve(lag, L, R);
    }
    template<class FUNC> pair<VAL, VAL> maximize(const FUNC &lag, VAL L, VAL R) {
        auto reverse_lag = [&](VAL lambda) -> pair<VAL, VAL> {
            auto [f, g] = lag(-lambda);
            return {-f, g};
        };
        auto [f, g] = solve(reverse_lag, L, R);
        return {-f, -g};
    }

    // find lo such that g(lo) >= R
    template<class FUNC> tuple<VAL, VAL, VAL> find_lower_lambda(const FUNC &lag, VAL K) {
        auto [f0, g0] = lag(VAL(0));
        VAL lo = 0, hi = 0, flo = f0, glo = g0;
        if (glo >= K) return tie(lo, flo, glo);
        lo = VAL(-1);
        tie(flo, glo) = lag(lo);
        while (glo < K) {
            assert(lo < hi);
            VAL dif = hi - lo;
            hi = lo, lo -= VAL(2) * dif;
            tie(flo, glo) = lag(lo);
        }
        return tie(lo, flo, glo);
    }

    // find hi such that g(hi) <= L
    template<class FUNC> tuple<VAL, VAL, VAL> find_upper_lambda(const FUNC &lag, VAL K) {
        auto [f0, g0] = lag(VAL(0));
        VAL lo = 0, hi = 0, fhi = f0, ghi = g0;
        if (ghi <= K) return tie(hi, fhi, ghi);
        hi = VAL(1);
        tie(fhi, ghi) = lag(hi);
        while (ghi > K) {
            assert(lo < hi);
            VAL dif = hi - lo;
            lo = hi, hi += VAL(2) * dif;
            tie(fhi, ghi) = lag(hi);
        }
        return tie(hi, fhi, ghi);
    }
};
Back to top page