Thuật toán QR tìm trị riêng

Bài 1: QR decomposition
Bài 2: Trị riêng và vector riêng
Bài 3: Bài toán tìm trị riêng
Bài 4: Thuật toán QR

Tôi dành buổi tối cuối cùng trước kì nghỉ lễ để hoàn tất nốt chuỗi bài về Trị riêng và vector riêng.

Như đã biết ở các bài trước, phương pháp kinh điển tìm trị riêng sử dụng phân tích QR, trong đó gồm 2 bước: bước đầu tiên biến ma trận ban đầu thành dạng Hessenger, bước thứ 2 biến ma trận Hessenberg về dạng tam giác trên, các trị riêng sẽ xuất hiện trên đường chéo của ma trận tam giác cuối cùng. Nếu ma trận ban đầu là ma trận Hermitian thì ma trận Hessenberg trở thành hermitian Hessenberg, tức là tridiagonal, và ma trận tam giác cuối cùng sẽ là ma trận đường chéo.

Bước thứ nhất biến về dạng Hessenberg đã được trình bày trong phần 3 của bài trước. Bài này sẽ nói về bước 2: biến ma trận Hessenberg thành ma trận tam giác. Để làm được việc đó, người ta dùng Thuật toán QR. Sở dĩ gọi là Thuật  toán QR là vì thuật toán này dựa trên phân tích QR của ma trận.

1. Thuật toán QR

Thuật toán QR “ngây thơ” nhìn đơn giản đến mức đáng ngạc nhiên:

Thuật toán 1: Thuật toán QR “ngây thơ”

A^{(0)}=A

for k = 1, 2…

Q^{(k)}R^{(k)} = A^{(k-1)}

A^{(k)} = R^{(k)}Q^{(k)}

Tất cả những gì cần làm là lấy phân tích QR của A, sau đó nhân kết quả theo thứ tự ngược lại. Dưới một số giả sử, thuật toán này hội tụ về dạng phân tích Schur của A: ma trận tam giác nếu A bất kì, và ma trận đường chéo nếu A là hermitian.

Khi lần đầu tiên nhìn thấy thuật toán này, thật khó tin là nó có thể hoạt động được. Do đó tôi đã thử viết một đoạn mã MATLAB ngắn sau để kiểm tra:


function [ V ] = eigen( A )
%EIGEN Summary of this function goes here
%   Detailed explanation goes here

At = A;
[sx, sy] = size(At);
while(abs(At(1, sx)) > eps || abs(At(sy, 1)) > eps)
[Q, R] = qr(At);
At = R*Q;
end
V = At;
end

Và thử nghiệm với 1 ma trận Hermitian làm ví dụ:


>> a = [2 1; 1 2];

>> v = eigen(a)

v =

3.0000    0.0000
0.0000    1.0000

Thật đáng ngạc nhiên, thuật toán hội tụ chỉ sau khoảng 10 lần lặp, và kết quả ma trận đường chéo gồm 2 phần tử 3, 1 chính là 2 trị riêng của ma trận.

Tất nhiên thuật toán QR được dùng trong thực tế không đơn giản như vậy. Thuật toán kinh điển hay được dùng có thay đổi như sau:

1. Trước khi bắt đầu vòng lặp, A được rút gọn về dạng Hessenberg, như trong bài trước.

2. Thay vì phân tích A^{(k)} tại mỗi bước, ta sẽ phân tích A^{(k)}-\mu^{(k)}I, với \mu^{(k)} là một ước lượng xấp xỉ nào đó cho trị riêng.

3. Khi tìm được trị riêng mới, ta chia nhỏ bài toán bằng cách chia A^{(k)} thành các ma trận con nhỏ hơn.

Cụ thể, thuật toán QR dùng trong thực tế có giải thuật như sau (để đơn giản, ta giả sử A là ma trận Hessenberg)

Thuật toán 2: Thuật toán QR thực tế

\left(Q^{(0)}\right)^\mathrm{T}A^{(0)}Q^{(0)} = A           A^{(0)} là dạng tridiagonal của A

for k=1,2,…

Chọn shift \mu^{(k)}           Chẳng hạn có thể chọn \mu^{(k)}=A^{(k-1)}_{mm}

Q^{(k)}R^{(k)}=A^{(k-1)}-\mu^{(k)}I

A^{(k)} = R^{(k)}Q^{(k)}+\mu^{(k)}I

Nếu có phần tử ngoài đường chéo A^{(k)}_{j,j+1} gần bằng 0 thì:

gán A^{(k)}_{j,j+1} = A^{(k)}_{j+1,j} = 0 để được

\left[\begin{array}{cc}A_1&\mathbf{0}\\\mathbf{0}&A_2\end{array}\right] = A^{(k)}

