Hatena::Grouptopcoder

skyaozoraの日記

2018-12-20

ConvexHullTrickを分割統治法で殴る

09:31

Competitive Programming (1) Advent Calendar 2018の20日目の記事です。この記事では、最近流行り(?)のConvexHullTrickを少し別の角度から見て行こうと思います。

では、まずこの問題を見て見ましょう

CF189DIV1C Kalila and Dimna in the Logging Industry

問題概要をきわめてざっくり説明すると、長さn(<=10^5)の2つの数列a,bが与えられて、dp[i]=min(j<i)dp[j]+a[i]*b[j]というDPを解けばいいです。ただし、aは狭義単調増加でbは狭義単調減少です。</ppp>

ConvexHullTrickによる解法

この問題は、2次元平面上に直線がどんどん追加されていき、あるx座標においてy座標が最も小さくなる直線を求めるということを繰り返すという問題ということができます。そして、「追加クエリにおける直線の傾きが単調」「最小値クエリにおけるxが単調」の2つを満たしているので、ConvexHullTrick(最近はCHTと略されることが多いらしいですね・・・最初に見たときは何のことか全く分かりませんでした)を用いてO(n)で解くことができます。ConvexHullTrickの詳しいことは例えばこちらの記事をご覧ください。ConvexHullTrickを用いた解答コードは以下のようになります。

#include<bits/stdc++.h>
using namespace std;
typedef long long lint;
#define REP(i,a,b) for(int i=a;i<b;i++)
#define rep(i,n) REP(i,0,n)
lint a[100010],b[100010],dp[100010];
int deq[100010];
lint cal(int i,int j){return dp[j]+b[j]*a[i];}
double cro(int i,int j){return (0.0+dp[j]-dp[i])/(b[i]-b[j]);}
int main()
{
	int s=0,t=1,n;
	cin>>n;
	rep(i,n) scanf("%I64d",&a[i]);
	rep(i,n) scanf("%I64d",&b[i]);
	dp[0]=0;deq[0]=0;
	REP(i,1,n){
		while(s+1<t && cal(i,deq[s])>cal(i,deq[s+1])) s++;
		dp[i]=cal(i,deq[s]);
		while(s+1<t && cro(deq[t-1],i)<=cro(deq[t-1],deq[t-2])) t--;
		deq[t++]=i;
	}
	cout<<dp[n-1]<<endl;
}

これ、コードの行数は少なくて極めてシンプルなコードに見えますよね?でも自分はこのコードを書く時に意外と苦労した記憶があります。配列を使ってdequeの中身を操作するのですが、どのタイミングでiteratorの値を動かせばいいのか、中身が十分に入ってなくて値の比較ができないというのはどこで判定すればいいのか・・・などなど、コードにも-1や-2をどこで入れるのかなど細かいことに気を使う印象です。

分割統治によって最小値たちを効率よく求める方法

さて、今の問題はいったん置いといて次の問題について考えてみましょう

n×nの2次元配列cが与えられる。各行ごとに最小値をとる列を全て求めよ。ただし、各行ごとに最小値をとる列の場所は広義単調増加になっている

つまり各行ごとに最小値をとる列を赤く塗るとこのようになっているということですね。

f:id:skyaozora:20181220083627p:image

これはもちろん愚直にやればO(n^2)ですが、この性質を利用してより効率的に求めることができます。

まず、真ん中の行の最小値をとる列の場所を調べます。そうすると単調性より、その場所の(図の上での)右上と左下には最小値をとる場所は存在しないことが分かります。

f:id:skyaozora:20181220083624p:image

次に、左上と右下のブロックそれぞれにおいて同様に真ん中の行の最小値をとる場所を調べて・・・ということを再帰的に繰り返すと、全体でO(nlogn)で求めることができます。より詳しい解説はこちらの記事をご覧下さい。

これを使うと、例えば2次元DPで、配列dp[i]の各値が配列dp[i-1]の各値から遷移され、かつその遷移元の位置が遷移先の位置に対して広義単調であることが分かる場合(図で表すと下のようになります。これを私は勝手に「遷移が交差しない」状態と呼んでいます。)、dp[i-1]からdp[i]への遷移は愚直にやればO(n^2)のところを上に述べたテクニックを使うことでO(nlogn)で全て行うことができます。

f:id:skyaozora:20181220083629p:image

ConvexHullTrickの代わりに分割統治を使う

