Bài toán tìm trị riêng của ma trận

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

Trị riêng và vector riêng có định nghĩa khá đơn giản, lại có nhiều tính chất “đẹp”, tuy nhiên không phải vì thế mà các thuật toán tìm trị riêng, vector riêng cũng dễ hiểu như vậy. Ngược lại các phương pháp hiệu quả và stable để tính trị riêng và vector riêng thường tương đối không tự nhiên. Có rất nhiều thuật toán đã được giới thiệu, và đây cũng là một trong những chủ đề thú vị trong Đại số tuyến tính, đặc biệt là khi nhìn từ góc độ ứng dụng và cài đặt trên máy tính.

Như đã biết trong bài trước, trị riêng của ma trận A là nghiệm của phương trình đa thức đặc trưng p_A\left(z\right)=det\left(zI - A\right) = 0. Để ý rằng nếu A có kích thước m\times m thì p_A\left(z\right) cũng là đa thức bậc m. Ngược lại cho đa thức bậc m (với hệ số của bậc m bằng 1) sau:

p\left(z\right)=z^m+a_{m-1}z^{m-1}+\ldots + a_1z + a_0      (1)

thì người ta cũng chứng minh được các nghiệm của đa thức (1) là các trị riêng của ma trận kích thước m\times m (gọi là ma trận liên hợpcompanion matrix – của p)

A = \left[\begin{array}{cccccc}0&\;&\;&\;&\;&-a_0\\1&0&\;&\;&\;&-a_1\\\;&1&0&\;&\;&-a_2\\\;&\;&1&\ddots&\;&\vdots\\\;&\;&\;&\ddots&0&-a_{m-2}\\\;&\;&\;&\;&1&-a_{m-1}\end{array}\right]

Thực tế có thể xem tìm trị riêng là bài toán tổng quát của việc tìm nghiệm cho đa thức. Tuy nhiên từ đây cũng dẫn ra một khó khăn lớn, đó là ta không tìm được công thức nghiệm tường minh cho đa thức bậc lớn hơn 4. Năm 1824, Abel đã chứng minh được định lí sau:

Định lí 1. Với m\geq 5 bất kì, tồn tại đa thức p(z) bậc m với hệ số hữu tỉ có một nghiệm thực p(r) = 0 mà r không thể biểu diễn bằng cách sử dụng số hữu tỉ, phép cộng, trừ, nhân, chia và căn bậc k.

Định lí này chỉ ra rằng ngay cả khi máy tính có thể làm việc với độ chính xác tuyệt đối, thì cũng không thể có thuật toán nào cho ra nghiệm chính xác của đa thức bất kì sau một số bước hữu hạn (trong khi một trong những tính chất căn bản của thuật toán là phải dừng được). Kết luận này cũng có thể áp dụng cho bài toán tổng quát hơn là tìm trị riêng của ma trận vuông: với ma trận có số chiều lớn hơn 4 thì không thể dùng đa thức đặc trưng để tìm trị riêng. Tuy vậy điều này không có nghĩa là không thể có phương pháp tìm trị riêng một cách hiệu quả. Nó chỉ ngụ ý rằng không thể có thuật toán tìm trị riêng cho ra nghiệm chính xác sau một số bước hữu hạn.

Hai trong số những bài toán lớn của Đại số tuyến tính là giải hệ phương trình tuyến tính Ax = b và tìm trị riêng của ma trận. Tuy nhiên với bài toán giải hệ phương trình tuyến tính thì các phương pháp như Householder reflector hay Gaussian elimination có thể cho ra kết quả chính xác sau một số bước hữu hạn. Ngược lại ta nhấn mạnh tính chất sau:

Mọi phương pháp tìm trị riêng đều lặp vô hạn

Mục tiêu của mọi phương pháp tìm trị riêng là tạo dãy bộ số hội tụ nhanh về trị riêng chính xác. Mặc dù yêu cầu này có vẻ dễ khiến ta nản, nhưng điều đáng mừng là các thuật toán trong lĩnh vực này đều hội tụ tương đối nhanh, và thông thường có thể dừng sau khoảng O(m) bước. Do đó mặc dù bài toán tìm trị riêng là không giải được một cách chính xác về mặt lí thuyết, nhưng trên thực tế lời giải của bài toán này có sai số rất nhỏ so với trị riêng thật sự.

