Bài toán Least Squares

Least squares (Bình phương tối thiểu) là phương pháp kinh điển dùng để giải xấp xỉ các hệ phương trình mà số phương trình nhiều hơn số ẩn. Ta gọi các hệ như vậy là overdetermined – dịch là quá chặt? Ứng dụng quan trọng nhất của phương pháp này là dùng để khớp dữ liệu (data fitting), giống như trong ứng dụng ban đầu của nó được Gauss đề xuất vào khoảng năm 1800. Ý tưởng chính là tìm xấp xỉ nghiệm bằng cách tối thiểu hóa độ lỗi nào đó, thông thường độ lỗi này được tính bằng một hàm bậc 2, do đó phương pháp mới có tên là bình phương tối thiểu.

Trong trường hợp phi tuyến, tức là độ lỗi không phải là hàm tuyến tính theo các ẩn, thì bài toán này có thể giải bằng phương pháp Newton-Rhapson (đã được đề cập trong một bài trước), hoặc thuật toán Levenberg-Marquardt.

Bài này chỉ xét đến trường hợp tuyến tính. Theo ngôn ngữ của Đại số tuyến tính, bài toán tìm nghiệm của một hệ overdetermined có thể xem là việc đi giải hệ phương trình tuyến tính Ax=b mà ma trận A có số dòng nhiều hơn số cột.

1. Bài toán Least squares

Xét hệ phương trình tuyến tính n ẩn nhưng có m > n phương trình. Nói cách khác, yêu cầu là ta phải tìm vector x\in \mathbb{C}^n sao cho thỏa Ax=b với A \in \mathbb{C}^{m\times n}, b\in \mathbb{C}^m.

Bài toán như vậy nói chung là vô nghiệm. Để Ax=b thì b phải là tổ hợp tuyến tính của các cột trong A, nói cách khác b phải thuộc không gian range(A) tạo bởi các cột của A. Tuy nhiên b nằm trong không gian m chiều, trong khi số chiều lớn nhất có thể của range(A) chỉ là n, nên b chỉ thuộc range(A) trong một số trường hợp đặc biệt (do ta cố ý chọn b một cách đặc biệt).

Do đó trong trường hợp tổng quát, thay vì tìm nghiệm, ta cố gắng cực tiểu hóa vector thặng dư (residual – chữ thặng dư nghe hơi ghê, trong bài này và các bài sau sẽ gọi chung là vector lỗi):

r = b - Ax \in \mathbb{C}^m            (1)

Vector này có thể được làm cho càng nhỏ càng tốt bằng cách chọn x thích hợp, tuy nhiên nói chung thì r sẽ không thể bằng zero. Để đo mức độ nhỏ của r, thông thường ta sẽ chọn một loại chuẩn nào đó. Trong trường hợp chuẩn L2 (Euclide), bài toán trở thành:

Cho A\in \mathbb{C}^{m\times n}, m\geq n, b\in \mathbb{C}^m,
tìm x\in \mathbb{C}^n sao cho cực tiểu \parallel b-Ax\parallel_2.          (2)

Đây là mô hình hóa của bài toán least squares tuyến tính. Chuẩn L2 được chọn trong tài liệu kinh điển nhờ vào những đặc tính hình học của nó: chuẩn L2 tương ứng với khoảng cách Euclide, nên ý nghĩa hình học của (2) là ta tìm x sao cho Ax là điểm trong không gian range(A) gần với b nhất. Hãy để ý nhận xét quan trọng này vì nó sẽ dẫn đến phương pháp giải bài toán Least squares, như sẽ đề cập trong phần 2.

Tuy nhiên một trong những lí do quan trọng khác là đạo hàm của hàm bậc 2 là tuyến tính. Để giải bài toán tối ưu ta sẽ đặt đạo hàm này bằng 0, do đó tuyến tính là một đặc tính quan trọng khiến chuẩn L2 được chọn. Hơn nữa dùng chuẩn L2 thì các phương trình liên quan sẽ đơn giản hơn, như trong các phần sau sẽ thấy.

2. Phương trình chuẩn