で、この分割統治法がConvexHullTrickとどういう関係があるのかというと、「追加クエリにおける直線の傾きが単調」「最小値クエリにおけるxが単調」の2つを満たしている場合、クエリに対して最小値をとる直線が追加された順番は単調になるという性質が成り立つからです。実際にdequeに直線が追加されたり捨てられたりする順番を考えれば成り立つのが分かると思います。また今回の問題に関しても、dp[j]+a[i]*b[j]が最小となるjがiに対して広義単調増加になる(あるi=i1に対して最小となるjをj1とした時、i=i2>i1に対して最小となるjがj1より小さいはずがない)のは、aが狭義単調増加でbが狭義単調減少なことから分かります。よって、前の節で述べた「遷移が交差しない」という条件を満たします。

しかし、今回の問題では遷移前の状態というのがまさに遷移によって求まるものなので、前節のような方法で求めることはできません。では駄目かというと、ここでもまた分割統治法を使うことで解決できます。

さらなる分割統治法で最適解を求める

結論から言うと、1次元DPで「遷移が交差しない」という条件が満たされている場合、以下のアルゴリズムを用いれば範囲[0,n)の全ての値の最適解を求めることができます。

1.範囲[0,m)の最適解を求める

2.範囲[m,n)の範囲[0,m)から遷移される部分に関する最適解を求める

3.範囲[m,n)の最適解を求める

ここでm=n/2とします。手順1,3に関しては再帰的に行い、手順2に関しては上に述べた分割統治法を行います。このアルゴリズムの詳しい解説に関してはこちらの記事をご覧下さい。手順2の計算量がO(nlogn)である場合、全体でO(nlog^2n)で全ての値の最適解を求めることができます。このアルゴリズムを用いた最初の問題の解答コードは以下のようになります。

#include<bits/stdc++.h>
using namespace std;
typedef long long lint;
#define REP(i,a,b) for(int i=a;i<b;i++)
#define rep(i,n) REP(i,0,n)
lint inf=1e18;
lint a[100010],b[100010],dp[100010];
lint cal(int i,int j){return dp[j]+b[j]*a[i];}
void cal2(int l1,int h1,int l2,int h2){
	if(l1>=h1) return;
	int mi=(l1+h1)/2,it=l2;
	REP(i,l2,h2){
		if(cal(mi,i)<cal(mi,it)) it=i;
	}
	dp[mi]=min(dp[mi],cal(mi,it));
	cal2(l1,mi,l2,it+1);
	cal2(mi+1,h1,it,h2);
}
void rec(int lo,int hi){
	if(lo+2>hi) return;
	int mi=(lo+hi)/2;
	rec(lo,mi);
	cal2(mi,hi,lo,mi);
	rec(mi,hi);
}
int main()
{
	int s=0,t=1,n;
	cin>>n;
	rep(i,n) scanf("%I64d",&a[i]);
	rep(i,n) scanf("%I64d",&b[i]);
	dp[0]=0;
	REP(i,1,n+10) dp[i]=inf;
	rec(0,n);
	cout<<dp[n-1]<<endl;
}

解答の長さ自体は最初のより長いですが、本当に値を渡して、その範囲の探索を行うだけなので少なくとも個人的にはこちらのコードの方が頭を使わずに書けました。また計算量もO(n)からO(nlog^2n)へと増えているのですが、感覚としてこのアルゴリズムは定数倍が極めて早く、O(nlogn)と同等か下手したらそれより早いイメージがあります。実際、実行時間も最初のコードが78msなのに対しこのコードは108msで済んでいます。

もちろん、この方法はConvexHullTrickが使える場合に限らず、遷移元と遷移先の関係に関して単調性が満たされてさえいれば使えるのですが、そのような単調性を満たす自然な設定にしようとすると大体ConvexHullTrickも使える形になってる場合が多いと思います。なので、単調性を利用した解法を思いついたがそれだとTLEするといった場合は、式変形などを経る事でConvexHullTrickが使える形に帰着できないか考えてみるといいと思います。

最後に

結局話としては想定解よりも計算量の大きい解法を説明しただけで、この話は知らなくても解ける問題というのが減るわけではないと思います。でも実装のアプローチを複数持っておくこと自体は悪くないと思いますし、問題の背景にこういう性質が隠れてるのを知るのもそれはそれで面白いんじゃないかと思います。この記事が少しでも皆さんの競プロ人生の役に立てば幸いです。

参考文献

Dynamic Programming Optimizations

Convex-Hull Trick

Totally Monotone Matrix Searching (SMAWK algorithm)

動的計画法高速化テクニック(オフライン・オンライン変換)

2018-07-04

AOJ2563 The J-th Number

22:28

問題概要

英語だけど短いし問題文を読んで

続きを読む

2018-07-03

AOJ2290 Attack the Moles

16:57

問題概要

問題文