và áp dụng thuật toán QR cho A_1A_2.

Thuật toán 2 chính là thuật toán kinh điển tìm trị riêng cho ma trận trong những năm 1960. Mãi cho đến những năm 1990, sau khi thuật toán chia để trị được Cuppen công bố thì thuật toán  QR mới kết thúc vai trò của nó. Thuật toán chia để trị của Cuppen có chi phí tính toán chỉ bằng 1/2 thuật toán QR trong trường hợp cần tìm cả trị riêng và vector riêng, do đó nó đã trở thành thuật toán được cài đặt rộng rãi trong các thư viện toán học hiện nay như BLAS hay CUBLAS. (Tuy nhiên chuỗi bài này lại kết thúc ở thuật tóa QR là do nó đã bắt đầu bằng phân tích QR, và để cho trọn vẹn thì thuật toán QR cũng nên được điểm qua).

Một lí do nữa để trình bày thuật toán QR là do đặc tính thú vị của nó. Mặc dù nhìn có vẻ đơn giản nhưng để giải thích tại sao thuật toán này lại hội tụ về trị riêng và vector riêng thì thật không hề đơn giản. Mặc dù đã cố gắng bổ sung 2 bài về phép chiếu và bài toán  least squares – 2 công cụ để hiểu thuật toán QR – nhưng cuối cùng thì cá nhân tôi cũng cảm thấy thuật toán này quá phức tạp. Do vậy bài này sẽ dừng lại ở đây, mà không đi sâu vào tìm hiểu bản chất bên trong, những ý tưởng (rất đẹp) bên dưới thuật toán này. Hơn nữa đứng từ góc độ của một người làm ứng dụng, tôi cũng không có ý định đi quá sâu vào những chi tiết cụ thể như thế này. Tất nhiên thuật toán  QR có thể khiến những người đam mê vẻ đẹp của Toán học phải tò mò và hứng thú, nhưng tốt nhất là tôi sẽ dành nó cho sau này, khi có dịp.

Thực tế thì khi ra đời vào những năm 1960, thuật toán QR đã được coi là “viên ngọc quý” của cộng đồng numerical analysis, có lẽ là vì đặc tính thú vị của nó: bề ngoài rất đơn giản nhưng ẩn chứa nhiều tư tưởng đẹp của Đại số tuyến tính.

Khép lại chuỗi bài về phân tích QR và trị riêng của ma trận, còn rất nhiều nội dung mà tôi muốn viết trong kì nghỉ lễ. Tuy nhiên hãy cứ xem tình hình như thế nào…

Advertisements

8 comments

  1. Hi anh,

    Cảm ơn anh vì đã có những bài viết thú vị, mà đọc kỹ chưa chắc em đã hiểu được hết. Thực ra, ban đầu lúc search đến đây, em chỉ muốn có một hàm tính eigenvalues đơn giản nhất có thể (ko dùng các lib blas/lapack..) để đáp ứng yêu cầu là: chuyển hàm eig() từ Matlab sang C/C++. Có small chance nào a có thể giúp e dc ko ah?

    Em nghĩ đoạn mã ngắn trên đây của anh “có vấn đề”, vì nó dùng hàm qr(), cũng là một hàm built-in của Matlab, như vậy, v = eigen(a); chạy nhanh & giống v=eig(a) thì cũng chưa khẳng định được. (ví dụ eig() của Matlab cũng dùng qr() của Matlab, hoặc ngược lại). Liệu em có nghĩ sai?

    Cảm ơn anh

    1. Hi Thơ,
      Thuật toán “đơn giản” để tìm trị riêng là không thể, do đó không nên kì vọng một hàm như eig() trong C/C++. Bạn có thể dùng Intel MKL với các hàm liên quan được mô tả trong trang này: http://software.intel.com/sites/products/documentation/hpc/compilerpro/en-us/cpp/win/mkl/refman/lse/lse_intro.html
      Tùy vào ứng dụng cụ thể của bạn, có thể sẽ không cần dùng hết tất cả các hàm trong đó, nhưng ít nhất bạn phải dùng 2 hay 3 hàm (rút gọn về Hessenberg, tìm phân tích Schur v.v…) thì mới có được kết quả cuối cùng. Thực tế chi phí cho việc tìm trị riêng nhỏ hơn nhiều so với chi phí cho các hàm biến đổi ma trận.

      Ngoài ra còn có thuật toán của Cuppen tìm trị riêng của ma trận Hessenberg, tuy nhiên mình chưa tìm được cài đặt C/C++ nào của thuật toán này.

      Đoạn mã trong bài chỉ là vài dòng minh họa ý tưởng của thuật toán QR đơn giản, nó không hoạt động được trong thực tế. Mình cũng không hề nói là nó nhanh hơn hàm eig() của matlab

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s