Phương pháp gradient liên hợp giải hệ phương trình đại số tuyến tính - Bùi Thị Thanh Xuân

Bùi Thị Thanh Xuân và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ 116 (02): 117 - 121 117 PHƢƠNG PHÁP GRADIENT LIÊN HỢP GIẢI HỆ PHƢƠNG TRÌNH ĐẠI SỐ TUYẾN TÍNH Bùi Thị Thanh Xuân*, Dƣơng Thị Nhung Trường Đại học Công nghệ thông tin và Truyền thông – ĐH Thái Nguyên TÓM TẮT Phƣơng pháp Gradient liên hợp đƣợc Hestenses và Stiefel nêu ra đầu tiên vào những năm 1950 để giải hệ phƣơng trình đại số tuyến tính (PTĐSTT). Vì việc giải một hệ phƣơng trình tuyến tính tƣơng đƣơng với việc tìm cực tiểu của một hàm toàn phƣơng xác định dƣơng nên năm 1960 Fletcher – Reeves đã cải biên và phát triển nó thành phƣơng pháp Gradient liên hợp cho cực tiểu không ràng buộc. Phƣơng pháp Gradient liên hợp (CG) chỉ cần tới đạo hàm bậc nhất nhƣng làm tăng hiệu quả và độ tin cậy của thuật toán. Trên cơ sở đó, bài báo nghiên cứu thuật toán CG và cài đặt thử nghiệm trên Matlab. Từ khóa: hệ phương trình đại số tuyến tính, phương pháp lặp, phương pháp Gradient liên hợp, CG, dạng toàn phương. TỔNG QUAN* Xét hệ phƣơng trình đại số tuyến tính Ax b trong đó A là ma trận vuông cấp n, x và b là véc tơ n chiều 11 12 1 1 1 21 22 2 2 2 1 2 ... ... ... ... ... ... ... ... ... n n n n nn n n a a a x b a a a x b a a a x b Hệ phƣơng trình đại số tuyến tính xuất hiện trong rất nhiều lĩnh vực nhƣ trong kinh tế, thống kê, hệ thống điện, xử lý ảnh, tối ƣu hóa, giải số các phƣơng trình vi phân,... với kích thƣớc của bài toán n có thể là 2 hoặc đến hàng chục triệu. Do đó một yêu cầu cần thiết là cần có các phƣơng pháp hiệu quả để giải hệ đại số tuyến tính nói trên. Các nhà Toán học đã nghiên cứu các phƣơng pháp giải hệ ĐSTT và phân loại thành 2 nhóm phƣơng pháp giải: phƣơng pháp trực tiếp (phƣơng pháp cho ta nghiệm đúng của hệ sau một số hữu hạn các phép tính) và phƣơng pháp lặp (phƣơng pháp xây dựng một dãy vô hạn các xấp xỉ k x mà giới hạn của nó là nghiệm đúng của hệ) Các nhà toán học cũng đƣa ra sự hạn chế của các phƣơng pháp trực tiếp. Chẳng hạn nhƣ phƣơng pháp Cramer về mặt lý thuyết là giải đƣợc hệ với n tùy ý, trong thực tế nếu ta có hệ * Tel: 0906 062458, Email: bttxuan@ictu.edu.vn n ẩn số phƣơng pháp Cramer cần n!(n+1)(n-1) phép tính số học, khi đó máy tính hiện đại không thể tính đƣợc với n =20! Hay với phƣơng pháp Gauss đòi hỏi 3 2 3 n phép tính và hơn nữa thuật toán Gauss không ổn định. Nhƣ vậy để giải số các bài toán có kích thƣớc lớn cần dùng đến các phƣơng pháp lặp.Các phƣơng pháp lặp có thể kể đến nhƣ phƣơng pháp Jacobi, Gauss-Seidel, phƣơng pháp độ dốc lớn nhất (Steepest descent), ...đặc biệt là phƣơng pháp Gradient liên hợp (Conjugate Gradient method) tỏ ra hiệu quả khi ma trận hệ số là ma trận đối xứng, xác định dƣơng. Trong không gian nR ta đƣa vào một tích vô hƣớng mới đƣợc xác định bởi công thức: , , A x y Ax y , trong đó .,. là tích vô hƣớng thông thƣờng cho bởi 1 , n i i i x y x y . Tích vô hƣớng mới này đƣợc gọi là tích năng lượng và chuẩn cảm sinh bởi nó đƣợc gọi là chuẩn năng lượng. Nhƣ vậy , A x Ax x Xét quá trình lặp ( 1) ( ) ( ) , 0,1,2... k k kx xB Ax b k (*) với 0 x bất kỳ. Định lý Samarski: Giả sử A là ma trận đối xứng xác định dƣơng và B là ma trận xác định Bùi Thị Thanh Xuân và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ 116 (02): 117 - 121 118 dƣơng. Khi đó nếu 2 B A thì phƣơng pháp lặp (*) hội tụ theo chuẩn năng lƣợng của A, tức là lim 0 k k A x x . PHƢƠNG PHÁP ĐỘ DỐC LỚN NHẤT (STEEPEST DESCENT) Cho hệ phƣơng trình đại số tuyến tính Ax b với giả thiết A là ma trận đối xứng xác định dƣơng 0, 0Tx Ax x . Bài toán trên tƣơng đƣơng với bài toán cực tiểu hóa dạng toàn phƣơng: min nx x  , trong đó dạng toàn phƣơng 1 2 T Tx x Ax b x Gradient của dạng toàn phƣơng: 1 2 , ,.., T n x x x x x x x . Gradient là trƣờng véc tơ hƣớng theo hƣớng giảm nhanh nhất của x tại điểm x đã cho 1 1 2 2 Tx A x Ax b . Nếu A đối xứng thì x Ax b . Cho gradient bằng 0 ta thu đƣợc nghiệm của phƣơng trình Ax b . Nghiệm của phƣơng trình là điểm tới hạn của hàm . Nếu A vừa đối xứng vừa xác định dƣơng thì nghiệm của phƣơng trình Ax b là điểm cực tiểu toàn cục của x . Do x đạt cực trị khi gradient 0x Ax b nên bài toán tìm cực trị tƣơng đƣơng với việc giải hệ phƣơng trình đại số tuyến tính Ax b . Ta biết rằng gradient là hƣớng hàm tăng nhanh nhất. Nhƣ thế muốn đi đến cực tiểu ta cho x, tính gradient và tìm theo hƣớng ngƣợc lại cho đến khi hàm không giảm nữa. Phƣơng pháp độ dốc lớn nhất (steepest desceent) thực hiện nhƣ sau: - Xuất phát từ 0x tạo ra dãy 1 2, ,...x x tiến gần đến nghiệm x. - Tại mỗi bƣớc chọn hƣớng giảm nhanh nhất – hƣớng ngƣợc lại với gradient của tại ix : i ix b Ax Sai số i ie x x . Độ lệch i ir b Ax . Nhƣ vậy ;i i i ir Ae r x độ lệch là hƣớng giảm nhanh nhất của hàm x . Thuật toán Steepest Descent Cho 0x Tính 0 0r b Ax Lặp k =1,2,... T k k k T k k r r r Ar 1k k k kx x r 1 1k kr b Ax cho tới khi hội tụ. Phƣơng pháp Steepest descent lặp theo hƣớng của các bƣớc trƣớc nên hội tụ không nhanh, sẽ hiệu quả hơn nếu trong mỗi bƣớc của thuật toán khi ta trƣợt theo hƣớng nào đó thì ta trƣợt đến chỗ thuật toán dừng tại bƣớc lặp cuối cùng. PHƢƠNG PHÁP GRADIENT LIÊN HỢP Phƣơng pháp Gradient liên hợp là phƣơng pháp tìm cực tiểu hóa của dạng toàn phƣơng: * arg min 1 2 nx R T T x x x x Ax b x (1) trong đó n nA R là ma trận đối xứng xác định dƣơng và nb R là véc tơ cho trƣớc. Rõ ràng ta có: 2,x Ax b A (2) Từ (2) ta thấy rằng x* cũng chính là nghiệm của phƣơng trình tuyến tính Ax b . Từ định lý Taylor cho g t y tz chúng ta thu đƣợc với mọi t  và với mọi ny  và nz  ta có: 2 2 T Tty tz y t y z z Az (3) Phƣơng pháp tìm kiếm theo hƣớng: Cho một xấp xỉ jx của phƣơng án tối ƣu x* và một véc tơ hƣớng jp , phƣơng pháp CG sẽ đƣa ra một xấp xỉ tiếp theo 1jx thông qua 2 bƣớc: Bùi Thị Thanh Xuân và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ 116 (02): 117 - 121 119 Tìm argminj j jx p Đặt 1j j j jx x p Giả sử rằng 0x là xấp xỉ ban đầu, sau đó, áp dụng k bƣớc lặp sẽ trả về kết quả là một dãy lặp 1 0 k j j x . Từ (2) và (3) chúng ta tìm đƣợc T j j j T j j p r p Ap trong đó j j jr x Ax b , ở đây jr thƣờng đƣợc gọi là độ lệch (residual). Định nghĩa 1. Tập hƣớng 1 0 k j j p là tập hƣớng liên hợp nếu , 0Ti j j iAp p p Ap i j Ký hiệu: 0 1 1 0, ,.., W 0k nW span p p p Siêu phẳng: 0 0W w ,w W n k k k k kU x z R z x 0 0U x Bổ đề 1. Cho ip liên hợp, 1j j j jx x p . Cho j jr Ax b , 0x cho trƣớc. argminj j jf x p . Khi đó , , ( )i i i ip r p f y y U Chứng minh: Nhận xét rằng i ix U . Lấy Wi i iy U x y và , 0i iA x y p , , , , 0 , , ( ) i i i i i i i i i i i i p r y p Ax b y p Ax b Ay b p A x y p r p y y U Định lý 1. Cho ip liên hợp. Khi đó arg min 1 * j j y U x y j k Chứng minh: Ta chứng minh định lý bằng phƣơng pháp quy nạp 1 11: argmin y U k x y . Giả sử (*) đúng với k i , tức là arg min 1 j j y U x y j i . Cần chứng minh (*) đúng với 1k i , tức là cần chứng minh 1 1 argmin i i y U x y trong đó 1i i i ix x p . Theo định nghĩa 1iU , 1ix U có thể viết ix y p trong đó , iy U Áp dụng (3) và bổ đề 1 ta có: 2 2 2 2 i T T i i i T T i i i x y p y p y p Ap y p y p Ap Ta có: 1 2 min min min 2ni i T T i i i i x U y U x y p r p Ap  . Vế phải đạt cực tiểu khi iy x và T i i i T i i p r p Ap . Vế trái đạt cực tiểu chính xác khi 1i i i ix x p . Nhƣ vậy để xây dựng thuật toán Gradient liên hợp ta phải: Xây dựng công thức truy toán sinh các hƣớng liên hợp ip Rút gọn một số công thức tính toán Sử dụng các hƣớng liên hợp tính toán các ix Bổ đề 2. Cho 0 0 0 0,p r r Ax b 1 0 , , k k j k k j j j j Ar p p r p Ap p Kết luận: , 0, 0m jAp p m j k Chứng minh bằng phƣơng pháp quy nạp: Với 1k kiểm tra trực tiếp ta thấy 0 1, 0Ap p . Giả sử với k i , hệ véc tơ 0 i j j p là liên hợp từng đôi một. Chúng ta cần chứng minh rằng 1, 0, 0m iAp p m i . Đặt m i , sau đó ta có: 1 1 1 0 , , , , , i i j m i m k m j j j j Ar p Ap p Ap r Ap p Ap p Bùi Thị Thanh Xuân và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ 116 (02): 117 - 121 120 1 1 , , , 0 , i m m k m m m m Ar p Ap r Ap p Ap p Bổ đề 3. Giả sử 0 k j j p đƣợc xây dựng theo Bổ đề 2. Khi đó 0 1 1, ,..,k kW span r r r , 0, 0j mr r j m k , , 0j k k kr p r r j k kp thỏa mãn 1 1 1 1 1 , , , k k k k k k k k k r r p r p r r Thuật toán CONJUGATE GRADIENT - Cho xấp xỉ ban đầu 0x - Đặt 0 0 0 0,r Ax b p r k = 0 For k = 1: k max 1 1 1 1 k k k k k k k k k k k k k k k k k k k k k k p r p Ap x x p r r Ap p Ar p Ap p r p end Ta thấy rằng phƣơng pháp CG nếu không có sai số tính toán thì sẽ dừng sau không quá n bƣớc và thu đƣợc nghiệm đúng. Tuy nhiên trong quá trình giải có sai số tính toán nên phƣơng pháp CG coi nhƣ phƣơng pháp lặp để tìm nghiệm gần đúng. Đánh giá sai số của phƣơng pháp Gradient liên hợp: 0 1 * 2 * 1 l l A A x x x x với 1 nA là số điều kiện của A Thử nghiệm trên MATLAB n = 100; R = randn(n); for npot = 0.1:0.4:2 A = (R'*R)^npot; lamb = eig(A); kappa = max(lamb)/min(lamb); xsol = rand(n,1); b = A*xsol; Norm2 = sqrt(xsol'*xsol); NormA = sqrt(xsol'*A*xsol); x = zeros(size(b)); g = A*x-b; p = -g; for i = 1:n Ap = A*p; alfa = -(p'*g)./(p'*Ap); x = x + alfa*p; g = g + alfa*Ap; beta = (g'*Ap)./(p'*Ap); p = -g + beta*p; err2(i) = sqrt((x-xsol)'*(x-xsol))/Norm2; errA(i) = sqrt((x-xsol)'*A*(x-xsol))/NormA; end subplot(1,2,1); plot(lamb,zeros(size(lamb)),'o'); Tittel = ['Eigenvalues, n=' num2str(n)]; title(Tittel); subplot(1,2,2); semilogy(1:n, err2, 1:n, errA); legend('chuan 2', 'chuan A'); xlabel('So lan lap'); ylabel('Sai so'); Tittel = ['npot = ',num2str(npot),'\kappa = ',num2str(kappa)]; title(Tittel); pause; end Hình 1. Ma trận có điều kiện tốt, hội tụ nhanh 0.5 1 1.5 2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Eigenvalues, ndim=100 0 50 100 10 -16 10 -14 10 -12 10 -10 10 -8 10 -6 10 -4 10 -2 10 0 So lan lap sa i s o npot = 0.1 = 3.5136 chuan 2 chuan A Bùi Thị Thanh Xuân và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ 116 (02): 117 - 121 121 Hình 2. Ma trận có điều kiện trung bình Hình 3. Ma trận với điều kiện xấu, tốc độ hội tụ chậm KẾT LUẬN Trong bài báo này chúng tôi đã giới thiệu đến một số phƣơng pháp giải hệ đại số tuyến tính, trên cơ sở đó nghiên cứu đến một phƣơng pháp hiệu quả giải hệ ĐSTT, đó chính là phƣơng pháp Gradient liên hợp (CG). Ngoài ra chúng tôi cũng đã mô phỏng thuật toán CG bằng Matlab để so sánh sự hiệu quả của phƣơng pháp trong các điều kiện của ma trận hệ số A. TÀI LIỆU THAM KHẢO 1. A. van der Sluis and H. A. van derVorst (1986), The Rate of Convergence of Conjugate Gradients, Numerische Mathematik 48, no. 5, 543–560. 2. Magnus R. Hestenes and Eduard Stiefel (1952), Methods of Conjugate Gradients for Solving Linear Systems, Journal of Research of the National Bureau of Standards 49, 409–436. 3. R. Fletcher and M. J. D. Powell (1963), A Rapidly Convergent Descent Method for Minimization, Computer Journal 6,163–168. 4. Đặng Quang Á (2009), Giáo trình Phương pháp số, Nxb Đại học Thái Nguyên. 5. Nguyễn Hoàng Hải (2006), Lập trình Matlab và ứng dụng, Nxb Khoa học kỹ thuật SUMMARY THE CONJUGATE GRADIENT METHOD FOR SOLVING THE SYSTEM OF LINEAR EQUATIONS Bui Thi Thanh Xuan * , Duong Thi Nhung College of Information & Communication Technology – TNU The conjugate gradient method (CG) has been first introduced in 1952 by M. R. Hestenes and E. Stiefel. It has been the subject of intense analysis for more than 50 years. The method was originally consider to be direct method for linear equations, but its favorable properties as an iterative method was soon realized, and it was later generalized to more general optimization problems. It provides a very effective way to optimize large, deterministic systems by gradient descent. In this paper, we introduce about CG and experiment with CG method using Matlab. Key words: iterative method, linear equations, conjugate gradient method, CG, quadratic functions. Ngày nhận bài:23/09/2013; Ngày phản biện:04/10/2013; Ngày duyệt đăng: 26/02/2014 Phản biện khoa học: TS. Trương Hà Hải – Trường ĐH Công nghệ Thông tin & Truyền thông - ĐHTN * Tel: 0906 062458, Email: bttxuan@ictu.edu.vn 0 5 10 15 20 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Eigenvalues, ndim=100 0 50 100 10 -15 10 -10 10 -5 10 0 So lan lap sa i s o npot = 0.5 = 535.4972 chuan 2 chuan A 0 1000 2000 3000 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Eigenvalues, ndim=100 0 50 100 10 -3 10 -2 10 -1 10 0 So lan lap sa i s o npot = 1.3 = 12438519.0126 chuan 2 chuan A