Phân tích LU

Bách khoa toàn thư mở Wikipedia
Bước tới: menu, tìm kiếm

Trong đại số tuyến tính, phân tích LU (LU decomposition, LU factorization) là phương pháp phân tích ma trận thành tích của một ma trận tam giác dưới và một ma trận tam giác trên. Phép phân tích này thường được dùng trong giải tích số để giải hệ phương trình tuyến tính hoặc tính định thức của ma trận.

Định nghĩa[sửa | sửa mã nguồn]

Gọi A là một ma trận vuông. Phân tích LU của A là cách viết A thành tích của 2 ma trận có dạng

A = LU, \,

trong đó LU lần lượt là các ma trận tam giác dưới và tam giác trên có cùng kích thước với A. Ví dụ với ma trận 3\times 3:


        \begin{bmatrix}
           a_{11} & a_{12} & a_{13} \\
           a_{21} & a_{22} & a_{23} \\
           a_{31} & a_{32} & a_{33} \\
        \end{bmatrix} =
      \begin{bmatrix}
           l_{11} & 0 & 0 \\
           l_{21} & l_{22} & 0 \\
           l_{31} & l_{32} & l_{33} \\
        \end{bmatrix}
        \begin{bmatrix}
           u_{11} & u_{12} & u_{13} \\
           0 & u_{22} & u_{23} \\
           0 & 0 & u_{33} \\
        \end{bmatrix}.

Phép phân tích LDU là cách phân tích có dạng:

 A = LDU, \,

với D là một ma trận chéo, LU là là các ma trận tam giác đơn vị, nghĩa là tất cả các phần tử trên đường chéo của LU đều bằng một.

Phép Phân tích LUP (còn gọi là LU decomposition with partial pivoting) là cách phân tích có dạng

 PA = LU, \,

với LU vẫn tương ứng là ma trận tam giác dưới và trên, và P là một ma trận hoán vị, nghĩa là P chỉ gồm không và một và chỉ có duy nhất một phần tử 1 trên mỗi dòng và cột.

Phép LU decomposition with full pivoting (Trefethen and Bau) là cách phân tích có dạng

 PAQ = LU, \,

Trong phần này chúng ta yêu cầu A là ma trận vuông, nhưng những phép phân tích này có thể được tổng quát hóa cho ma trận bất kì. Trong trường hợp đó, LP là các ma trận vuông có số dòng bằng số dòng của A, trong khi U có kích thước giống như A. Ma trận tam giác trên khi đó được hiểu là chứa toàn giá trị 0 ở dưới đường chéo chính bắt đầu từ góc trên trái.

Sự tồn tại và tính duy nhất[sửa | sửa mã nguồn]

Một ma trận khả nghịch cấp n thỏa phép phân tích LU nếu và chỉ nếu tất cả các ma trận con chính cấp k, k = 0, 1,..., n-1 của nó đều khả nghịch. Phép phân tích là duy nhất nếu ta yêu cầu các phần tử đường chéo của L (hoặc U) đều bằng 1. Ma trận cũng có phép phân tích LDU duy nhất với cùng điều kiện này.

Ma trận suy biến vẫn có thể có phân tích LU. Thực tế, một ma trận vuông hạng k vẫn có phân tích LU nếu k ma trận con chính đầu tiên của nó khả nghịch. Tuy nhiên điều ngược lại không chắc.

Người ta đã tìm được điều kiện cần và đủ để một ma trận không nhất thiết khả nghịch trong trường bất kì có phân tích LU. Điều kiện này được biểu diễn bằng hạng của một số ma trận con nào đó. Phép khử Gauss cũng được mở rộng cho trường hợp tổng quát này Okunev & Johnson 1997.

Mọi ma trận khả nghịch A - bất kể vuông hay không - đều có phân tích LUP.

Các thuật toán[sửa | sửa mã nguồn]