Phương trình chuẩn (normal equation) là công cụ nền tảng để giải bài toán Least squares. Tư tưởng của nó được minh họa trong Hình 1. Mục tiêu của ta là tìm một điểm y=Ax trong không gian range(A) sao cho nó gần b nhất, khi đó thì chuẩn L2 của độ lỗi r=b-Ax sẽ là nhỏ nhất. Rõ ràng điều này sẽ xảy ra khi y là hình chiếu trực giao của b lên không gian range(A), nói cách khác r=b-Ax phải trực giao với không gian range(A). Như đã biết trong bài trước, một phép chiếu trực giao lên range(A) có thể biểu diễn bằng ma trận vuông P, do đó để y=Ax là hình chiếu trực giao của b trên range(A) thì Ax = Pb. Ta nhấn mạnh bằng định lí sau:

Định lí 1: Cho ma trận A\in \mathbb{C}^{m\times n}, m\geq n và vector b\in \mathbb{C}^m. Một vector x\in\mathbb{C}^n làm tối thiểu lỗi \parallel r\parallel_2=\parallel b-Ax \parallel_2, do đó giải bài toán Least squares trong (2), nếu và chỉ nếu r \perp range(A), nghĩa là:

A^*r = 0,     (3)

và có thể viết lại thành:

A^*b=A^*Ax,      (4)

hay tương đương với:

Pb = Ax.       (5)

Trong đó P\in\mathbb{C}^{m\times m} là phép chiếu trực giao lên range(A).

Ma trận chiếu P có thể được thiết lập như đã nói ở phần 5 trong bài trước. Ta thấy (4) là hệ phương trình n\times n không suy biến nếu A hạng đủ, nói cách khác ta tìm được nghiệm x duy nhất nếu A là ma trận hạng đủ. Hệ (4) còn được gọi là phương trình chuẩn (normal equation).

3. Ma trận giả nghịch đảo

Ta vừa thấy nếu A là ma trận hạng đủ thì bài toán Least squares có nghiệm duy nhất x = \left(A^*A\right)^{-1}A^*b. Khi đó ma trận \left(A^*A\right)^{-1}A^* gọi là ma trận giả nghịch đảo (pseudoinverse) của A, kí hiệu A^+:

A^+=\left(A^*A\right)^{-1}A^*\in\mathbb{C}^{n\times m}     (6)

Ma trận này ánh xạ vector bm phần tử thành vector xn phần tử, do đó kích thước của nó là n\times m: số cột nhiều hơn số dòng.

Như vậy có thể tóm tắt bài toán Least squares trong trường hợp A là ma trận hạng đủ như sau: để giải bài toán này, ta tìm cách tính một hoặc cả 2 vector sau:

x = A^+b,    (7)

và        y=Pb         (8)

trong đó A^+ là ma trận giả nghịch đảo của A, và P là ma trận chiếu trực giao lên không gian range(A). Các phần tiếp theo lần lượt giới thiệu 3 thuật toán để giải quyết bài toán này.

4. Phương trình chuẩn

Cách cổ điển để giải bài toán Least squares là giải phương trình chuẩn (4, 7). Nếu A hạng đủ thì đó là một hệ phương trình tuyến tính n ẩn với n phương trình, hơn nữa A^+=\left(A^*A\right)^{-1}A^* là ma trận vuông, Hermitian xác định dương. Phương pháp chuẩn mực để giải hệ này là dùng phân tích Cholesky (sẽ được nhắc đến trong một bài khác). Phân tích Cholesky sẽ xây dựng phân tích A^*A=R^*R, trong đó R là ma trận tam giác trên. Phương trình chuẩn trở thành:

R^*Rx = A^*b

Thuật toán 1: Giải Least squares bằng phương trình chuẩn

  1. Xây dựng ma trận A^*A và vector A^*b.
  2. Tính phân tích Cholesky A^*A=R^*R.
  3. Giải hệ tam giác dưới R^*w=A^*b để tìm w.
  4. Giải hệ tam giác trên Rx = w để tìm x.

5. Phân tích QR

Phân tích QR cũng thường được dùng để giải bài toán Least squares. Trước tiên ta dùng phân tích QR rút gọn: A=\hat{Q}\hat{R}. Khi đó do các cột của \hat{Q} là các vector trực chuẩn và là tập sinh của không gian range(A) nên ma trận P của phép chiếu trực giao lên không gian range(A) có thể được viết thành (xem phần 4 của bài trước) P=\hat{Q}\hat{Q}^*. Vậy ta có:

y = Pb = \hat{Q}\hat{Q}^*b     (9)

