Phân tích ma trận: Phân tích LU

Trong bài trước, ta đã dùng phương pháp Newton để giải phương trình, hệ phương trình tổng quát. Tuy nhiên trong trường hợp hệ phương trình là tuyến tính thì có thể sử dụng một số phương pháp tốt hơn để tìm nghiệm.

Giả sử hệ phương trình cần tìm nghiệm được biểu diễn dưới dạng ma trận \mathbf{Ax} = \mathbf{b}, trong đó \mathbf{A} là ma trận m \times n\mathbf{x}, \mathbf{b} lần lượt là các vector cột có n, m phần tử.

Phương pháp kinh điển để giải hệ này là sử dụng phép khử Gaussian (thực ra còn một phương pháp “kinh điển” hơn nhiều là định thức Cramer, tuy nhiên phương pháp này không thực tế cho các hệ có kích thước lớn). Mặc dù phép khử Gaussian có thể được cài đặt thành thuật toán, nhưng chi phí sẽ khá lớn do phải thực hiện liên tục các phép biến đổi dòng trên ma trận. Một biến thể của phép khử Gaussian là LU Decomposition thường được sử dụng hơn, và sẽ được trình bày ở đây.

Trong một số bài toán thực tế như 3d reconstruction hay giải các phương trình vi phân riêng phần (PDE)…, các ma trận thường có kích thước cực lớn, trong khi số lượng các phần tử khác 0 lại không nhiều, gọi là ma trận thưa. Khi đó một số phương pháp khác như Cholesky decomposition, Gradient descent, Conjugate gradient, Biconjugate gradient, LSQR… được dùng để giải các phương trình loại này. Các phương pháp này cũng được trình bày trong phần 2.

Một lưu ý rằng mặc dù bài viết này lấy tên là “Hệ phương trình tuyến tính” – với mục tiêu là trình bày các phương pháp giải hệ phương trình tuyến tính – nhưng các phương pháp trình bày trong bài không nhất thiết chỉ dùng để giải hệ phương trình tuyến tính. Chúng có nhiều ứng dụng thú vị khác sẽ đề cập trong mỗi phần tương ứng.

1. Phân tích LU (LU Decomposition)

Cho ma trận vuông A kích thước n\times n, Phân tích LU của  A là cách viết lại A thành tích

\displaystyle 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. Phân tích LUP, LDU là một số biến thể của phân tích LU.

Đây là phương pháp cơ bản, nên thay vì viết lại, người viết đã dịch lại bài viết Phân tích LU trên wikipedia từ bản tiếng Anh, đồng thời thêm vào mã giả của thuật toán LUP. Đây là thuật toán hiện đang được sử dụng rộng rãi.

Trong matlab, việc cài đặt thuật toán LUP là rất đơn giản:


function [ L, U, P ] = mylup( A )

  [m, n] = size(A);
  if (m ~= n)
    error('mylup:input', 'Input matrix should be square.');
  end

  C = zeros(n, n+1);

  % initialize
  C(1:n, 1:n) = A;
  C(1:n, n+1) = (1:n)';

  % main algorithm
  for j=1:n-1
    if (sum(C(1:n,j)) == 0)
      error('mylup:input', 'The solution is not determined.');
    end
    [~, i] = max(abs(C(1:n, j)));

    if (i ~= j)
      tt = C(i, 1:n+1);
      C(i, 1:n+1) = C(j, 1:n+1);
      C(j, 1:n+1) = tt;
    end

    for i = j+1:n
      t = C(i, j) / C(j, j);
      C(i, j) = t;
      C(i, j+1:n) = C(i, j+1:n) - t*C(j, j+1:n);
    end
  end

  % refine the results
  PP = C(1:n, n+1);
  C = C(1:n, 1:n);
  L = tril(C, -1);
  U = triu(C);
  P = zeros(n, n);

  for i=1:n
    L(i,i)=1;
    P(i, PP(i)) = 1;
  end
end

Trong MATLAB, phân tích LU đã được cài đặt sẵn trong hàm lu() với rất nhiều tùy chọn. Việc “re-invert the wheel” sẽ giúp ta hiểu rõ hơn bản chất của thuật toán này và có thể cài đặt lại trên các môi trường khác, thay vì quá lệ thuộc vào các thư viện có sẵn.

Ứng dụng trong việc giải hệ phương trình tuyến tính Ax = LUx = b, phân tích LU được dùng thông qua 2 bước:

1. Giải phương trình Ly = b để tìm y.

2. Giải phương trình Ux = y để tìm x.

Tại mỗi bước, các ma trận đều là tam giác nên việc giải là đơn giản. Ta chỉ tốn chi phí một lần cho việc phân tích A thành LU. Rõ ràng phương pháp này chỉ tỏ ra hiệu quả hơn so với phương pháp khử Gassian cổ điển khi ta có nhu cầu giải Ax=b nhiều lần với các vector b khác nhau.

Ngoài ra có thể sử dụng phương pháp này để tính ma trận nghịch đảo, tìm định thức v.v… Các ứng dụng này đã được đề cập trên wikipedia.

2. Cholesky decomposition

3. QR decomposition

<T.B.D>

Advertisements

3 comments

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