日本語だし書くの面倒だから問題文を見て

続きを読む

2018-07-02

AOJ2560 Point Distance

23:10

問題概要

問題文

2次元平面の各格子点上に9個以下の点がある。(本質だけに絞ると)距離が√0、√1、√2・・・になる点のペアがそれぞれ何個ずつあるかを答えよ。座標のサイズは縦横1024以下。

続きを読む

2017-12-12

「インラインDP」というテクニックに関して

22:17

12月14日,式に一部間違いがあったので修正しました

競技プログラミングAdventCalender2017の12日目の記事です。この記事では、最近コンテストでしばしば見かける「インラインDP」というテクニックに付いて説明します。まぁ「インラインDP」と呼んでいるのはおそらく自分だけで、他の人たちは「segtreeを使ったDP高速化のやつ」とか「実家」(典型テクニックだから、という意味で)とか呼んでいますが、おそらくこの呼び方が自分では一番中身をあらわしていると思うので、それで行こうと思います。

さて、インラインDPの例題として次の問題の解法を考えていきましょう。

ARC073F問題 Many Moves

まずこの問題を解くためには、

dp[何個目までのクエリを処理したか][1つ目のコマの場所][2つ目のコマの場所]:=それを達成する最小の秒数

という状態を持ってDPすれば解けます。しかし、i個目のクエリを処理した直後は、必ず片方のコマはx_iに置かれているので、

dp[何個目までのクエリを処理したか][もう片方のコマの場所]:=それを達成する最小の秒数

という状態を持てばいいです。これは以下のような2次元DPで、i個目までのクエリを処理した、というもののiを0からNまで加算しつつ埋めていけば解けます。しかしこれでは状態量が最大NQとなり、実行時間的にもメモリ的にも到底収まりません。

f:id:skyaozora:20171212085509p:image

そこで、計算量を改善するために、まずは遷移を考えてみましょう。dp[i][j]、つまりi番目までのクエリを処理して2つのコマがそれぞれ位置x_iと位置jにある場合、次のi+1番目のクエリを処理するために、

・位置x_iにあるコマを位置x_{i+1}に動かす

・位置jにあるコマを位置x_{i+1}に動かす

という2通りの方法があり、それぞれに対応するDPの更新式は

dp[i+1\][j\]=min(dp[i+1\][j\]\:,\:dp[i\][j\]+|x_i-x_{i+1}|)

dp[i+1\][x_i\]=min(dp[i+1\][x_i\]\:,\:dp[i\][j\]+|j-x_{i+1}|)

となります。

ここで、更新式の遷移元ではなく遷移先に注目すると(界隈でよく言われている呼び方で言うと、配るDPから貰うDPに変換すると)

dp[i+1\][j\]=dp[i\][j\]+|x_i-x_{i+1}| (j!=x_{i})

dp[i+1\][j\]=\min_{k=1,...,N}dp[i\][k\]+|k-x_{i+1}| (j=x_{i})

となります。すると、dp[i+1]という1次元配列(つまりDPテーブルのi+1行目)は、dp[i]という1次元配列(つまりDPテーブルのi行目)と比較して、(全体に|x_i-x_{i+1}|という定数を足した上で)x_{j+1}番目の要素しか変化していないということがわかります。このように、DPテーブルのi行目とi+1行目が1箇所しか異なっていない場合、インラインDPというテクニックを用いて計算量を減らすことが出来ます。

どうするかというと、まずDPテーブルを愚直に2次元で持つのではなく、1次元で持ちます。そして最初はその1次元配列をDPテーブルの1行目で初期化し、次にその配列の1箇所のみを更新してDPテーブルの2行目を表し、その配列の1箇所をのみ更新してDPテーブルの3行目を表し・・・ということを繰り返していくことで、更新回数がO(Q)回、空間計算量もO(N)でDPテーブルをQ行目まで計算することが出来ます。これがインラインDP(私が勝手に名付けただけですが)というテクニックです。

f:id:skyaozora:20171212085506p:image

ただし、その1回の更新では

\min_{k=1,...,N}dp[k\]+|k-x_{i+1}|

というものを求めなければいけないので、これを愚直に計算すると1回の更新につきO(N)かかってしまい、全体の計算量がO(NQ)になってしまいます。しかし、これも競プロの問題でよく使われるテクニックで解決できます。これにはまず各k=1,...,Nに対して

dpl[k\]=dp[k\]+N-k

dpr[k\]=dp[k\]+k

を満たす配列を持っておきます。すると、

\min_{k=1,...,N}dp[k\]+|k-x_{i+1}|

