Các phương pháp bậc hai hội tụ tốt hơn so với các phương pháp bậc một và được sử dụng nhiều trong các phần mềm giải bài toán tối ưu. Tuy nhiên, trong một số ứng dụng như trong máy học, mạng truyền thông thì phương pháp bậc hai không được áp dụng nhiều như các phương pháp bậc một do khối lượng tính toán trong một lần cập nhật phức tạp. Ví dụ như phương pháp Newton phải tính nghịch đảo của đạo hàm bậc hai. Đồng thời, với các bài toán kích thước lớn, phương pháp bậc hai khó thực hiện tính toán phân tán và song song. Hai phần tiếp theo đây sẽ trình bày hai phương pháp bậc hai phổ biến là Newton và log barrier.
Xét bài toán tối ưu \[\begin{aligned} \min_x f[x],\end{aligned}\] trong đó \[f[x]\] là hàm lồi và có đạo hàm bậc hai trên mọi điểm. Cập nhật theo phương pháp Newton có dạng sau: \[\begin{aligned} x^{[k]} = x^{[k – 1]} – [\nabla^2 f[x^{[k – 1]}]]{-1} \nabla f[x{[k – 1]}]. \end{aligned}\]
Nhắc lại, cập nhật theo gradient descent là chọn điểm \[x^{[k]}\] để cực tiểu hàm toàn phương xấp xỉ hàm \[f[x]\] trong lân cận điểm \[x^{[k-1]}\]: \[\begin{aligned} f[y] \approx f[x^{[k-1]}] + \nabla^Tf[x^{[k-1]}] [y – x] + \frac{1}{2t}\|y – x^{[k-1]}\|2_2,\end{aligned}\] và \[\begin{aligned} x{[k]} = arg\min_y f[x^{[k-1]}] + \nabla^Tf[x^{[k-1]}] [y – x] + \frac{1}{2t}\|y – x^{[k-1]}\|2_2.\end{aligned}\] Hàm toàn phương trong phương pháp Newton là xấp xỉ Taylor, tốt hơn trong GD: \[\begin{aligned} f[y] \approx f[x{[k-1]}] + \nabla^Tf[x^{[k-1]}] [y – x] + \frac{1}{2}[y-x^{[k-1]}]T \nabla^2 f[x{[k-1]}][y-x^{[k-1]}],\end{aligned}\] do đó, \[\begin{aligned} x^{[k]} = arg\min_y f[x^{[k-1]}] + \nabla^Tf[x^{[k-1]}] [y – x] + \frac{1}{2}[y-x^{[k-1]}]T \nabla^2 f[x{[k-1]}][y-x^{[k-1]}].\end{aligned}\]
Ví dụ áp dụng phương pháp GD và Newton để tìm cực tiểu hàm số: \[\begin{aligned} f[x_1, x_2] = \frac{10x_1^2 + x_2^2}{2} + 5\log[1 + e^{-x_1-x_2}].\end{aligned}\] Gradient và Hessian được cho bởi: \[\begin{aligned} \nabla f[x_1, x_2] = \begin{bmatrix} 10x_1 – \frac{5e^{-x_1-x_2}}{1 + e^{-x_1-x_2}}\\ x_2 – \frac{5e^{-x_1-x_2}}{1 + e^{-x_1-x_2}} \end{bmatrix} \qquad \nabla^2 f[x_1, x_2] = \begin{bmatrix} 10 + \frac{e^{-x_1-x_2}}{[1 + e^{-x_1-x_2}]2} & \frac{e{-x_1-x_2}}{[1 + e^{-x_1-x_2}]2}\\ \frac{e{-x_1-x_2}}{[1 + e^{-x_1-x_2}]2} & 1 + \frac{e{-x_1-x_2}}{[1 + e^{-x_1-x_2}]^2} \end{bmatrix} \end{aligned}\]
Hình 1 cho thấy phương pháp Newton có số vòng lặp nhỏ hơn GD rất nhiều. Đặc biệt tìm cực tiểu hàm toàn phương, phương pháp Newton hội tụ chỉ trong 1 bước lặp.
Lưu ý là phương pháp Newton thuần túy không dùng step size. Trong thực tế, ta có thể thực hiện cập nhật Newton với step size:\[\begin{aligned} x^{[k]} = x^{[k – 1]} – t_k[\nabla^2 f[x^{[k – 1]}]]{-1} \nabla f[x{[k – 1]}], \end{aligned}\]trong đó step size \[t_k\] có thể được xác định bằng backtracking line search.
Bạn đọc có thể viết chương trình áp dụng phương pháp Newton cho linear regression hay logistic regression dễ dàng dựa theo những code Python ví dụ đã trình bày trước đây.
Trong chương trình phổ thông, khoảng lớp 5 học sinh được học cách giải phương trình bậc nhất 1 ẩn [qua một lớp bài toán được biết với cái tên nổi tiếng là tìm x]. Khoảng cấp 2, HS được học cách giải phương trình bậc hai một ẩn và hệ phương trình bậc nhất hai ẩn.
Tốt nghiệp cấp 3, nếu đi học tiếp Cử nhân Công nghệ thông tin [như người viết bài này], HS sẽ được học tiếp cách giải hệ phương trình tuyến tính nhiều ẩn dạng [với là ma trận ] trong môn Đại số tuyến tính. Giải thuật lúc đó được dạy thì người viết không còn nhớ rõ, nhưng đại khái là tìm cách sử dụng các phép biến đổi dòng trên ma trận một cách khéo léo để tìm nghiệm. Mặc dù phương pháp đó có thể cài đặt được thành thuật toán, nhưng vì dựa vào trick [và “quan sát” theo kiểu của con người] khá nhiều nên không thực sự straightforward [sau vài ngày thì người viết nhớ ra tên phương pháp là Phép khử Gaussian]
Người viết có dịp tìm hiểu phương pháp Newton để giải hệ phương trình phi tuyến lần đầu tiên trong quá trình thực hiện luận văn. Sau đó gần đây lại có cơ hội làm việc với các phương pháp khác như Conjugate Gradient, LSQR …
Một điều đáng ngạc nhiên là, theo quá trình học của bản thân người viết, những phương pháp giải xấp xỉ/chính xác các phương trình/hệ phương trình tuyến tính/phi tuyến đã không được hệ thống hoá một cách đầy đủ, cả trong chương trình Toán phổ thông lẫn Đại học không chuyên ngành Toán [mặc dù tư tưởng chung của các phương pháp này khá tự nhiên]. Vì vậy chuỗi bài này muốn hệ thống lại các phương pháp trên. Bài đầu tiên sẽ nói về phương pháp Newton-Raphson. Các thuật toán khác và một số vấn đề liên quan sẽ được trình bày trong các bài sau. Các phương pháp này là bước đi ban đầu vào optimization, lĩnh vực mà người viết mong muốn được tìm hiểu sâu hơn.
Phương pháp Newton-Raphson giải một [hệ] phương trình thực bằng cách xấp xỉ dần giá trị của nghiệm thông qua một dãy số [hội tụ về giá trị thật của nghiệm?]. Chúng ta sẽ bắt đầu với trường hợp hàm một biến, sau đó mở rộng cho trường hợp nhiều biến.
1. Trường hợp một biến
Cho hàm một biến , mục tiêu là tìm giá trị của thỏa phương trình . [Tất nhiên trong trường hợp là hàm bậc nhất theo x thì không còn gì phải bàn, ở đây ta muốn xét đến trường hợp tổng quát].
Phương pháp Newton xuất phát từ việc xấp xỉ giá trị của đạo hàm tại một điểm:
Do ta đang muốn tìm nghiệm sao cho , do đó với thì phương trình trên trở thành:
Tổng quát, phương pháp Newton lặp liên tục để tìm nghiệm ngày càng gần giá trị thực, với công thức nghiệm là:
với giá trị khởi tạo nào đó.
Hình 1 cho một cái nhìn trực quan hơn về lời giải của phương pháp này. Với càng lớn, giá trị của càng gần về vị trí giao với trục hoành của đồ thị hàm .
Tuy nhiên cũng giống như nhiều phương pháp tính xấp xỉ khác, khả năng tốc độ hội tụ của phương pháp Newton phụ thuộc nhiều vào giá trị “đoán mò” ban đầu . Giá trị càng gần nghiệm thực của phương trình thì phương pháp này càng hội tụ nhanh chóng. Trong trường hợp xấu [khi vô nghiệm, hoặc khi giá trị đoán mò không cùng nằm trên một đoạn biến thiên của hàm số] thì phương pháp sẽ khó hội tụ nhanh. Trong trường hợp phương trình không có nghiệm, phương pháp này cũng không thể phát hiện.
Một điểm hạn chế nữa là do yêu cầu tính tích phân của hàm số, phương pháp này sẽ khó khăn trong trường hợp việc tính đạo hàm trở nên phức tạp. Trong trường hợp này, người ta có thể xấp xỉ việc tính đạo hàm bằng cách tìm phương trình đường thẳng đi qua hai điểm lân cận của hàm số tại vị trí đang xét [để xấp xỉ hệ số góc của tiếp tuyến – vốn có thể tính chính xác bằng đạo hàm tại điểm đó].
Do sự phức tạp trong việc tính đạo hàm nên việc cài đặt một chương trình tổng quát cho thuật toán Newton-Raphson là tương đối khó. Tuy nhiên trong trường hợp là đa thức theo thì có thể dễ dàng cài đặt phương pháp này trong MATLAB như sau
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
function [ x ] = nrsolve[ ce, res, maxLoop ]
[ceX, ceY] = size[ce];
if [ceX > 1 && ceY == 1]
ce = ce';
[ceX, ceY] = size[ce];
`function`3 `function`4
`function`6`function`7`function`8
[ x ] = nrsolve[ ce, res, maxLoop ]`0[ x ] = nrsolve[ ce, res, maxLoop ]`1
`[ x ] = nrsolve[ ce, res, maxLoop ]`3
`[ x ] = nrsolve[ ce, res, maxLoop ]`5
`[ x ] = nrsolve[ ce, res, maxLoop ]`7
`[ x ] = nrsolve[ ce, res, maxLoop ]`9
` `1 ``2
if ``5
`6 `7
`[ x ] = nrsolve[ ce, res, maxLoop ]`3
`[ x ] = nrsolve[ ce, res, maxLoop ]`3
`[ceX, ceY] = size[ce];`3
`[ceX, ceY] = size[ce];`5
`[ceX, ceY] = size[ce];`7
`[ceX, ceY] = size[ce];`9
` `1
` `3 ``4
`[ceX, ceY] = size[ce];`7
`[ceX, ceY] = size[ce];`9
`if`0
`[ x ] = nrsolve[ ce, res, maxLoop ]`3
`if`4
`[ x ] = nrsolve[ ce, res, maxLoop ]`3
function `if`7
`if`9
`1 [ceX > 1 && ceY == 1]`2
if `[ceX > 1 && ceY == 1]`5
`6[ceX > 1 && ceY == 1]`7
`[ x ] = nrsolve[ ce, res, maxLoop ]`3
`[ x ] = nrsolve[ ce, res, maxLoop ]`3
`[ x ] = nrsolve[ ce, res, maxLoop ]`3
Trong chương trình trên, là vector chứa hệ số của đa thức. Chẳng hạn với đa thức thì . Tham số là độ lỗi cho phép, chương trình dừng khi sai số nhỏ hơn ngưỡng này. Trong trường hợp không tìm được nghiệm, chương trình dừng sau lần lặp.
Sử dụng hàm này, ta có thể tính nhanh nghiệm của phương trình [cũng chính là tìm xấp xỉ giá trị ] như sau:
1
2
3
``3
``4
` `6
Một điều khá thú vị là sử dụng cặp hàm , ta thấy hàm chạy không chậm hơn phép lũy thừa trong MATLAB nhiều lắm khi tính [xét tới việc ta đã “đoán mò” giá trị rất thô thiển].
1
2
3
4
5
6
7
8
9
10
11
12
13
14
``7
``8
``9
`ce = ce';`0
`ce = ce';`1
``9
``4
`ce = ce';`4`ce = ce';`5
`ce = ce';`6
``4
`ce = ce';`4`ce = ce';`5
`function`00
2. Trường hợp nhiều biến
Giả sử ta phải giải hệ phương trình biến sau:
Một trong những phương pháp giải hệ phương trình này là phương pháp Newton nhiều biến [Multivariate Newton-Raphson Method – MNRM]. Mở rộng từ phương pháp Newton cho hàm một biến, ta cũng bắt đầu từ việc đạo hàm từng phần mỗi phương trình :
Giống như trước, ta xấp xỉ như sau:
[1]
với là giá trị của nghiệm ở lần lặp thứ k. Ta mong muốn , do đó hệ phương trình [1] có thể viết dưới dạng ma trận:
[2]
trong đó là ma trận Jacobi kích thước tại bước lặp thứ
[3]
là vector thặng dư kích thước tại bước lặp thứ :
[4]
và cuối cùng
[5]
Do đó nghiệm của hệ tại bước thứ là
[6]
Sử dụng những phương trình này, thuật toán Newton tìm lời giải cho hệ phương trình rất đơn giản:
- Đoán mò nghiệm 2. Tính ma trận Jacobi và vector thặng dư.3. Giải phương trình [2] bằng các phương pháp trong đại số tuyến tính4. Cập nhật giá trị theo [6].5. Lặp lại bước 2 nếu lời giải chưa hội tụ.
3. Tổng kết
Phương pháp Newton-Raphson trong bài là một phương pháp hết sức đơn giản. Nó nằm trong họ các phương pháp Household mà nếu có dịp sẽ quay trở lại sau. Do hạn chế trong việc phải tính được đạo hàm của hàm số, phương pháp này thực sự có ý nghĩa khi biết trước dạng hàm. Trong trường hợp không biết trước dạng hàm thì cũng có thể áp dụng phương pháp này bằng cách tính xấp xỉ giá trị đạo hàm tại từng điểm [một trong những cách xấp xỉ là tính hiệu với “đủ gần” với ]. Phương pháp xấp xỉ đạo hàm này đặc biệt có ý nghĩa khi cài đặt trên máy tính.
Trong MATLAB, việc tìm nghiệm của hệ các phương trình được cài đặt sẵn trong hàm fslove[], fzero[]… Các hàm này cài đặt một số thuật toán khác để không phải tính đạo hàm như thuật toán Newton, hoặc tìm cách xấp xỉ giá trị đạo hàm. Đặc biệt hàm fsolve[] không cần giá trị đoán mò ban đầu.