Do y\in range(A) nên hệ Ax = y có lời giải xác định. Thay phân tích QR và phương trình (9) vào hệ này, ta được:

\hat{Q}\hat{R}x = \hat{Q}\hat{Q}^*b     (10)

Nhân trái 2 vế với \hat{Q}^* ta được:

\hat{R}x = \hat{Q}^*b    (11)

Phương trình này đã là một hệ tam giác trên nên có thể giải đơn giản bằng phương pháp thế.

(Để ý rằng nếu ta nhân trái 2 vế của (11) với \hat{R}^{-1} thì sẽ được ma trận giả nghịch đảo tính bằng phân tích QR: A^+=\hat{R}^{-1}\hat{Q}^*).

Thuật toán 2: Giải Least squares bằng phân tích QR

  1. Tính phân tích QR rút gọn A=\hat{Q}\hat{R}.
  2. Tính vector \hat{Q}^*b.
  3. Giải hệ tam giác trên \hat{R}x = \hat{Q}^*b để tìm x.

6. SVD

Thay vì phân tích QR, ta có thể tính phân tích SVD rút gọn của A: A=\hat{U}\hat{\Sigma}V^*, khi đó ma trận chiếu P được viết thành P=\hat{U}\hat{U}^*, suy ra:

y = Pb = \hat{U}\hat{U}^*b   (12)

Tương tự như phần trước, ta cũng có các phương trình giống như (10) và (11):

\hat{U}\hat{\Sigma}V^*x=\hat{U}\hat{U}^*b   (13)

\hat{\Sigma}V^*x=\hat{U}^*b    (14)

Nhân trái 2 vế của (14) với V\hat{\Sigma}^{-1} cũng cho ta công thức tính ma trận giả nghịch đảo theo SVD: A^+ = V\hat{\Sigma}^{-1}\hat{U}^*.

Thuật toán 3: Giải Least squares bằng phân tích SVD

  1. Tính phân tích SVD rút gọn A=\hat{U}\hat{\Sigma}V^*.
  2. Tính vector \hat{U}^*b.
  3. Giải hệ đường chéo \hat{\Sigma}w = \hat{U}^*b để tìm w.
  4. Tính x bằng phương trình x = Vw.

7. Kết luận

Ba thuật toán trên tỏ ra hữu dụng trong các trường hợp khác nhau. Khi tốc độ là yêu cầu cao thì Thuật toán 1 nên được sử dụng, do việc tính phân tích Cholesky tốn chi phí ít hơn so với QR hay SVD. Tuy nhiên thuật toán này không thực sự stable do lỗi làm tròn số. Trong trường hợp đó Thuật toán 2 nên được dùng, và có thể dùng trong hầu hết các trường hợp. Thuật toán 3 dựa trên SVD là lí tưởng cho các tình huống “xấu” của ma trận A.

Tóm lại trong bài này chúng ta đã dựa vào phép chiếu trong Đại số tuyến tính để đi tìm lời giải cho bài toán Least squares. Ba thuật toán dựa trên phân tích Cholesky, QR và SVD cũng đã được đưa ra. Cài đặt các thuật toán này rõ ràng không yêu cầu kĩ năng gì đặc biệt, do đó có thể được thực hiện dễ dàng.

Chúng ta đã nhắc tới phân tích Cholesky và SVD. Hi vọng trong các bài tới sẽ có dịp quay lại với 2 phép phân tích ma trận này.

Sau khi đọc xong bài này và bài trước về phép chiếu, ta đã có đủ kiến thức để tiếp tục tìm hiểu về thuật toán QR. Bài tiếp theo sẽ trình bày thuật toán QR tìm trị riêng cho ma trận, kết thúc chuỗi bài về trị riêng và vector riêng.

Advertisements

