Trong phần này mình giới thiệu kĩ thuật bao lồi (convex hull trick) để cải tiến thuật toán quy hoạch động với thời gian $O(n^2)$ xuống còn $O(n\log n)$, thậm chí với một số trường hợp xuống còn $O(n)$, cho một lớp các bài toán. Trước hết chúng ta hãy xem xét bài toán quy hoạch động tổng quát sau:
$C[i] = \begin{cases} B[i], & \mbox{if } i = 1 \\ \min_{j < i}\{C[j] + A[i]B[j] \} , & \mbox{otherwise } \end{cases} \qquad (1)$
Tìm giá trị $C[n]$.
Từ công thức $(1)$, ta có thể dễ dàng thiết kế giải thuật $O(n^2)$ để giải bài toán như sau:
$C[1] \leftarrow B[1]$
for $i \leftarrow 2$ to $n$
tmp $\leftarrow +\infty$
for $j \leftarrow 1$ to $i-1$
tmp $\leftarrow \min\{$tmp, $C[j] + A[i]B[j]\}$
$C[i] \leftarrow$ tmp
return $C[n]$
Sau đây mình sẽ giới thiệu kĩ thuật bao lồi để giảm thời gian tính của thuật toán xuống còn $O(n \log n)$.
Kĩ thuật bao lồi
Trước hết chúng ta hãy tạm quên bài toán 1 để xét bài toán hình học sau:
$q_\ell = \min_{1 \leq i \leq n} a_i\cdot p_\ell + b_i$
Với mỗi điểm $p_\ell$, bằng cách duyệt qua tất cả các đường thẳng trong $D$, ta có thể tính được $q_\ell$ trong thời gian $O(n)$. Do đó tổng thời gian để tìm tất cả các điểm $q_1,q_2 \ldots, q_k$ là $O(kn)$. Tuy nhiên, ta có thể tìm các điểm này nhanh hơn:
Định lí 1 đúng cả với trường hợp tập các đường thẳng $D$ có đường thẳng với $a_i < 0$. Tuy nhiên, trong các ứng dụng với bài toán quy hoạch động, các $a_i$ luôn không âm. Ngoài ra, điều kiện $a_i \geq 0$ cũng làm cho việc phân tích thuật toán trở nên đơn giản hơn.
Trước hết chúng ta hãy xem ví dụ với $D = \{-x+y = 3, -x+2y = 4,$ $-x+6y = 16, -4x + 5y = 20\}$. Ta có thể thấy với mỗi điểm $(p,0) \in Ox$, giá trị nhỏ nhất tương ứng $(0,q)$ thỏa mãn $(p,q)$ thuộc phần gạch đỏ. Phần gạch đỏ đó gọi là bao lồi (do đó kĩ thuật này còn gọi là kĩ thuật bao lồi). Như vậy, dựa vào hình vẽ trên, đường thẳng $-4x + 5y = 20$ trong $D$ là dư thừa (tung độ của điểm có hoành độ $p$ trên đường thằng này luôn lớn hơn tung độ của điểm có cùng hoành độ trên một đường thẳng khác của $D$).
Dựa vào ví dụ trên, ta thấy nếu chúng ta có một cách biểu diễn bao lồi phù hợp, chúng ta có thể nhanh chóng xác định được tung độ $q$ nhỏ nhất ứng với hoành độ $p$. Ta có thể biểu diễn bao lồi bằng tập các interval và một đường thẳng tương ứng mỗi interval. Với ví dụ trên, tập các interval là $I_1 = [-\infty, -2],I_2 = [-2, 2], I_3 = [2, +\infty]$ và các đường thẳng tương ứng lần lượt là $-x+y=3, -x + 2y = 4, -x+6y = 16$. Định lý sau cho ta biết số interval trong biểu diễn bao lồi luôn nhỏ hơn hoặc bằng số đường thẳng.
Để đơn giản, ta giả sử trong $D$ không có hai đường thẳng nào song song. Nếu có hai đường thẳng song song, ta chỉ cần thay đổi một chút trong giả mã dưới đây mà không thay đổi thời gian tính của thuật toán. Thuật toán tìm bao lồi của $D$ như sau :
sort $\{d_1,d_2,\ldots, d_n\}$ by decreasing order of slope
Stack $S \leftarrow \emptyset$
$S$.push($d_1$)
$I_1 \leftarrow [-\infty, +\infty]$
$m \leftarrow 1$
for $i \leftarrow 2$ to $n$
$d \leftarrow S$.peek() [[examine the top element]]
$x_p \leftarrow$ FindIntersectionX($d,d_i$)
while $x_p \leq \mathrm{left}(I_m)$ [[found a redundant line]]
$S.$pop()
$m \leftarrow m-1$
$d \leftarrow S.$peek()
$x_p \leftarrow$ FindIntersectionX($d_i, d$)
$S.$push($d_i$)
$I_m \leftarrow [\mathrm{left}(I_m), x_p]$
$I_{m+1} \leftarrow [x_p, +\infty]$
associate $I_{m+1}$ with $d_i$
$m \leftarrow m+1$
Code của giả mã bằng C (code đầy đủ được cho ở cuối bài):
typedef struct{
double a;
double b;
} double_pair;
struct Stack_node {
int lineId;
struct Stack_node *next;
};
typedef struct Stack_node stack_node;
stack_node *S; // the stack
double_pair D[MAX_SIZE]; // the set of lines
double_pair I[MAX_SIZE]; // the set of intervals
double Q[MAX_SIZE]; // the set of query points
int ALines[MAX_SIZE];// the set of lines associated with intervals
int find_hull(int n){
qsort(D, n, sizeof(*D), slope_compare);
S = malloc( sizeof(stack_node));
S->lineId = -1; // head of the stack
push(0); // push d_1 to the stack;
I[0].a = -(double)INFTY;
I[0].b = (double)INFTY;
int m = 0;
ALines[m] = 0;
int i = 0;
stack_node *d;
double x = 0.0;
for(i = 1; i <n ;i++){
d = peek();
x = find_intersection_x(d->lineId, i);
while(x <= I[m].a){ // found a redundant line
pop(); // remove the redundant line
m--;
d = peek();
x = find_intersection_x(d->lineId, i);
}
push(i); // push d_i to the stack
I[m].b = x;
I[m+1].a = x;
I[m+1].b = (double)INFTY;
ALines[m+1] = i;
m++;
}
return m+1; // the number of intervals
}
double find_intersection_x(int x, int y){
double_pair d_a = D[x];
double_pair d_b = D[y];
return (d_b.b - d_a.b)/(d_a.a - d_b.a);
}
Thủ tục FindIntersectionX($d_1,d_2$) tìm hoành độ giao điểm của hai đường thẳng $d_1,d_2$ (chi tiết coi như bài tập cho bạn đọc). Ta có thể thấy vòng for của thuật toán tìm bao lồi có thời gian chạy $O(n)$ vì mỗi đường thẳng sẽ được xem xét tối đa hai lần: khi đưa vào stack và khi lấy ra khỏi stack. Nếu một đường thẳng bị lấy ra khỏi stack, nó sẽ không bao giờ được kiểm tra lại nữa. Do đó, thời gian của thuật toán tìm bao lồi chủ yếu dành để sắp xếp các đường thẳng theo slope và mất $O(n\log n)$.
Chứng minh Lemma 1 Giả sử chúng ta đã tìm được bao lồi của $D$, với mỗi điểm $p_\ell$ trên trục hoành, sử dụng tìm kiếm nhị phân để tìm interval $I \in {I_1,I_2,\ldots, I_m}$ sao cho $p_\ell \in I$. Giả sử $(d) \quad y= ax+b$ là đường thẳng tương ứng với $i$, ta sẽ có $q_\ell = a\cdot p_\ell +b$. Do đó, ta có thể tính $q_\ell$ trong thời gian $O(\log n)$. Nếu ta có $k$ điểm $p_1,p_2,\ldots, p_k$, ta có thể tìm $q_1,q_1,\ldots, q_k$ trong thời gian $O(k\log n)$. Chi thiết thuật toán như sau:
FindHull($d_1,d_2,\ldots, d_n$)
for $i \leftarrow 1$ to $k$
$I \leftarrow$ IntervalBinSearch($\{I_1,I_2,\ldots, I_m\}, p_i$)
$d \leftarrow$ the line associated with $I$
output $q_i = a\cdot p_i + b$ [[assuming $(d) \quad y= ax + b$]]
if $m = 1$
return $I_m$
else
$\ell = \lfloor m/2 \rfloor$
if $p \in I_\ell$
return $I_\ell$
else if $p > \mathrm{right}(I_\ell)$
return IntervalBinSearch($\{I_{\ell +1},\ldots, I_m\}, p$)
else
return IntervalBinSearch($\{I_1,\ldots, I_{\ell -1}\}, p$)
Code của giả mã bằng C:
void query(double points[], int k, int n){
int m = find_hull(n);
int i = 0, q = 0;
for(i = 0; i < k; i++){
q = interval_search(0, m-1, points[i]);
Q[i] = D[ALines[q]].a*points[i] + D[ALines[q]].b;
}
}
int interval_search(int x, int y, double p){
if(y <= x){
return x;
} else {
int mid = (x+y)/2;
if((I[mid].a <= p) && (p <= I[mid].b)){
return mid;
} else if (p > I[mid].b){
return interval_search(mid+1, y, p);
} else {
return interval_search(x, mid-1, p);
}
}
}
Quay trở lại bài toán quy hoạch động chúng ta đang xét. Nếu ta đặt $(d_j) \quad y_j = B[j]x + C[j],$ $1 \leq j \leq n$ thì $C[i]$ chính là tung độ $q$ nhỏ nhất thuộc một trong các đường thẳng $d_1,d_2, \ldots, d_{i-1}$ tương ứng với hoành độ $p = A[i]$. Do đó, chúng ta có thể áp dụng ý tưởng của kĩ thuật bao lồi vào bài toán quy hoạch động này.
Một số trường hợp đặc biệt
Trường hợp 1: thuật toán $O(n\log n)$
Trước tiên chúng ta xét trường hợp mảng $B[1,2,\ldots,n]$ thỏa mãn tính chất sau:
$B[1] \geq B[2] \geq \ldots \geq B[n]$
Giả sử $B[1] \geq B[2] \geq \ldots \geq B[n]$ tuy khá mạnh nhưng rất nhiêu bài toán chúng ta gặp sẽ thảo mãn điều kiện này. Với giả sử này, slope cuả các đường thẳng $d_1,d_2,\ldots, d_n$ đã được sắp xếp theo thử tự giảm dần. Do đó ta tiết kiệm được bước sắp xếp trong thuật tóan FindHull($d_1,d_2,\ldots, d_n$). Để giả mã đơn giản, ta giả sử trong $B[i] \not= B[i+1]$ với mọi $i$. Nếu tồn tại $i$ sao cho $B[i] = B[i+1]$, ta chỉ cần thay đổi một chút trong giả mã dưới đây mà không thay đổi thời gian tính của thuật toán. Thuật toán chi tiết như sau :
$d_1 \leftarrow B[1]x + C[1]$
Stack $S \leftarrow \emptyset$
$S$.push($d_1$)
$I_1 \leftarrow [-\infty, +\infty]$
$m \leftarrow 1$
for $i \leftarrow 2$ to $n$
$I \leftarrow$ IntervalBinSearch($\{I_1,I_2,\ldots, I_m\}, A[i]$)
$d \leftarrow$ the line associated with $I$
$C[i] \leftarrow a\cdot A[i] + b$ [[assuming $(d) \quad y= ax + b$]]
$d_i \leftarrow B[i]x + C[i]$
$d \leftarrow S$.peek() [[examine the top element]]
$x_p \leftarrow$ FindIntersectionX($d,d_i$)
while $x_p \leq \mathrm{left}(I_m)$ [[found a redundant line]]
$S.$pop()
$m \leftarrow m-1$
$d \leftarrow S.$peek()
$x_p \leftarrow$ FindIntersectionX($d_i, d$)
$S.$push($d_i$)
$I_m \leftarrow [\mathrm{left}(I_m), x_p]$
$I_{m+1} \leftarrow [x_p, +\infty]$
associate $I_{m+1}$ with $d_i$
$m \leftarrow m+1$
return $C[n]$
Phân tích thời gian Mỗi bước tìm kiếm nhị phân ta phải trả thời gian $O(\log m)$. Do đó tổng thời gian của thuật toán là $O(n\log m) = O(n\log n)$. Thực tế, số đường thằng thuộc bao lồi thường nhỏ hơn rất nhiều số đường thẳng đầu vào, do đó, thuật toán này khá nhanh.
Trường hợp 2: thuật toán $O(n)$
Trong trường hợp này, chúng ta giả sử cả hai mảng $A[1,2,\ldots,n]$, $B[1,2,\ldots,n]$ thỏa mãn:
$A[1] \leq A[2] \leq \ldots \leq A[n] \\ B[1] \geq B[2] \geq \ldots \geq B[n]$
Bài toán chặt cây xét dưới đây sẽ thoả mãn trường hợp này. Vì $A[i-1] \leq A[i]$, khi tìm kiếm interval chứa $A[i]$, chúng ta không cần phải xét các interval nằm bên trái $A[i-1]$. Do đó, thay vì thực hiện tìm kiếm nhịn phân trong thuật toán FastDynamic, ta tìm kiếm tuyến tính với điểm tìm kiếm bắt đầu là interval chứa $A[i-1]$. Trong giả mã dưới đây, ta giả sử $B[i] \not= B[i+1]$ với mọi $i$.
$d_1 \leftarrow B[1]x + C[1]$
Stack $S \leftarrow \emptyset$
$S$.push($d_1$)
$I_1 \leftarrow [-\infty, +\infty]$
associate $I_1$ with $A[1]$
$m \leftarrow 1$
for $i \leftarrow 2$ to $n$
$I_j \leftarrow$ the interval associated with $A[i-1]$
$I\leftarrow \emptyset$
$\ell \leftarrow j-1$
while $I = \emptyset$
$\ell \leftarrow \ell +1$
if $A[i] \in I[\ell]$
$I \leftarrow I[\ell]$
$d \leftarrow$ the line associated with $I$
$C[i] \leftarrow a\cdot A[i] + b$ [[assuming $(d) \quad y= ax + b$]]
$d_i \leftarrow B[i]x + C[i]$
$d \leftarrow S$.peek() [[examine the top element]]
$x_p \leftarrow$ FindIntersectionX($d,d_i$)
while $x_p \leq \mathrm{left}(I_m)$ [[found a redundant line]]
$S.$pop()
$m \leftarrow m-1$
$d \leftarrow S.$peek()
$x_p \leftarrow$ FindIntersectionX($d_i, d$)
$S.$push($d_i$)
$I_m \leftarrow [\mathrm{left}(I_m), x_p]$
$I_{m+1} \leftarrow [x_p, +\infty]$
if $\ell < m$ [[the line associated with $I_\ell$ is not redundant]]
associate $I_\ell$ with $A[i]$
else
if $A[i] \in I_m$
associate $I_m$ with $A[i]$
else
associate $I_{m+1}$ with $A[i]$
associate $I_{m+1}$ with $d_i$
$m \leftarrow m+1$
return $C[n]$
Phân tích thời gian Ta chỉ cần phân tích thời gian của thủ tục tìm kiếm interval $I_\ell$ sao cho $A[i]\in I_\ell$ trong mỗi vòng lắp. Trong trường hợp tồi nhất, chúng ta phải duyệt qua toàn bộ tập các interval. Tuy nhiên, nếu chúng ta đã duyệt qua $K$ intervals, thì các vòng lặp sau, chúng ta không cần phải xét các interval này và các đường thằng tương ứng với chúng nữa. Do đó, tổng thời gian tìm kiếm interval của thuật toán chỉ là $O(n)$ và đó cũng là thời gian của cả thuật toán.
Trường hợp tổng quát
Trong trường hợp tổng quát, chúng ta có thể áp dụng thủ tục IntervalBinSearch($\{I_1,I_2,\ldots, I_m\}, A[i]$) để tìm kiếm interval chứa $A[i]$. Do đó chúng ta mất $O(\log n)$ cho việc tìm kiếm interval cho mỗi vòng lặp của thuật toán. Vấn đề còn lại chỉ là làm sao để cập nhật bao lồi trong mỗi bước khi ta thêm đường thằng $(d) \quad B[i]x + C[i]$. Ở đây mình chỉ đưa ra ý tưởng để thực hiện cập nhật bao lồi trong thời gian trung bình $O(\log n)$ mỗi bước. Do đó, tổng thời gian của thuật toán là $O(n\log n)$.
Ý tưởng của thuật toán dựa trên thuật toán Sweep Line (mình sẽ giới thiệu trong bài viết sau). Ta sẽ luôn luôn lưu các đường thẳng trong bao lồi theo thứ tự tăng dần của slope sau mỗi bước. Gọi $d_1,d_2,\ldots, d_{m}$ là cá đường thẳng trong bao lồi sau bước $i-1$. Để đơn giản ta giả sử không có hai đường thẳng nào song song. Trong trường hợp có hai đường song song, ta chỉ cần thay đổi thuậtt toán một chút. Bây giờ ta có thêm đường thằng $(d_i) \quad y = a_ix + b_i$, để cập nhật bao lồi trước hết ta thực hiện tìm kiếm nhị phân để tìm ra đường thẳng $(d_j) \quad y = a_jy + b_j$ sao cho $a_j > a_i > a_{j+1}$. Trường hợp $j = m$, ta cập nhật như trong trường hợp đặc biệt 1. Do đó, ta có thể giả sử $j < m$. Gọi $(x_p, y_p)$ là giao điểm của $d_j$ và $d_{j+1}$. Xét đường thằng $(d) \quad y = a_ix + b_i'$ có cùng slope $a_i$ với $d_i$ và đi qua $(x_p,y_p)$. Nếu $b_i' \leq b_i$ (đường thằng $d_i$ có màu xanh trong hình dưới đây), đường thằng $d_i$ là dư thừa và do đó, bao lồi không thay đổi khi ta thêm đường thằng $d_i$. Nếu $b_i' > b_i$, ta lần lượt "dịch" đường thằng $d$ về phía $d_i$ cho đến khi $d$ trùng với $d_i$. Nếu có sự "thay đổi" của bao lồi (các đường thẳng màu nâu), ta cập nhật lại bao lồi. Ta có thể nhận thấy trong trường hợp tồi nhất, ta phải dịch đường thằng $d$ $O(m)$ lần. Tuy nhiên, nếu ta dịch $K$ lần, thì mỗi lần dịch ta sẽ "loại bớt" được $K$ đường thẳng dư thừa và các đường thẳng này sẽ không được xét đến trong các vòng lặp tiếp theo. Do đó, thời gian trung bình mỗi bước lặp vẫn là $O(\log n)$
Ví dụ áp dụng
Bài toán sau được lấy từ codeforces. Bạn đọc có thể áp dụng phương pháp mình giới thiệu ở đây để nộp bài:
Ví dụ: $n = 6, A[1,2,\ldots,6] = \{1,2,3,10,20,30\}$, $B[1,2,\ldots, 6] = \{6,5,4,3,2,0\}$. Nếu bạn chặt cây theo thứ tự
$1 \rightarrow 2 \rightarrow 6 \rightarrow 3 \rightarrow 4 \rightarrow 5$
chi phí bạn phải trả là $2\times 6 + 30\times 5 + 3 \times 0 + 10 \times 0 + 20 \times 0 = 162$. Nếu bạn chặt theo thứ tự:
$1 \rightarrow 3 \rightarrow 6 \rightarrow 2 \rightarrow 4 \rightarrow 5$
chí phí bạn phải trả là $3 \times 6 + 30\times 4 + 4 \times 0 + 10 \times 0 + 20 \times 0 = 138$
Từ ví dụ trên, ta có nhận xét sau:
Do đó, bài toán trên có thể được quy về bài toán: tìm chi phí nhỏ nhất để chặt cây thứ $n$. Gọi $OPT = \{1 = i_1,i_2,\ldots, i_k\}$ là một thứ tự chặt cây tối ưu với $i_k = n$. Ta có bổ đề sau:
$1 = i_1 \leq i_2 \leq \ldots \leq i_k = n$
Chứng minh Giả sử bổ đề là sai, suy ra tồn tại $a$ nhỏ nhất sao cho sao cho $ 1 = i_1 \leq i_2 \leq i_a \geq i_{a+1}$. Xét dãy $OPT\setminus\{i_{a+1}\} = \{1 = i_1, i_2, \ldots, i_a, i_{a+2}\ldots, i_k = n \}$. Chi phí chặt cây theo thứ tự của dãy này nhỏ hơn $OPT$ (bạn đọc tự kiểm chứng), do đó, trái với giả thiết dãy $OPT$ là thứ tự chặt cây tối ưu.
Dựa vào Lemma 1, ta có thể phát triển thuật toán quy hoạch động như sau: Gọi $C[i]$ là chi phí nhỏ nhất để chặt cây thứ $i$ "sử dụng" các cây $1,2,\ldots, i-1$. Ta có công thức sau:
$C[i] = \begin{cases} B[i], & \mbox{if } i = 1 \\ \min_{j < i}\{C[j] + A[i]B[j] \} + B[j], & \mbox{otherwise } \end{cases} \qquad (1)$
Dễ thấy, chi phí để chặt cây thứ $n$ là $C[n]$. Để ý kĩ các bạn sẽ thấy công thức quy hoạch động để tính $C[n]$ với hai mảng $A,B$ là các mảng đã săp xếp tương ứng với trường hợp đặc biệt 2. Do đó, bạn có thể giải bài toán này trong thời gian $O(n)$.
Code bằng $C$ của thuật toán được cho ở cuối bài.
Code: bao lồi, tree-cutting, code đã accept trên codeforces,
Tham khảo
[1] http://wcipeg.com/wiki/Convex_hull_trick
[2] Brucker, P. (1995). Efficient algorithms for some path partitioning problems. Discrete Applied Mathematics, 62(1-3), 77-85.
[3] Codeforces forum: http://codeforces.com/blog/entry/8219
Recent Comments