Implementation of various algorithms
#include "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;
VAL は符号付き整数型とするtemplate<class FUNC> pair<VAL, VAL> solve(const FUNC &lag, VAL K)
lag(lambda): $x^{\mathrm{opt}} = \mathrm{argmin}_{x} (f(x) + \lambda g(x))$ としたときのペア値 $(f(x^{\mathrm{opt}}), g(x^{\mathrm{opt}}))$ を返す関数
template<class FUNC> pair<VAL, VAL> solve(const FUNC &lag, VAL L, VAL R)
インターフェースについては同上。
Alien Trick は、とくに Monge 単一始点 d-辺最短路問題への応用 (Alien DP) がよく問われる。以下のページに、Alien DP における使用例を載せている。
#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);
}
};