Phân tích LU về cơ bản là một dạng của phép khử Gaussian. Ta chuyển ma trận A thành ma trận tam giác trên U bằng cách khử các phần tử bên dưới đường chéo chính. Thuật toán Doolittle khử từng cột một bắt đầu từ bên trái, bằng cách nhân bên trái A với các ma trận tam giác dưới. Kết quả của thuật toán này là một ma trận tam giác dưới đơn vị và một ma trận tam giác trên. Thuật toán Crout hơi khác ở chỗ tạo thành một ma trận tam giác dưới và một ma trận tam giác trên đơn vị.

Phân tích LU sử dụng các thuật toán này yêu cầu khoảng 2n3 / 3 phép tính dấu chấm động.

Thuật toán Doolittle[sửa | sửa mã nguồn]

Cho một ma trận N × N


A= (a_{n,n})

ta định nghĩa

 A^{(0)} := A

và lặp với n = 1,...,N-1 như sau.

Khử các phần tử bên dưới đường chéo chính của cột thứ n của A(n-1) bằng cách cộng vào dòng thứ i của ma trận này với dòng thứ n và nhân thêm hệ số

l_{i,n} := -\frac{a_{i,n}^{(n-1)}}{a_{n,n}^{(n-1)}}

với i = n+1,\ldots,N. Nói cách khác, ta nhân bên trái A(n-1) với ma trận tam giác dưới


L_n =
\begin{pmatrix}
     1 &        &           &         &         & 0 \\
       & \ddots &           &         &         &   \\
       &        &         1 &         &         &   \\
       &        & l_{n+1,n} &  \ddots &         &   \\
       &        &    \vdots &         &  \ddots &   \\
     0 &        &   l_{N,n} &         &         & 1 \\
\end{pmatrix}.

Đặt

 A^{(n)} := L_n A^{(n-1)}.

Sau N-1 bước, ta đã khử tất cả các phần tử bên dưới đường chéo chính, và nhận được ma trận tam giác trên A(N-1). Phép phân tích LU được xác định bằng nhận xét rằng


A = L_{1}^{-1} L_{1} A^{(0)}
= L_{1}^{-1} A^{(1)} = L_{1}^{-1} L_{2}^{-1} L_{2} A^{(1)} = 
L_{1}^{-1}L_{2}^{-1} A^{(2)} =\ldots = L_{1}^{-1} \ldots L_{N-1}^{-1} A^{(N-1)}.

Kí hiệu ma trận tam giác trên A(N-1)U, và L=L_{1}^{-1} \ldots L_{N-1}^{-1}. Vì nghịch đảo của ma trận tam giác dưới 'Ln cũng là ma trận tam giác dưới, và tích hai ma trận tam giác dưới cũng là một ma trận tam giác dưới nên L là ma trận tam giác dưới cần tìm. Hơn nữa, nhận xét rằng


L =
\begin{pmatrix}
       1 &        &            &        &            & 0 \\
-l_{2,1} & \ddots &            &        &            &   \\
         &        &          1 &        &            &   \\
  \vdots &        & -l_{n+1,n} & \ddots &            &   \\
         &        &     \vdots &        &       1    &   \\
-l_{N,1} &        & -l_{N,n}   &        & -l_{N,N-1} & 1 \\
\end{pmatrix}.

Vậy ta có A=LU.

Rõ ràng là để thuật toán hoạt động được, cần phải đảm bảo a_{n,n}^{(n-1)}\not=0 tại mỗi bước (xem công thức l_{i,n}). Nếu giả sử này không đúng ở một bước nào đó, ta có thể hoán vị dòng thứ n với một dòng khác bên dưới nó để tiếp tục thuật toán. Đây là lí do mà phép phân tích LU tổng quát tương tự với phép phân tích P^{-1}A = L U .

Thuật toán Crout và LUP[sửa | sửa mã nguồn]