1. Schur factorization

Trong bài trước, ta còn nhớ rằng mọi ma trận đều có phân tích Schur A = QTQ^*, trong đó Q là ma trận unitary và T là tam giác trên. Do T là ma trận tương đương (similarity - để tránh nhầm lẫn khi phát biểu các thuật ngữ này bằng tiếng Việt) với A nên các trị riêng của A cũng là trị riêng của T. Đồng thời do T là tam giác trên nên các trị riêng của T nằm trên đường chéo. Điều này dẫn đến hệ quả là hầu hết các phương pháp tìm trị riêng đều cố gắng thực hiện phân tích Schur bằng cách thực hiện một loạt các biến đổi X\mapsto Q_j^*XQ_j (Q_j là các ma trận unitary) sao cho tích

\underset{Q^*}{\underbrace{Q_j^*\ldots Q_2^*Q_1^*}}A\underset{Q}{\underbrace{Q_1Q_2\ldots Q_j}}         (2)

hội tụ về ma trận tam giác T khi j \rightarrow \infty.

Nếu A là ma trận thực nhưng không đối xứng thì nói chung sẽ có thể có các trị riêng phức, và dạng phân tích Schur của nó cũng sẽ là ma trận phức. Điều này rất quan trọng khi cài đặt các thuật toán tổng quát để tìm trị riêng, do người lập trình phải lường sẵn trường hợp lưu trữ và xử lí số phức.

Ngược lại nếu  A là hermitian (A = A^*) thì Q_j^*\ldots Q_2^*Q_1^*AQ_1Q_2\ldots Q_j cũng là hermitian, (2) sẽ hội tụ về ma trận vừa hermitian vừa là tam giác, do đó là ma trận đường chéo. Điều đó có nghĩa là cùng một thuật toán có thể thực hiện tam giác hóa unitary (unitarily triagularization) cho ma trận tổng quát và thực hiện chéo hóa unitary (unitarily diagonalization) cho ma trận hermitian.

2. Hai bước tìm trị riêng

Bất kể ma trận A có hermitian hay không, thông thường quá trình tính (2) được chia làm 2 bước. Bước đầu tiên biến đổi trực tiếp ma trận ban đầu về dạng Hessenberg trên (nghĩa là chứa zero ở tất cả các vị trí bên dưới đường chéo con). Bước thứ 2 là chiến lược lặp để biến ma trận Hessenberg nói trên về dạng tam giác.

\underset{A\neq A^{*}}{\left[\begin{array}{ccccc}  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\end{array}\right]}\quad\underrightarrow{\text{Phase 1}}\quad\underset{H}{\left[\begin{array}{ccccc}  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\\  \; & \times & \times & \times & \times\\  \; & \; & \times & \times & \times\\  \; & \; & \; & \times & \times\end{array}\right]}\quad\underrightarrow{\text{Phase 2}}\quad\underset{T}{\left[\begin{array}{ccccc}  \times & \times & \times & \times & \times\\  \; & \times & \times & \times & \times\\  \; & \; & \times & \times & \times\\  \; & \; & \; & \times & \times\\  \; & \; & \; & \; & \times\end{array}\right]}

Bước 1 là một phép biến đổi trực tiếp, độ phức tạp tính toán khoảng O\left(m^3\right). Bước thứ 2 trên lí thuyết thì phải lặp vô tận, nhưng trong thực tế thì thường có thể dừng sau khoảng O(m) bước, mỗi bước làm việc với ma trận Hessenberg nên có độ phức tạp khoảng O(m^2), tổng cộng độ phức tạp tính toán là O(m^3). Sơ đồ này cho thấy tầm quan trọng của bước 1. Nếu không có bước 1 thì ở bước 2 ta phải làm việc với ma trận đủ, khi đó độ phức tạp sẽ là O(m^4), hoặc hơn do khi đó có thể sẽ khó hội tụ về ma trận tam giác hơn.