24 comments

  1. Em chào anh!
    E hiện đang là sinh viên cao học tại Korea University ngành Địa hóa Môi trường nước. Không hiểu lý do gì mà hồi chiều ông giáo sư của e yêu cầu em tìm hiểu về thuật toán Levenberg-Marquardt. Em đang tìm tài liệu về nó. Qua tìm hiểu e được biết nó là 1 thuật toán được ứng dụng khá nhiều. Tuy nhiên có vẻ nó không nằm trong chuyên ngành mà e đang học. Anh có thể giải đáp 1 cách ngắn gọn và dễ hiểu là cái thuật toán này là gì và ứng dụng của nó hiện nay là gì được không ạ? Em cảm ơn anh

    1. Hi Kiên,
      Mục đích của thuật toán Levenberg-Marquardt là để ước lượng tham số cho hàm số nào đó.
      Cụ thể, trong nhiều trường hợp, đôi khi ta được cho một tập giá trị \left\{x_1, x_2,\dots x_n\right\} và các giá trị \left\{y_1,y_2,\ldots y_n \right\} tương ứng. Ta cũng biết x_iy_i liên hệ với nhau bằng công thức y_i = f\left(x_i, \beta\right), trong đó f là hàm theo x_i và tham số \beta.
      Yêu cầu đặt ra là từ tập \left\{x_1, x_2,\dots x_n\right\}\left\{y_1,y_2,\ldots y_n \right\}, cho trước công thức của f\left(x_i, \beta\right), làm thế nào xác định giá trị tốt nhất của tham số \beta. “Tốt nhất”, trong trường hợp này, được hiểu là tìm \beta sao cho độ sai khác giữa y_if\left(x_i, \beta\right) là nhỏ nhất. Levenberg-Marquardt là một trong những thuật toán để làm việc này.

      Ví dụ: cho tập \left\{x_1 = 1, x_2 = 4, x_3 = 5, x_4 = 3\right\}\left\{y_1 = 100,y_2 = 140, y3 = 150, y_4 = 130 \right\}, và ta biết x_iy_i liên hệ bằng công thức y_i = \sin\left(\beta x_i\right).
      Điều này có nghĩa là ta tin rằng:
      y_1 = \sin\left(\beta x_1\right)\;\Leftrightarrow\; 100 =\sin\left(1*\beta\right)
      y_2 = \sin\left(\beta x_2\right)\;\Leftrightarrow\; 140 =\sin\left(4*\beta\right)

      y_4 = \sin\left(\beta x_4\right)\;\Leftrightarrow\; 130 =\sin\left(3*\beta\right)

      Khi đó ta có thể dùng Levenberg-Marquardt để tìm \beta sao cho nó thỏa tất cả 4 phương trình này. Trong trường hợp không thể thỏa được, thì thuật toán sẽ tìm \beta tốt nhất sao cho y_i\sin\left(\beta x_i\right) xấp xỉ gần bằng nhau nhất có thể được.

      Đó là công dụng của thuật toán Levenberg-Marquardt. Mình không rõ trong lĩnh vực của bạn nó được dùng làm gì nên không thể nói “ứng dụng” cụ thể của nó được. Tuy nhiên thuật toán này khá kinh điển và đã được cài đặt trong hầu hết các phần mềm tính toán số học.

      1. Em cảm ơn anh! Về cơ bản là e đã hiểu mục đích của thuật toán này. Hôm sau e đã trao đổi với giáo sư về mục đích của ông ấy khi yêu cầu e tìm hiểu về nó. Theo e được biết thì ông ấy có 1 bộ số liệu về chất lượng nước ngầm (kết quả phân tích) khá lớn (theo dõi trong 2 năm từ 38 giếng khoan theo các độ sâu khác nhau). Ông ấy muốn xử lý bộ số liệu này để có cái nhìn tổng quát về những biến đổi theo không gian và thời gian của chất lượng nước ngầm ở khu vực này. Một lần nữa cảm ơn anh rất nhiều

  2. Em chào anh chị!

    Em đang làm NCS, lĩnh vực Quản lý đất đai. GS hướng dẫn có yêu cầu em tìm hiểu về phương pháp Least Squares và phương pháp PCA (Principal Component Analysis)

    Mục đích là nghiên cứu các yếu tố như vị trí, hình thế, diện tích, an ninh, môi trường,… (n biến) ảnh hưởng đến giá đất và từ đó tìm ra mô hình xác định giá đất sát với giá thị trường,..

    Em đọc ở đây thấy anh chị giới thiệu về phương pháp Levenberg-Marquardt, e cũng không biết phương pháp này có thể ứng dụng cho đề tài em nghiên cứu không? mong anh chị góp ý và gợi ý giúp em!

    Trân thành cám ơn anh chị!

    1. Hi Ngọc Anh,
      Một trong những phương pháp đơn giản nhất để “dự đoán giá đất” là linear/polynomial regression. Bạn có thể tìm hiểu thêm về phương pháp này trên internet. Trong khi ta không biết dạng hàm của mô hình thì regression có vẻ phù hợp (và hay được dùng) hơn các phương pháp như kiểu Levenberg-Marquardt.

  3. Em cám ơn anh!

    Em sẽ tìm hiểu Phương pháp linear/polynomial regression. Nếu có gì không hiểu cần tham khảo ý kiến mong anh dành chút thời gian giúp đỡ em!

    Cám ơn a nhiều!

  4. em chào Thầy, em đang làm luận văn về các thuật toán tìm giá trị riêng của ma trận. Thầy có thể cho em xin tài liệu về một số thuật toán như phân tích LU, QR, SVD được không Thầy

  5. Em đang học về matrix computations, nhưng nhiều thứ em không hiểu được rõ, thầy có thể cho em xin tài liệu về môn này được không. Em xin cảm ơn. Em đang học giáo trình matrix computations của Golub.

  6. Câu đầu tiên nên là “số phương trình nhiều hơn số ẩn” mới đúng phải không ạ?

  7. Có một chỗ em thắc mắc như sau:
    Nếu gọi ma trận P là chiếu tử của phép chiếu vectơ b lên không gian range(A) thì P là ma trận vuông cấp m, b là vectơ m chiều. Khi đó hình chiếu Pb của vectơ b vào không gian range(A) vẫn là một vectơ m chiều (do phép nhân ma trận P(mxm)xb(mx1)), trong khi số chiều của range(A) là n ???
    Mong nhận hồi âm sớm từ Thầy ạ.

    1. Em đang nhầm lẫn giữa chiếu tử (projection matrix) và phép nhân ma trận thông thường (transformation matrix).

      Cho x \in \mathrm{R}^{n}P \in \mathrm{R}^{m\times n} thì y = Px \in \mathrm{R}^m, và ma trận P được gọi là linear transformation. Vector x thuộc không gian n chiều được “biến đổi” vào không gian m chiều.

      Trường hợp đặc biệt, nếu P\in \mathrm{R}^{n\times n}P^2 = P thì P gọi là chiếu tử.

      Có lẽ vì trong bài anh dùng chữ “chiếu” nên em nhầm. Nhưng anh không biết từ tiếng Việt nào phù hợp cho “transformation” cả.

  8. Vâng ạ,
    Em hiểu chiếu tử P là một ma trận chiếu (theo vật lý thì vai trò của nó như một nguồn sáng đúng không ạ).
    Theo như bài trước Anh chỉ ra thì ma trận chiếu luôn là ma trận vuông. Trong bài viết của Anh thì ma trận chiếu P phải là ma trận vuông cấp m đúng không ạ? Khi đó thì kết quả của phép chiếu P.b (thực chất là một phép nhân ma trận) với b là vector thuộc không gian m chiều vẫn là một vector trong không gian m chiều đúng không ạ? Trong khi range(A) là không gian n chiều? Đây là vấn đề em thắc mắc ạ. Mong hồi âm sớm từ Anh ạ.

  9. Vâng ạ, em hiểu sai về khái niệm không gian cột ạ. Cảm ơn Anh nhiều ạ.
    À, anh có bài viết nào chi tiết về SVD không ạ. Anh share hộ em với ạ.

  10. Anh Vũ cho em hỏi một chút ạ. Đối với phương pháp Cholesky, yêu cầu của ma trận A trong không gian phức phải là ma trận Hermitian, hạng đủ và xác định dương. Vậy làm cách nào để kiểm tra một ma trận A (kích thước lớn) là xác định dương mà không phải tính từng định thúc con góc trái không ạ. Với giả thiết A là ma trận Hermitian có dạng A = B*.B. Em cảm ơn Anh ạ.

    1. Hi em,
      Nói chung thì cho một ma trận Hermitian A \in R^{n\times n}, các mệnh đề sau là tương đương:
      – A là xác định dương
      – Tất cả trị riêng của A đều dương
      – Tất cả định thức con góc trái đều dương

      Ngoài ra còn 1 số mệnh đề nữa mà a quên rồi.

      Vậy nên anh nghĩ em có thể kiểm tra trị riêng của A.

      Nếu A có dạng đặc biệt A = B*B thì a nghĩ có thể suy ra 1 vài tính chất nào đó.

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