Một số kĩ thuật cải tiến quy hoạch động III- Advanced Dynamic Programming

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:

Problem 1: Cho một bài toán với n phần tử khác nhau 1,2,\ldots,n. Phần tử thứ i có hai tính chất A[i], B[i]. Hàm mục tiêu của bài toán này có thể được biểu diễn thông qua bảng quy hoạch động một chiều C[1,2,\ldots, n] thỏa mãn hệ thức 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:

SlowDynamic(A[1,2,\ldots, n], B[1,2,\ldots, n]):
    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:

Problem 2: Cho một tập n đường thẳng D = \{d_1,d_2,\ldots, d_n\} trong đó (d_i) \quad y_i = a_ix + b_i, a_i \geq 0k điểm trên trục Ox: p_1, p_2, \ldots, p_k . Tìm q_1,q_2,\ldots, q_k trong đó mỗi giá trị q_\ell được định nghĩa như 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_kO(kn). Tuy nhiên, ta có thể tìm các điểm này nhanh hơn:

Lemma 1: Tồn tại một thuật toán tìm tất cả các điểm q_1,q_2 \ldots, q_k trong thời gian O((k+n)\log 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). lines-fig
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.

Lemma 2: Bao lồi của tập các đường thẳng D=\{d_1,d_2,\ldots, d_n\} cho trước có thể biểu diễn bằng m interval I_1,I_2, \ldots, I_mm đường thẳng tương ứng với mỗi đoạn d_1,d_2, \ldots, d_m với m \leq n. Hơn nữa, ta có thể tìm được biểu diễn đó trong thời gian O(n\log n)

 
Để đơ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 :

FindHull(d_1,d_2,\ldots, d_n):
    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:

Query(\{d_1,d_2,\ldots, d_n\}, \{p_1,p_2,\ldots, p_k\}):
    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]]

 

IntervalBinSearch(\{I_1,I_2,\ldots, I_m\}, p):
    if m = 1
    return I_m
    else
        \ell = \lfloor m/2 \rfloor
        if p \in I_\ell
            return I_\ell
        else if p  data-recalc-dims= \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:

Special Case 1:

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 :

FastDynamic(A[1,2,\ldots, n], B[1,2,\ldots, n]):
    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:

Special Case 2:

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.

FastFastDynamic(A[1,2,\ldots, n], B[1,2,\ldots, n]):
    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  data-recalc-dims= 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_jd_{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'  data-recalc-dims= 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)
fully-dynamic-fig

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:

Problem 1: n cái cây 1,2,\ldots,n trong đó cây thứ i có chiều dài A[i] \in \mathbb{N}. Bạn phải cắt tất cả các cây thành các đoạn có chiều dài 1. Bạn có một cái cưa máy (rởm), mỗi lần chỉ chặt được một đơn vị chiều dài, và mỗi lần sử dụng thì lại phải sạc pin. Chi phí để sạc pin phụ thuộc vào chi phí của cây đã bị cắt hoàn toàn (một cây bị cắt hoàn toàn có chiều dài 0). Nếu chỉ số lớn nhất của cây bị cắt hoàn toàn là i thì chi phí sạc là B[i], và khi bạn đã chọn một cây để cắt, bạn phải cắt nó hoàn toàn. Ban đầu cái cưa máy được sạc đầy pin. Giả sử a_1 = 1 \leq a_2 \leq \ldots \leq a_n b_1 \geq b_2 \geq \ldots \geq b_n = 0. Tìm một cách cưa cây với chi phí nhỏ nhất.

 
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:

Observation 1: Chi phí để chặt các cây sau khi đã chặt cây thứ n là 0.

 

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:

Lemma 1:

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ứ nC[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

Facebook Comments

Tags: , ,

Reply

Your email address will not be published. Required fields are marked *