Nếu  A là hermitian thì cách tiếp cận 2 bước này càng tỏ ra hiệu quả. Sau bước 1 thì ma trận trở thành hermitian Hessenberg, tức là tridiagonal. Sau bước 2 thì ma trận là hermitian triangular, do đó là ma trận đường chéo.

\underset{A=A^{*}}{\left[\begin{array}{ccccc}  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\end{array}\right]}\quad\underrightarrow{\text{Phase 1}}\quad\underset{T}{\left[\begin{array}{ccccc}  \times & \times & \; & \; & \;\\  \times & \times & \times & \; & \;\\  \; & \times & \times & \times & \;\\  \; & \; & \times & \times & \times\\  \; & \; & \; & \times & \times\end{array}\right]}\quad\underrightarrow{\text{Phase 2}}\quad\underset{D}{\left[\begin{array}{ccccc}  \times & \; & \; & \; & \;\\  \; & \times & \; & \; & \;\\  \; & \; & \times & \; & \;\\  \; & \; & \; & \times & \;\\  \; & \; & \; & \; & \times\end{array}\right]}

Trong trường hợp này, bước 2 chỉ cần khoảng O(m) phép tính cho mỗi lần lặp, do đó độ phức tạp tính toán cho cả bước 2 sẽ chỉ còn O(m^2).

3. Bước 1: Rút gọn về dạng Hessenberg

Ý tưởng ở đây là đơn giản. Ta áp dụng Householder reflector giống như khi tính QR decomposition để biến đổi ma trận ban đầu về dạng Hessenberg. Tất nhiên ta không thể áp dụng một cách ngây thơ: nhân bên trái A với Q_j^* rồi nhân bên phải với Q_j với cùng một ma trận unitary. Như thế thì toàn bộ các giá trị zero tạo ra do kết quả của phép nhân bên trái sẽ bị “phá hủy” bởi phép nhân bên phải. Rõ ràng ta cần một chiến lược thông minh hơn.

Ý tưởng như sau. Ở bước đầu tiên ta chọn Householder reflector Q_1^* sao cho nó không thay đổi dòng đầu tiên của A. Khi nhân bên trái A, nó chỉ thực hiện tổ hợp tuyến tính các dòng từ 2 đến m của A, và làm xuất hiện 0 ở các phần tử từ thứ 3 đến m trên cột đầu tiên. Sau đó khi nhân bên phải của Q_1^*A với Q_1, nó sẽ chỉ tổ hợp tuyến tính các cột từ 2 đến m, không làm thay đổi cột đầu tiên và sẽ không ảnh hưởng đến các số 0 ở cột đầu tiên vừa tạo ra.

\underset{A}{\left[\begin{array}{ccccc}  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\end{array}\right]}\quad\underrightarrow{\quad Q_{1}^{*}.\quad}\quad\underset{Q_{1}^{*}A}{\left[\begin{array}{ccccc}  \times & \times & \times & \times & \times\\  \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  0 & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  0 & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  0 & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\end{array}\right]}\quad\underrightarrow{\quad.Q_{1}\quad}\quad\underset{Q_{1}^{*}AQ_{1}}{\left[\begin{array}{ccccc}  \times & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \times & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \; & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \; & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \; & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\end{array}\right]}

Ý tưởng này được lặp lại cho các cột kế tiếp. Chẳng hạn bước thứ 2 sẽ không làm thay đổi 2 cột, 2 dòng đầu tiên:

\underset{Q_{1}^{*}AQ_{1}}{\left[\begin{array}{ccccc}  \times & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \times & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \; & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \; & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \; & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\end{array}\right]}\quad\underrightarrow{\quad Q_{2}^{*}.\quad}\quad\underset{Q_{2}^{*}Q_{1}^{*}AQ_{1}}{\left[\begin{array}{ccccc}  \times & \times & \times & \times & \times\\  \times & \times & \times & \times & \times\\  \; & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \; & \mathbf{0} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \; & \mathbf{0} & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\end{array}\right]}\quad\underrightarrow{\quad.Q_{2}\quad}\quad\underset{Q_{2}^{*}Q_{1}^{*}AQ_{1}Q_{2}}{\left[\begin{array}{ccccc}  \times & \times & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \times & \times & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \; & \times & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \; & \; & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\\  \; & \; & \boldsymbol{\times} & \boldsymbol{\times} & \boldsymbol{\times}\end{array}\right]}