Thuật toán phân tích LUP của Cormen et al. là trường hợp tổng quát của phép phân tích Crout. Phần này trình bày thuật toán phân tích LUP.

  1. Nếu A có một phần tử khác không trong dòng đầu tiên, chọn ma trận hoán vị P_1 sao cho A P_1 có phần tử khác không ở góc trên trái. Ngược lại, chọn P_1ma trận đơn vị. Đặt A_1 = A P_1.
  2. Gọi A_2 là ma trận có được từ A_1 bằng cách bỏ đi cột đầu tiên và dòng đầu tiên của nó. Phân tích A_2 = L_2 U_2 P_2 theo đệ quy. Tạo L từ L_2 bằng cách thêm một dòng 0 vào phía trên và thêm cột đầu tiên của A_1 vào bên trái.
  3. Tạo U_3 từ U_2 bằng cách thêm một dòng không vào phía trên và một cột không vào bên trái, sau đó thay phần tử ở góc trên trái (đang bằng 0) thành 1. Tạo P_3 từ P_2 theo cách tương tự và định nghĩa A_3 = A_1 / P_3 = A P_1 / P_3. Gọi P là nghịch đảo của P_1 / P_3.
  4. Lúc này, A_3 bằng L U_3, (có thể) ngoại trừ dòng đầu tiên. Nếu dòng đầu tiên của A là 0, thì A_3 =  L U_3 vì cả hai đều có dòng đầu tiên là 0, và theo đó A = L U_3 P. Ngược lại, A_3L U_3 có cùng phần tử khác 0 ở góc trên trái, và A_3 = L U_3 U_1 với ma trận vuông tam giác trên U_1 với các phần tử đường chéo bằng 1 (U_1 khử các phần tử của L U_3 và thêm các phần tử của A_3). Khi đó A = L U_3 U_1 P là phép phân tích cần tìm.

Đoạn mã giả sau thể hiện từng bước của thuật toán phân tích LUP:

// A: input - ma trận ban đầu, kích thước n x n.
// n: input - kích thước của A.
// C: output - ma trận kích thước (n x n+1), ...
//              với n cột đầu tiên chứa L và U, ...
//              cột cuối cùng thể hiện ma trận hoán vị P.

FUNCTION LUP(n, A; C)
 
    // khởi tạo C
    FOR i=1 TO n DO
        C[i;n+1] = i             // khởi tạo p
        FOR j=1 TO n DO
            C[i,j] = A[i,j]
        END
    END

    FOR j=1 TO n-1 DO           // với mỗi cột j
        Chọn phần tử khác 0 lớn nhất (phần tử neo) trong cột j.
        Gọi i là dòng chứa phần tử neo đó.

        IF không tìm được i THEN
            NGỪNG THUẬT TOÁN    // Không có lời giải duy nhất
        END

        IF i ~= j THEN
            Đảo dòng i và j
        END

        FOR i = j + 1 TO n DO   // với mỗi dòng bên dưới dòng thứ j
            t = C[i,j]/C[j,j]   // thừa số
            C[i, j] = t

            FOR k = j+1 TO n DO                // với mỗi cột bên phải cột j
                C[i,k] = C[i,k] - t*C[j,k]
            END
        END
    END
END

Độ phức tạp lí thuyết[sửa | sửa mã nguồn]

Nếu nhân hai ma trận bậc n có độ phức tạp M(n), trong đó M(n)≥na với a>2 nào đó, thì phép phân tích LU có thể được tính trong thời gian O(M(n)).[1] Nghĩa là, chẳng hạn dựa trên thuật toán Coppersmith–Winograd, ta có thể có thuật toán phân tích LU với độ phức tạp O(n2.376).

Ví dụ đơn giản[sửa | sửa mã nguồn]


        \begin{bmatrix}
           4 & 3 \\
           6 & 3 \\
        \end{bmatrix} =
      \begin{bmatrix}
           l_{11} & 0 \\
           l_{21} & l_{22} \\
        \end{bmatrix}
        \begin{bmatrix}
           u_{11} & u_{12} \\
           0 & u_{22} \\
        \end{bmatrix}.

Một cách đơn giản để tìm phân tích LU của ma trận là giải hệ phương trình tuyến tính của các phép nhân ma trận tương ứng. Biết rằng:

l_{11} \cdot u_{11} + 0 \cdot 0 = 4
l_{11} \cdot u_{12} + 0 \cdot u_{22} = 3
l_{21}\cdot u_{11} + l_{22} \cdot 0 = 6
l_{21}\cdot u_{12} + l_{22} \cdot u_{22} = 3.

Hệ này có vô số nghiệm. Trong trường hợp đó bất kì hai phần tử khác 0 nào của LU đều có thể được xem là tham số, và có thể chọn bất kì giá trị khác 0 nào. Do đó để tìm phân tích LU duy nhất, ta cần thêm một số giới hạn cho LU. Chẳng hạn ta có thể yêu cầu ma trận tam giác dưới L là ma trận đơn vị (nghĩa là các phần tử đường chéo chính của nó đều bằng 1). Khi đó hệ trở thành:

l_{21} = 1.5
u_{11} = 4
u_{12} = 3
u_{22} = -1.5.

Thay các giá trị này vào ma trận, ta được:


        \begin{bmatrix}
           4 & 3 \\
           6 & 3 \\
        \end{bmatrix} =
      \begin{bmatrix}
           1 & 0 \\
           1.5 & 1 \\
        \end{bmatrix}
        \begin{bmatrix}
           4 & 3 \\
           0 & -1.5 \\
        \end{bmatrix}.

Ứng dụng[sửa | sửa mã nguồn]

Giải hệ phương trình tuyến tính[sửa | sửa mã nguồn]

Cho phương trình

A x = L U x = b \,

ta muốn giải phương trình này với Ab cho trước. Khi đó nghiệm của phương trình được tính qua 2 bước:

  1. Đầu tiên, giải phương trình  Ly = b để tìm y
  2. Sau đó giải phương trình  Ux = y để tìm x.

Nhận xét rằng trong cả hai bước ta chỉ phải làm việc với các ma trận tam giác (trên và dưới). Các phương trình này có thể được giải đơn giản bằng các phép thế thay vì sử dụng phép khử Gauss (tuy nhiên ta vẫn cần một thuật toán tương tự như phép khử Gauss để tính phân tích LU). Do đó phân tích LU chỉ tỏ ra hiệu quả khi ta phải giải phương trình trên nhiều lần với các giá trị khác nhau của b; khi đó chỉ cần tính phân tích LU của A một lần và giải các ma trận tam giác với các giá trị khác nhau của b, thay vì phải sử dụng nhiều lần các phép khử Gauss.

Tính ma trận nghịch đảo[sửa | sửa mã nguồn]

Khi giải hệ phương trình, thường thì b được xem là vector có chiều dài bằng số dòng của A. Nếu thay vì vector b, ta có ma trận B, với B là ma trận kích thước n\times p, thì ta sẽ phải tìm ma trận X (cũng có kích thước n\times p):


A X = L U X = B.

Có thể sử dụng cùng phương pháp ở trên để giải cho mỗi cột của ma trận X. Với giả sử rằng B là ma trận đơn vị với kích thước n thì X khi đó là nghịch đảo của A.[2]

Tính định thức[sửa | sửa mã nguồn]

Các ma trận LU có thể được dùng để tính định thức của ma trận A rất hiệu quả vì det(A) = det(L) det(U) và định thức của các ma trận tam giác đơn giản là tích các phần tử trên đường chéo của nó. Đặc biệt, nếu L là ma trận tam giác đơn vị thì:

 \det(A) = \det(L) \det(U) = 1 \cdot \det(U) = \prod_{i=1}^n u_{ii}.

Tương tự với phân tích LUP. Định thức của ma trận hoán vị P là (−1)S, với S là số lượng phép hoán đổi dòng trong phép phân tích.

Xem thêm[sửa | sửa mã nguồn]

Ghi chú[sửa | sửa mã nguồn]

  1. ^ J.R. Bunch and J.E. Hopcroft, Triangular factorization and inversion by fast matrix multiplication, Mathematics of Computation, 28 (1974) 231–236.
  2. ^ Matrix Computations. 3rd Edition, 1996. p121.