=min(\min_{k=1,...,x_{i+1}-1}dpl[k\]-N+x_{i+1}\:,\:\min_{k=x_{i+1},...,N}dpr[k\]-x_{i+1})

となるので、配列dpl,dprそれぞれのある連続区間の最小値を求めればよくなります。連続区間の最小値は皆さんご存知のようにsegtreeを用いればO(logN)で求めることが出来るので、1次元配列をそのまま持つのではなく、

dpl[k\]=dp[k\]+N-k

dpr[k\]=dp[k\]+k

という補正を行ったdpl,dprをそれぞれsegtreeとして持つことで、全体の計算量O(QlogN)でこの問題を解くことが出来ました。このように、この「インラインDP」というテクニックを用いる場合、単純な1次元配列を用いるだけでは解けずsegtreeやBITといったデータ構造を用いる必要があることが多いです。そのため「segtreeを用いたDPの高速化テク」と呼んでる人たちが多いですが、それだともうちょっと指す範囲が広そうなのと、自分はこのテクニックで一番重要なのは

「1次元配列の一部だけを更新していくことによって、実際には2次元であるDPを上から下まで計算できる」

という点だと思っているので、そこを表せそうな「インラインDP」という名前をつけてみました。

最後にこの問題の解答コードを載せます

#include<bits/stdc++.h>
using namespace std;
#define pb push_back
#define pf push_front
typedef long long lint;
typedef complex<double> P;
#define mp make_pair
#define fi first
#define se second
typedef pair<int,int> pint;
#define All(s) s.begin(),s.end()
#define rAll(s) s.rbegin(),s.rend()
#define REP(i,a,b) for(int i=a;i<b;i++)
#define rep(i,n) REP(i,0,n)
#define N 262144
lint dat[2][N*2+10];
lint inf=1145141919810364364LL;
//[a,b)の最小値
//外からは(a,b,0,0,n)として呼ぶ
lint query(int a,int b,int id,int k=0,int l=0,int r=N){
	if(r<=a || b<=l) return inf;
	if(a<=l && r<=b) return dat[id][k];
	lint vl=query(a,b,id,k*2+1,l,(l+r)/2);
	lint vr=query(a,b,id,k*2+2,(l+r)/2,r);
	return min(vl,vr);
}
void update(int id,int k,lint a){
	k+=N-1;
	dat[id][k]=a;
	while(k>0){
		k=(k-1)/2;
		dat[id][k]=min(dat[id][k*2+1],dat[id][k*2+2]);
	}
	return;
}
int x[252521];
int main()
{
	lint sum=0,out=inf,out2=inf;int n,q,a,b;
	cin>>n>>q>>a>>b;x[0]=a;
	rep(i,2) rep(j,N*2+5) dat[i][j]=inf;
	update(0,b,n-b);update(1,b,b);
	rep(i,q) cin>>x[i+1];
	rep(i,q){
		lint dif=abs(x[i]-x[i+1]);
		sum+=dif;
		lint ne=min(query(0,x[i+1]+1,0)-(n-x[i+1]),query(x[i+1],n+1,1)-x[i+1]);
		update(0,x[i],ne-dif+n-x[i]);update(1,x[i],ne-dif+x[i]);
	}
	rep(i,n+1) out=min(out,dat[0][i+N-1]-(n-i));
	rep(i,n+1) out2=min(out2,dat[1][i+N-1]-i);
	assert(out==out2);
	cout<<out+sum<<endl;
}

さて、インラインDPの説明をするためにあえて無視をしていたのですが、この問題を解くためには「全体に|x_i-x_{i+1}|という定数を足した上で」1箇所を更新しなければいけませんでした。StarrySky木などを使えばSegtreeの全体に加算するということも出来ますが、このコードでは全体にいくつ足されたか、という情報をまた別に持っています(ただし、dp[x_{i+1}\]の最小値を更新する際にその分を考慮しなくてはいけません)これもよく使われるテクニックなので抑えておきましょう。

関連問題

さて、この問題以外にもこの「インラインDP」のテクニックを使って解くことの出来る問題を載せておきます(ここに乗せるということ自体が一種のネタバレで申し訳ございませんがご了承ください

ARC085F問題 NRE

ARC056D問題 サケノミ

AGC015E問題 Mr.Aoki Incubator

最後に

さて、この記事は参考になりましたでしょうか。式や図が分かりにくくて理解できないって場合はこの記事のコメント欄やtwitterまでお願いします。出来る限り対処したいと思います。この記事を読んでくれた皆さんが来年以降も競技プログラミングを楽しんでくれることを、そしてこの記事が少しでもお役に立てることを祈りつつ、良いお年を!