Lặp lại quá trình này m-2 lần, ma trận cuối cùng sẽ là tam giác

\underset{Q^*}{\underbrace{Q_{m-2}^*\ldots Q_2^*Q_1^*}}A\underset{Q}{\underbrace{Q_1Q_2\ldots Q_{m-2}}} = H

Dựa trên thuật toán tính QR decomposition, ở đây ta có thuật toán biến đổi ma trận về dạng Hessenberg:

Thuật toán 1: Thuật toán Householder rút gọn ma trận về dạng Hessenberg

for k = 1 to m-2

x=A(k+1:m,k)

v_k = sign(x_1)\parallel x\parallel_2e_1 + x

v_k = \frac{v_k}{\parallel v_k\parallel_2}

A(k+1:m,k:m)=A(k+1:m,k:m)-2v_k\left(v_k^*A(k+1:m,k:m)\right)

A(1:m,k+1:m)=A(1:m,k+1:m)-2\left(A(1:m,k+1:m)v_k\right)v_k^*

Giống như khi tính QR decomposition, ma trận Q không được tạo một cách tường minh, tuy nhiên các vector v_k có thể được dùng để thực hiện các tính toán cần thiết, hoặc tạo lại ma trận Q khi cần.

Độ phức tạp của thuật toán này vào khoảng \frac{10}{3}m^3. Trong trường hợp ma trận ban đầu là Hermitian thì ma trận kết quả sẽ là hermitian Hessenberg, do đó là tridiagonal. Khi đó do có nhiều phần tử bằng 0 nên độ phức tạp tính toán giảm xuống còn \frac{4}{3}m^3.

4. Cài đặt MATLAB

Thuật toán khá rõ ràng do đó có thể cài đặt bằng MATLAB mà không cần tốn quá nhiều công sức. Khi nào có thời gian tôi sẽ quay lại để bổ sung cài đặt MATLAB cho thuật toán này.

Một bài sau cuối của loạt bài này sẽ nói về thuật toán QR được áp dụng trong bước 2 của phương pháp này để tính trị riêng của ma trận.

About these ads

9 Responses to “Bài toán tìm trị riêng của ma trận”

  1. ban co the giai dum minh tri rieng va vecto rieng bang lap trinh c ko

  2. Hi,
    Trong C/C++ bạn có thể sử dụng các thư viện tính toán ma trận có sẵn như BLAS (hàm syev), Eigen (http://eigen.tuxfamily.org/dox/classEigen_1_1EigenSolver.html). BLAS được cài đặt trong GotoBLAS2 (http://www.tacc.utexas.edu/tacc-projects/gotoblas2/), Intel MKL (http://software.intel.com/en-us/articles/intel-mkl/) và nhiều thư viện có sẵn khác.
    Trong trường hợp bạn cần cài đặt lại thì có thể tự cài đặt theo thuật toán mình viết trong bài. Để tính được trị riêng thì còn cần thuật toán QR nữa. Bạn có thể tìm thêm trên mạng.
    Hi vọng trong thời gian tới mình có thời gian sẽ viết tiếp về thuật toán QR và cung cấp một cài đặt đầy đủ của thuật toán này.

  3. Hi anh, em không hiểu chỗ “chọn Householder reflector sao cho nó không thay đổi dòng đầu tiên của A”, phương pháp chọn là sao vậy ạ? Xin cám ơn anh ạ.

    • Hi Thịnh,
      Trong bài QR Decomposition (http://phvuresearch.wordpress.com/2011/03/19/qr-decomposition/), ta đã dùng phương pháp Householder để biến ma trận A về dạng tam giác trên (nghĩa là các phần tử bên dưới đường chéo chính đều bằng 0).

      Trong bước 1 của thuật toán này, ta cần biến đổi ma trận A về dạng Hessenberg trên (nghĩa là các phần tử bên dưới đường chéo con – subdiagonal – đều bằng 0). Như vậy một cách tự nhiên, ta nghĩ đến việc dùng pp Householder trong thuật toán này.

      Tuy nhiên điểm khác biệt là trong bài QR Decomposition, ta chỉ chọn các Q_1, Q_2\cdots Q_k và lần lượt nhân trái với A để Q_k\cdots Q_2Q_1A là tam giác trên. Nhưng trong bài này, ta cần phải thực hiện phân tích Schur: A = QTQ^* để tìm trị riêng, nên không những nhân bên trái với Q_i^* mà ta còn phải nhân bên phải với Q_i để Q_j^*\cdots Q_2^*Q_1^*AQ_1Q_2\cdots Q_j là ma trận tam giác.

      Vấn đề là do phải nhân bên phải, nên nếu áp dụng Q_i lên toàn bộ A, thì sau khi nhân bên phải, các số 0 vừa được tạo ra sẽ lại trở thành khác 0, như vậy không thể thành ma trận tam giác được. Do đó ta chỉ áp dụng Q_i lên m-i dòng cuối của A. Cụ thể trong bước đầu tiên, i = 1, ta chỉ áp dụng Householder reflector Q_1 lên m - 1 dòng cuối của A.

      Xét ví dụ sau:
      Cho ma trận
      A = \left[\begin{array}{ccc}3 & 2 & 1 \\ 4 & 5 & 3 \\ -2 & \text{0} & 1\end{array}\right]
      Trong thuật toán QR, ta sẽ chọn Q_1^* là (Q_1^* kí hiệu ma trận chuyển vị liên hợp của Q_1, trong trường hợp này là ma trận thực nên bạn có thể coi nó là ma trận chuyển vị của Q_1):
      Q_1^* = \left[\begin{array}{ccc} -0.5571 & -0.7428 & 0.3714 \\ -0.7428 & 0.6457 & 0.1772 \\ 0.3714 & 0.1772 & 0.9114\end{array}\right]

      Q_1^*A = \left[\begin{array}{ccc}\mathbf{-5.3852} & -4.8281 & -2.4140 \\ \mathbf{-0.0000} & 1.7428 & 1.3714 \\ \mathbf{0.0000} & 1.6286 & 1.8143\end{array}\right]

      nhưng nếu nhân bên phải với Q_1 thì:
      Q_1^*AQ_1 = \left[\begin{array}{ccc}  \mathbf{5.6897} & 0.4550 & -5.0556 \\ \mathbf{-0.7852} & 1.3682 & 1.5587 \\ \mathbf{-0.5359} & 1.3730 & 1.9421\end{array}\right]

      Do vậy trong thuật toán này, ta chỉ nhân với:
      F_1^* = \left[\begin{array}{ccc} 1.0000 & \text{0} & \text{0} \\ \text{0} & -0.8944 & 0.4472 \\ \text{0} & 0.4472 & 0.8944\end{array}\right]

      F_1^*A = \left[\begin{array}{ccc} 3 & 2 & 1 \\ -4.4721 & -4.4721 & -2.2361 \\ \text{0} & 2.2361 & 2.2361\end{array}\right]

      bạn có thể thử lại trong Matlab rằng F_1^*AF_1 cũng có phần tử góc dưới trái bằng 0.
      Do F_1^* chỉ thay đổi 2 dòng dưới của A (dòng đầu tiên của F_1^*A cũng chính là dòng đầu tiên của A), nên khi nhân phải với F_1, ma trận F_1^*AF_1 vẫn có cột đầu tiên bằng cột đầu tiên của F_1^*A. Đây là ý của mình khi nói “chọn Householder reflector sao cho nó không thay đổi dòng đầu tiên của A”.

      Vấn đề còn lại là làm sao tính F_1, F_2,\cdots. Nội dung này đã được giải thích kĩ trong phần 3 của bài QR Decomposition (http://phvuresearch.wordpress.com/2011/03/19/qr-decomposition/), bạn có thể xem lại.

      Vu Pham

Trackbacks

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

Follow

Get every new post delivered to your Inbox.

Join 61 other followers

%d bloggers like this: