Các thuật toán phân rã ( decomposition ) cho ma trận ( matrix ) là một trong những thuật toán quan trọng nhất trong machine learning nói chung và toán ứng dụng nói riêng. Những phép phân rã này không những giúp ta giảm độ phức tạp và chi phí khi thực hiện các phép tính với các ma trận lớn mà còn có những ứng dụng đặc biệt khác.
Ở bài viết này, chúng ta sẽ cùng tìm hiểu về QR Decomposition ( Phân rã QR ) là gì và chi tiết cách tính bằng Gram-Schmidt Process ( Chu trình Gram-Schmidt ) .
Khái niệm
QR Decomposition là một phương pháp phân rã một ma trận bất kì thành tích của hai ma trận Q \mathbf{Q} Q R \mathbf{R} R Q \mathbf{Q} Q ma trận trực chuẩn ( orthonormal matrix ) và R \mathbf{R} R ma trận tam giác trên ( upper triangular matrix ) .
Lưu ý
Thuật ngữ "Orthonormal Matrix" hiếm khi được sử dụng mà thay vào đó, những ma trận như vậy sẽ được gọi là "Special Orthogonal Matrix". Tuy nhiên, trong bài viết này, chúng ta sẽ sử dụng thuật ngữ "Orthonormal Matrix" để phân biệt với "Orthogonal Matrix" một cách tường minh và dễ hiểu hơn.
Nhắc lại
Ma trận
Ma trận là một tập hợp các vector đơn vị ( unit vector ) được sắp xếp theo hàng hoặc cột. Một ma trận A ∈ R m × n \mathbf{A} \in \mathbb{R}^{m \times n} A ∈ R m × n
A = [ a 1 a 2 … a n ] = [ a 11 a 12 … a 1 n a 21 a 22 … a 2 n ⋮ ⋮ ⋱ ⋮ a m 1 a m 2 … a m n ] \begin{align}
\mathbf{A} &= \begin{bmatrix} \mathbf{a}_1 & \mathbf{a}_2 & \dots & \mathbf{a}_n \end{bmatrix} \notag \\
&= \begin{bmatrix} \mathbf{a}_{11} & \mathbf{a}_{12} & \dots & \mathbf{a}_{1n} \\ \mathbf{a}_{21} & \mathbf{a}_{22} & \dots & \mathbf{a}_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{a}_{m1} & \mathbf{a}_{m2} & \dots & \mathbf{a}_{mn} \end{bmatrix} \\
\end{align} A = [ a 1 a 2 … a n ] = a 11 a 21 ⋮ a m 1 a 12 a 22 ⋮ a m 2 … … ⋱ … a 1 n a 2 n ⋮ a mn Các vector đơn vị đóng vai trò định nghĩa lại một không gian mới bằng tổ hợp tuyến tính ( linear combination ) của chúng. Ví dụ như không gian 3 chiều cơ bản được định nghĩa bởi 3 vector đơn vị i = ( 1 , 0 , 0 ) i = (1, 0, 0) i = ( 1 , 0 , 0 ) j = ( 0 , 1 , 0 ) j = (0, 1, 0) j = ( 0 , 1 , 0 ) k = ( 0 , 0 , 1 ) k = (0, 0, 1) k = ( 0 , 0 , 1 )
R 3 = [ i j k ] = [ 1 0 0 0 1 0 0 0 1 ] \begin{align}
R^3 &= \begin{bmatrix} i & j & k \end{bmatrix} \notag \\
&= \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \\
\end{align} R 3 = [ i j k ] = 1 0 0 0 1 0 0 0 1 Ý nghĩa hình học của ma trận là nó biểu diễn một biến đổi tuyến tính ( linear transformation ) , hay nói cách khác nó giống như một phép biến đổi một không gian này sang một không gian khác mà phép biến đổi này có thể là tổ hợp của các phép quay ( rotation ) , co giãn ( stretch ) , trượt ( shear ) ,...
Ma trận A ∈ R m × n \mathbf{A} \in \mathbb{R}^{m \times n} A ∈ R m × n R n \mathbb{R}^n R n R m \mathbb{R}^m R m
Phụ thuộc Tuyến tính và Độc lập Tuyến tính
Phụ thuộc Tuyến tính
Một tập hợp các vector v 1 , v 2 , … , v n \mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_n v 1 , v 2 , … , v n phụ thuộc tuyến tính ( linearly dependent ) nếu tồn tại nghiệm là một tập hợp các hệ số c 1 , c 2 , … , c n c_1, c_2, \dots, c_n c 1 , c 2 , … , c n 0 0 0
c 1 v 1 + c 2 v 2 + ⋯ + c n v n = 0 \begin{align}
c_1 \mathbf{v}_1 + c_2 \mathbf{v}_2 + \dots + c_n \mathbf{v}_n = 0
\end{align} c 1 v 1 + c 2 v 2 + ⋯ + c n v n = 0 Điều này có nghĩa là một trong các hệ số phải khác 0 0 0 c 1 ≠ 0 c_1 \neq 0 c 1 = 0
{ v 1 = − c 2 c 1 v 2 − ⋯ − c n c 1 v n khi n > 1 v 1 = 0 khi n = 1 \begin{align}
\begin{cases}
\mathbf{v}_1 = -\frac{c_2}{c_1} \mathbf{v}_2 - \dots - \frac{c_n}{c_1} \mathbf{v}_n &\text{khi} \enspace n > 1 \\
\mathbf{v}_1 = 0 &
\text{khi} \enspace n = 1 \\
\end{cases}
\end{align} { v 1 = − c 1 c 2 v 2 − ⋯ − c 1 c n v n v 1 = 0 khi n > 1 khi n = 1 Lúc này, ta có thể thấy v 1 \mathbf{v}_1 v 1 đã không tạo ra một chiều không gian mới mà nó chỉ nằm trong không gian của các vector khác. Nói cách khác, v 1 \mathbf{v}_1 v 1 đã phụ thuộc tuyến tính vào các vector còn lại.
Độc lập Tuyến tính
Ngược lại với phụ thuộc tuyến tính , một tập hợp các vector v 1 , v 2 , … , v n \mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_n v 1 , v 2 , … , v n độc lập tuyến tính ( linearly independent ) nếu phương trình (3) chỉ có một nghiệm duy nhất là c 1 = c 2 = ⋯ = c n = 0 c_1 = c_2 = \dots = c_n = 0 c 1 = c 2 = ⋯ = c n = 0
Vì lúc này các vector đều độc lập tuyến tính với nhau và phương trình chỉ bằng 0 0 0 0 0 0
Ma trận Trực giao và Ma trận Trực chuẩn
Ma trận Trực giao
Hai vector u \mathbf{u} u v \mathbf{v} v trực giao ( orthogonal ) khi và chỉ khi chúng vuông góc với nhau, tức là u T v = 0 \mathbf{u}^T \mathbf{v} = 0 u T v = 0
Một ma trận trực giao ( orthogonal matrix ) là một tập hợp các vector mà mọi cặp vector đơn vị trong đó đều trực giao với nhau.
Gọi A \mathbf{A} A là một ma trận trực giao khi đó:
A T A = I = A A T ⟹ A − 1 = A T \begin{align}
\mathbf{A}^T \mathbf{A} = \mathbf{I} = \mathbf{A} \mathbf{A}^T \implies \mathbf{A}^{-1} = \mathbf{A}^T
\end{align} A T A = I = A A T ⟹ A − 1 = A T Ma trận Trực chuẩn
Hai vector u \mathbf{u} u v \mathbf{v} v trực chuẩn ( orthonormal ) khi và chỉ khi chúng là trực giao và đồng thời có độ dài bằng 1 1 1 u T v = 0 \mathbf{u}^T \mathbf{v} = 0 u T v = 0 ∣ ∣ u ∣ ∣ = ∣ ∣ v ∣ ∣ = 1 ||\mathbf{u}|| = ||\mathbf{v}|| = 1 ∣∣ u ∣∣ = ∣∣ v ∣∣ = 1
Một ma trận trực chuẩn là một trường hợp đặc biệt của ma trận trực giao khi nó bị ràng buộc bởi điều kiện rằng mọi vector đơn vị trong đó đều có độ dài bằng 1 1 1
Vì các tính chất trên, nếu A \mathbf{A} A là một ma trận trực chuẩn thì nó sẽ thỏa mãn toàn bộ tính chất của một ma trận trực giao và định thức ( determinant ) d e t ( A ) = 1 det(\mathbf{A}) = 1 d e t ( A ) = 1
Phép chiếu Vector
Cho trước hai vector u \mathbf{u} u v \mathbf{v} v phép chiếu vector ( vector projection ) của v \mathbf{v} v u \mathbf{u} u
p r o j u ( v ) = ⟨ u , v ⟩ ⟨ u , u ⟩ u \begin{align}
proj_{\mathbf{u}} (\mathbf{v}) &= \frac{\langle \mathbf{u}, \mathbf{v} \rangle}{\langle \mathbf{u}, \mathbf{u} \rangle} \mathbf{u} \\
\end{align} p ro j u ( v ) = ⟨ u , u ⟩ ⟨ u , v ⟩ u Trong đó:
⟨ u , v ⟩ \langle \mathbf{u}, \mathbf{v} \rangle ⟨ u , v ⟩ tích vô hướng ( inner product ) của u \mathbf{u} u v \mathbf{v} v p r o j u ( v ) proj_{\mathbf{u}} (\mathbf{v}) p ro j u ( v ) : phép chiếu vector (hình chiếu) của v \mathbf{v} v u \mathbf{u} u
Cách tính
Có hai phương pháp để thực hiện bài toán này là Gram-Schmidt Process và Householder Reflection ( Phép phản xạ Householder ) .
Chúng ta sẽ được giới thiệu về Gram-Schmidt Process trong bài viết này.
Cho trước ma trận độc lập tuyến tính A ∈ R m × n \mathbf{A} \in \mathbb{R}^{m \times n} A ∈ R m × n có các vector đơn vị a 1 , a 2 , … , a n \mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_n a 1 , a 2 , … , a n
A = [ a 1 a 2 … a n ] = [ a 11 a 12 … a 1 n a 21 a 22 … a 2 n ⋮ ⋮ ⋱ ⋮ a m 1 a m 2 … a m n ] \begin{align}
\mathbf{A} &= \begin{bmatrix} \mathbf{a}_1 & \mathbf{a}_2 & \dots & \mathbf{a}_n \end{bmatrix} \notag \\
&= \begin{bmatrix} \mathbf{a}_{11} & \mathbf{a}_{12} & \dots & \mathbf{a}_{1n} \\ \mathbf{a}_{21} & \mathbf{a}_{22} & \dots & \mathbf{a}_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{a}_{m1} & \mathbf{a}_{m2} & \dots & \mathbf{a}_{mn} \end{bmatrix} \\
\end{align} A = [ a 1 a 2 … a n ] = a 11 a 21 ⋮ a m 1 a 12 a 22 ⋮ a m 2 … … ⋱ … a 1 n a 2 n ⋮ a mn Bước 1: Trực giao hóa các unit vector
Chúng ta sẽ dùng Gram-Schmidt Process để tìm ra hệ trực giao là ma trận V ∈ R m × n \mathbf{V} \in \mathbb{R}^{m \times n} V ∈ R m × n A \mathbf{A} A . Quá trình trực giao hóa từng vector của A \mathbf{A} A
v 1 = a 1 \begin{align}
\mathbf{v}_1 &= \mathbf{a}_1 \\
\end{align} v 1 = a 1
v 2 \mathbf{v}_2 v 2 v 1 \mathbf{v}_1 v 1 a 2 \mathbf{a}_2 a 2 v 1 \mathbf{v}_1 v 1 v 2 \mathbf{v}_2 v 2 v 1 \mathbf{v}_1 v 1
v 2 = a 2 − p r o j v 1 ( a 2 ) \begin{align}
\mathbf{v}_2 &= \mathbf{a}_2 - proj_{\mathbf{v}_1} (\mathbf{a}_2) \\
\end{align} v 2 = a 2 − p ro j v 1 ( a 2 )
Tương tự với v 3 \mathbf{v}_3 v 3 , tuy nhiên để đảm bảo vector này vuông góc với tất cả các vector trước đó, nó sẽ phức tạp hơn một chút vì ta sẽ phải trực giao nó dựa trên cả v 1 \mathbf{v}_1 v 1 v 2 \mathbf{v}_2 v 2 v 1 \mathbf{v}_1 v 1 v 2 \mathbf{v}_2 v 2
v 3 = a 3 − p r o j v 1 ( a 3 ) − p r o j v 2 ( a 3 ) \begin{align}
\mathbf{v}_3 &= \mathbf{a}_3 - proj_{\mathbf{v}_1} (\mathbf{a}_3) - proj_{\mathbf{v}_2} (\mathbf{a}_3) \\
\end{align} v 3 = a 3 − p ro j v 1 ( a 3 ) − p ro j v 2 ( a 3 )
Cuối cùng, ta có công thức tổng quát cho vector thứ k k k
v k = a k − ∑ j = 1 k − 1 p r o j v j ( a k ) \begin{align}
\mathbf{v}_k = \mathbf{a}_k - \sum_{j=1}^{k-1} proj_{\mathbf{v}_j} (\mathbf{a}_k)
\end{align} v k = a k − j = 1 ∑ k − 1 p ro j v j ( a k ) Bước 2: Trực chuẩn hóa không gian vừa tính được
Ta đã có ma trận V \mathbf{V} V là một ma trận trực giao , tuy nhiên nó chưa phải là một ma trận trực chuẩn . Để trực chuẩn hóa nó, ta sẽ chia mỗi vector cột của V \mathbf{V} V 1 1 1
Lúc này, ta tìm được ma trận Q \mathbf{Q} Q là từ các vector đã được chuẩn hóa của V \mathbf{V} V
q k = v k ∣ ∣ v k ∣ ∣ \begin{align}
\mathbf{q}_k &= \frac{\mathbf{v}_k}{||\mathbf{v}_k||} \\
\end{align} q k = ∣∣ v k ∣∣ v k Bước 3: Tìm ma trận R \mathbf{R} R
Ta có thể tìm được ma trận R \mathbf{R} R bằng cách nhân ma trận Q T \mathbf{Q}^T Q T A \mathbf{A} A
A = Q R Q T A = Q T Q R Q T A = I R R = Q T A \begin{align}
\mathbf{A} &= \mathbf{Q} \mathbf{R} \notag \\
\mathbf{Q}^T \mathbf{A} &= \mathbf{Q}^T \mathbf{Q} \mathbf{R} \notag \\
\mathbf{Q}^T \mathbf{A} &= \mathbf{I} \mathbf{R} \tag*{\text{xem lại (5)}} \\
\mathbf{R} &= \mathbf{Q}^T \mathbf{A} \\
\end{align} A Q T A Q T A R = QR = Q T QR = IR = Q T A xem lại (5) Ngoài ra, ta cũng có thể tìm được ma trận R \mathbf{R} R bằng cách tính tích vô hướng của các vector đơn vị từ ma trận Q \mathbf{Q} Q với các vector đơn vị từ ma trận A \mathbf{A} A
R = [ ⟨ q 1 , a 1 ⟩ ⟨ q 1 , a 2 ⟩ … ⟨ q 1 , a n ⟩ 0 ⟨ q 2 , a 2 ⟩ … ⟨ q 2 , a n ⟩ ⋮ ⋮ ⋱ ⋮ 0 0 … ⟨ q n , a n ⟩ ] \begin{align}
\mathbf{R} = \begin{bmatrix} \langle \mathbf{q}_1 , \mathbf{a}_1 \rangle & \langle \mathbf{q}_1 , \mathbf{a}_2 \rangle & \dots & \langle \mathbf{q}_1 , \mathbf{a}_n \rangle \\ 0 & \langle \mathbf{q}_2 , \mathbf{a}_2 \rangle & \dots & \langle \mathbf{q}_2 , \mathbf{a}_n \rangle \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \langle \mathbf{q}_n , \mathbf{a}_n \rangle \end{bmatrix}
\end{align} R = ⟨ q 1 , a 1 ⟩ 0 ⋮ 0 ⟨ q 1 , a 2 ⟩ ⟨ q 2 , a 2 ⟩ ⋮ 0 … … ⋱ … ⟨ q 1 , a n ⟩ ⟨ q 2 , a n ⟩ ⋮ ⟨ q n , a n ⟩ Từ đây ta có thể hiểu được lí do R \mathbf{R} R là một ma trận tam giác trên vì trong quá trình trực giao hóa, chúng ta cho vector v k \mathbf{v}_k v k trực giao với toàn bộ các vector i < k i < k i < k , dẫn đến các tích vô hướng của chúng bằng 0 0 0 ( 9 ) (9) ( 9 )
⟨ v k , a i ⟩ = 0 ⟺ ⟨ q k , a i ⟩ = 0 ∀ i < k \begin{align}
&\langle \mathbf{v}_k , \mathbf{a}_i \rangle = 0 \notag \\
\iff &\langle \mathbf{q}_k , \mathbf{a}_i \rangle = 0 \enspace \forall i < k
\end{align} ⟺ ⟨ v k , a i ⟩ = 0 ⟨ q k , a i ⟩ = 0 ∀ i < k Ví dụ
Hãy thực hiện QR Decomposition cho ma trận A \mathbf{A} A được tạo bởi 3 vector a 1 \mathbf{a}_1 a 1 a 2 \mathbf{a}_2 a 2 a 3 \mathbf{a}_3 a 3
a 1 = ( 1 , 3 , 2 ) a 2 = ( 2 , 2 , 1 ) a 3 = ( 3 , 1 , 2 ) \begin{align}
\mathbf{a}_1 &= (1, 3, 2) \notag \\
\mathbf{a}_2 &= (2, 2, 1) \notag \\
\mathbf{a}_3 &= (3, 1, 2) \\
\end{align} a 1 a 2 a 3 = ( 1 , 3 , 2 ) = ( 2 , 2 , 1 ) = ( 3 , 1 , 2 ) A = [ a 1 a 2 a 3 ] = [ 1 2 3 3 2 1 2 1 2 ] \begin{align}
\mathbf{A} &= \begin{bmatrix} \mathbf{a}_1 & \mathbf{a}_2 & \mathbf{a}_3 \end{bmatrix} \notag \\
&= \begin{bmatrix} 1 & 2 & 3 \\ 3 & 2 & 1 \\ 2 & 1 & 2 \end{bmatrix} \\
\end{align} A = [ a 1 a 2 a 3 ] = 1 3 2 2 2 1 3 1 2 v 1 = a 1 = ( 1 , 3 , 2 ) v 2 = a 2 − p r o j v 1 ( a 2 ) = ( 2 , 2 , 1 ) − ⟨ v 1 , a 2 ⟩ ⟨ v 1 , v 1 ⟩ v 1 = ( 2 , 2 , 1 ) − 2 × 1 + 2 × 3 + 1 × 2 1 2 + 3 2 + 2 2 ( 1 , 3 , 2 ) = ( 2 , 2 , 1 ) − 5 7 ( 1 , 3 , 2 ) = 1 7 ( 9 , − 1 , − 3 ) v 3 = a 3 − p r o j v 1 ( a 3 ) − p r o j v 2 ( a 3 ) = ( 3 , 1 , 2 ) − ⟨ v 1 , a 3 ⟩ ⟨ v 1 , v 1 ⟩ v 1 − ⟨ v 2 , a 3 ⟩ ⟨ v 2 , v 2 ⟩ v 2 = ( 3 , 1 , 2 ) − 3 × 1 + 1 × 3 + 2 × 2 1 2 + 3 2 + 2 2 ( 1 , 3 , 2 ) − 1 7 ( 3 × 9 + 1 × − 1 + 2 × − 3 ) 1 7 2 ( 9 2 + ( − 1 ) 2 + ( − 3 ) 2 ) 1 7 ( 9 , − 1 , − 3 ) = ( 3 , 1 , 2 ) − 5 7 ( 1 , 3 , 2 ) − 20 91 ( 9 , − 1 , − 3 ) = 4 13 ( 1 , − 3 , 4 ) \begin{align}
\mathbf{v}_1 &= \mathbf{a}_1 \notag \\
&= (1, 3, 2) \\
\mathbf{v}_2 &= \mathbf{a}_2 - proj_{\mathbf{v}_1} (\mathbf{a}_2) \notag \\
&= (2, 2, 1) - \frac{\langle \mathbf{v}_1, \mathbf{a}_2 \rangle}{\langle \mathbf{v}_1, \mathbf{v}_1 \rangle} \mathbf{v}_1 \notag \\
&= (2, 2, 1) - \frac{2 \times 1 + 2 \times 3 + 1 \times 2}{1^2 + 3^2 + 2^2} (1, 3, 2) \notag \\
&= (2, 2, 1) - \frac{5}{7} (1, 3, 2) \notag \\
&= \frac{1}{7} (9, -1, -3) \\
\mathbf{v}_3 &= \mathbf{a}_3 - proj_{\mathbf{v}_1} (\mathbf{a}_3) - proj_{\mathbf{v}_2} (\mathbf{a}_3) \notag \\
&= (3, 1, 2) - \frac{\langle \mathbf{v}_1, \mathbf{a}_3 \rangle}{\langle \mathbf{v}_1, \mathbf{v}_1 \rangle} \mathbf{v}_1 - \frac{\langle \mathbf{v}_2, \mathbf{a}_3 \rangle}{\langle \mathbf{v}_2, \mathbf{v}_2 \rangle} \mathbf{v}_2 \notag \\
&= (3, 1, 2) - \frac{3 \times 1 + 1 \times 3 + 2 \times 2}{1^2 + 3^2 + 2^2} (1, 3, 2) - \frac{\frac{1}{7} (3 \times 9 + 1 \times -1 + 2 \times -3)}{\frac{1}{7^2} (9^2 + (-1)^2 + (-3)^2)} \frac{1}{7} (9, -1, -3) \notag \\
&= (3, 1, 2) - \frac{5}{7} (1, 3, 2) - \frac{20}{91} (9, -1, -3) \notag \\
&= \frac{4}{13} (1, -3, 4) \\
\end{align} v 1 v 2 v 3 = a 1 = ( 1 , 3 , 2 ) = a 2 − p ro j v 1 ( a 2 ) = ( 2 , 2 , 1 ) − ⟨ v 1 , v 1 ⟩ ⟨ v 1 , a 2 ⟩ v 1 = ( 2 , 2 , 1 ) − 1 2 + 3 2 + 2 2 2 × 1 + 2 × 3 + 1 × 2 ( 1 , 3 , 2 ) = ( 2 , 2 , 1 ) − 7 5 ( 1 , 3 , 2 ) = 7 1 ( 9 , − 1 , − 3 ) = a 3 − p ro j v 1 ( a 3 ) − p ro j v 2 ( a 3 ) = ( 3 , 1 , 2 ) − ⟨ v 1 , v 1 ⟩ ⟨ v 1 , a 3 ⟩ v 1 − ⟨ v 2 , v 2 ⟩ ⟨ v 2 , a 3 ⟩ v 2 = ( 3 , 1 , 2 ) − 1 2 + 3 2 + 2 2 3 × 1 + 1 × 3 + 2 × 2 ( 1 , 3 , 2 ) − 7 2 1 ( 9 2 + ( − 1 ) 2 + ( − 3 ) 2 ) 7 1 ( 3 × 9 + 1 × − 1 + 2 × − 3 ) 7 1 ( 9 , − 1 , − 3 ) = ( 3 , 1 , 2 ) − 7 5 ( 1 , 3 , 2 ) − 91 20 ( 9 , − 1 , − 3 ) = 13 4 ( 1 , − 3 , 4 ) V = [ v 1 v 2 v 3 ] = [ 1 9 1 3 − 1 − 3 2 − 3 4 ] \begin{align}
\mathbf{V} &= \begin{bmatrix} \mathbf{v}_1 & \mathbf{v}_2 & \mathbf{v}_3 \end{bmatrix} \notag \\
&= \begin{bmatrix} 1 & 9 & 1 \\ 3 & -1 & -3 \\ 2 & -3 & 4 \end{bmatrix} \\
\end{align} V = [ v 1 v 2 v 3 ] = 1 3 2 9 − 1 − 3 1 − 3 4 q 1 = v 1 ∣ ∣ v 1 ∣ ∣ = 1 1 2 + 3 2 + 2 2 ( 1 , 3 , 2 ) = 1 14 ( 1 , 3 , 2 ) q 2 = v 2 ∣ ∣ v 2 ∣ ∣ = 1 9 2 + ( − 1 ) 2 + ( − 3 ) 2 ( 9 , − 1 , − 3 ) = 1 91 ( 9 , − 1 , − 3 ) q 3 = v 3 ∣ ∣ v 3 ∣ ∣ = 1 1 2 + ( − 3 ) 2 + 4 2 ( 1 , − 3 , 4 ) = 1 26 ( 1 , − 3 , 4 ) \begin{align}
\mathbf{q}_1 &= \frac{\mathbf{v}_1}{||\mathbf{v}_1||} \notag \\
&= \frac{1}{\sqrt{1^2 + 3^2 + 2^2}} (1, 3, 2) \notag \\
&= \frac{1}{\sqrt{14}} (1, 3, 2) \\
\mathbf{q}_2 &= \frac{\mathbf{v}_2}{||\mathbf{v}_2||} \notag \\
&= \frac{1}{\sqrt{9^2 + (-1)^2 + (-3)^2}} (9, -1, -3) \notag \\
&= \frac{1}{\sqrt{91}} (9, -1, -3) \\
\mathbf{q}_3 &= \frac{\mathbf{v}_3}{||\mathbf{v}_3||} \notag \\
&= \frac{1}{\sqrt{1^2 + (-3)^2 + 4^2}} (1, -3, 4) \notag \\
&= \frac{1}{\sqrt{26}} (1, -3, 4) \\
\end{align} q 1 q 2 q 3 = ∣∣ v 1 ∣∣ v 1 = 1 2 + 3 2 + 2 2 1 ( 1 , 3 , 2 ) = 14 1 ( 1 , 3 , 2 ) = ∣∣ v 2 ∣∣ v 2 = 9 2 + ( − 1 ) 2 + ( − 3 ) 2 1 ( 9 , − 1 , − 3 ) = 91 1 ( 9 , − 1 , − 3 ) = ∣∣ v 3 ∣∣ v 3 = 1 2 + ( − 3 ) 2 + 4 2 1 ( 1 , − 3 , 4 ) = 26 1 ( 1 , − 3 , 4 ) Q = [ q 1 q 2 q 3 ] = [ 1 14 9 91 1 26 3 14 − 1 91 − 3 26 2 14 − 3 91 4 26 ] ≈ [ 0.267 0.943 0.196 0.802 − 0.105 − 0.588 0.535 − 0.314 0.784 ] \begin{align}
\mathbf{Q} &= \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 & \mathbf{q}_3 \end{bmatrix} \notag \\
&= \begin{bmatrix} \frac{1}{\sqrt{14}} & \frac{9}{\sqrt{91}} & \frac{1}{\sqrt{26}} \\ \frac{3}{\sqrt{14}} & -\frac{1}{\sqrt{91}} & -\frac{3}{\sqrt{26}} \\ \frac{2}{\sqrt{14}} & -\frac{3}{\sqrt{91}} & \frac{4}{\sqrt{26}} \end{bmatrix} \notag \\
&\approx \begin{bmatrix} 0.267 & 0.943 & 0.196 \\ 0.802 & -0.105 & -0.588 \\ 0.535 & -0.314 & 0.784 \end{bmatrix} \\
\end{align} Q = [ q 1 q 2 q 3 ] = 14 1 14 3 14 2 91 9 − 91 1 − 91 3 26 1 − 26 3 26 4 ≈ 0.267 0.802 0.535 0.943 − 0.105 − 0.314 0.196 − 0.588 0.784
Tính tích Q T A \mathbf{Q}^T \mathbf{A} Q T A
⟨ q 1 , a 1 ⟩ = 1 14 ( 1 × 1 + 3 × 3 + 2 × 2 ) = 14 ⟨ q 1 , a 2 ⟩ = 1 14 ( 1 × 2 + 3 × 2 + 2 × 1 ) = 10 14 ⟨ q 1 , a 3 ⟩ = 1 14 ( 1 × 3 + 3 × 1 + 2 × 2 ) = 10 14 ⟨ q 2 , a 2 ⟩ = 3 91 ( 9 × 2 − 1 × 2 − 3 × 1 ) = 13 91 ⟨ q 2 , a 3 ⟩ = 3 91 ( 9 × 3 − 1 × 1 − 3 × 2 ) = 20 91 ⟨ q 3 , a 3 ⟩ = 1 26 ( 1 × 3 − 3 × 1 + 4 × 2 ) = 8 26 \begin{align}
\langle \mathbf{q}_1 , \mathbf{a}_1 \rangle &= \frac{1}{\sqrt{14}} (1 \times 1 + 3 \times 3 + 2 \times 2) = \sqrt{14} \\
\langle \mathbf{q}_1 , \mathbf{a}_2 \rangle &= \frac{1}{\sqrt{14}} (1 \times 2 + 3 \times 2 + 2 \times 1) = \frac{10}{\sqrt{14}} \\
\langle \mathbf{q}_1 , \mathbf{a}_3 \rangle &= \frac{1}{\sqrt{14}} (1 \times 3 + 3 \times 1 + 2 \times 2) = \frac{10}{\sqrt{14}} \\
\langle \mathbf{q}_2 , \mathbf{a}_2 \rangle &= \frac{3}{\sqrt{91}} (9 \times 2 - 1 \times 2 - 3 \times 1) = \frac{13}{\sqrt{91}} \\
\langle \mathbf{q}_2 , \mathbf{a}_3 \rangle &= \frac{3}{\sqrt{91}} (9 \times 3 - 1 \times 1 - 3 \times 2) = \frac{20}{\sqrt{91}} \\
\langle \mathbf{q}_3 , \mathbf{a}_3 \rangle &= \frac{1}{\sqrt{26}} (1 \times 3 - 3 \times 1 + 4 \times 2) = \frac{8}{\sqrt{26}} \\
\end{align} ⟨ q 1 , a 1 ⟩ ⟨ q 1 , a 2 ⟩ ⟨ q 1 , a 3 ⟩ ⟨ q 2 , a 2 ⟩ ⟨ q 2 , a 3 ⟩ ⟨ q 3 , a 3 ⟩ = 14 1 ( 1 × 1 + 3 × 3 + 2 × 2 ) = 14 = 14 1 ( 1 × 2 + 3 × 2 + 2 × 1 ) = 14 10 = 14 1 ( 1 × 3 + 3 × 1 + 2 × 2 ) = 14 10 = 91 3 ( 9 × 2 − 1 × 2 − 3 × 1 ) = 91 13 = 91 3 ( 9 × 3 − 1 × 1 − 3 × 2 ) = 91 20 = 26 1 ( 1 × 3 − 3 × 1 + 4 × 2 ) = 26 8 R = [ ⟨ q 1 , a 1 ⟩ ⟨ q 1 , a 2 ⟩ ⟨ q 1 , a 3 ⟩ 0 ⟨ q 2 , a 2 ⟩ ⟨ q 2 , a 3 ⟩ 0 0 ⟨ q 3 , a 3 ⟩ ] = [ 14 10 14 10 14 0 13 91 20 91 0 0 8 26 ] ≈ [ 3.74165739 2.67261242 2.67261242 0 1.36277029 2.09656967 0 0 1.56892908 ] \begin{align}
\mathbf{R} &= \begin{bmatrix} \langle \mathbf{q}_1 , \mathbf{a}_1 \rangle & \langle \mathbf{q}_1 , \mathbf{a}_2 \rangle & \langle \mathbf{q}_1 , \mathbf{a}_3 \rangle \\ 0 & \langle \mathbf{q}_2 , \mathbf{a}_2 \rangle & \langle \mathbf{q}_2 , \mathbf{a}_3 \rangle \\ 0 & 0 & \langle \mathbf{q}_3 , \mathbf{a}_3 \rangle \end{bmatrix} \notag \\
&= \begin{bmatrix} \sqrt{14} & \frac{10}{\sqrt{14}} & \frac{10}{\sqrt{14}} \\ 0 & \frac{13}{\sqrt{91}} & \frac{20}{\sqrt{91}} \\ 0 & 0 & \frac{8}{\sqrt{26}} \end{bmatrix} \notag \\
&\approx \begin{bmatrix} 3.74165739 & 2.67261242 & 2.67261242 \\ 0 & 1.36277029 & 2.09656967 \\ 0 & 0 & 1.56892908 \end{bmatrix} \\
\end{align} R = ⟨ q 1 , a 1 ⟩ 0 0 ⟨ q 1 , a 2 ⟩ ⟨ q 2 , a 2 ⟩ 0 ⟨ q 1 , a 3 ⟩ ⟨ q 2 , a 3 ⟩ ⟨ q 3 , a 3 ⟩ = 14 0 0 14 10 91 13 0 14 10 91 20 26 8 ≈ 3.74165739 0 0 2.67261242 1.36277029 0 2.67261242 2.09656967 1.56892908
Như vậy, ta đã tìm được ma trận Q \mathbf{Q} Q R \mathbf{R} R
Q R = [ 1 14 9 91 1 26 3 14 − 1 91 − 3 26 2 14 − 3 91 4 26 ] [ 14 10 14 10 14 0 13 91 20 91 0 0 8 26 ] = [ 1 2 3 3 2 1 2 1 2 ] = A \begin{align}
\mathbf{Q} \mathbf{R} &= \begin{bmatrix} \frac{1}{\sqrt{14}} & \frac{9}{\sqrt{91}} & \frac{1}{\sqrt{26}} \\ \frac{3}{\sqrt{14}} & -\frac{1}{\sqrt{91}} & -\frac{3}{\sqrt{26}} \\ \frac{2}{\sqrt{14}} & -\frac{3}{\sqrt{91}} & \frac{4}{\sqrt{26}} \end{bmatrix} \begin{bmatrix} \sqrt{14} & \frac{10}{\sqrt{14}} & \frac{10}{\sqrt{14}} \\ 0 & \frac{13}{\sqrt{91}} & \frac{20}{\sqrt{91}} \\ 0 & 0 & \frac{8}{\sqrt{26}} \end{bmatrix} \notag \\
&= \begin{bmatrix} 1 & 2 & 3 \\ 3 & 2 & 1 \\ 2 & 1 & 2 \end{bmatrix} \notag \\
&= \mathbf{A} \\
\end{align} QR = 14 1 14 3 14 2 91 9 − 91 1 − 91 3 26 1 − 26 3 26 4 14 0 0 14 10 91 13 0 14 10 91 20 26 8 = 1 3 2 2 2 1 3 1 2 = A
Kiểm tra phép tính tại đây: Matrix Calculator
Triển khai code Python
Toàn bộ code có thể xem chi tiết tại: snowyfield1906/ai-general-research/least_squares/least_squares.py .
Thuật toán chính
def qr_decomposition ( A: np. array) - > ( np. array, np. array) :
m, n = A. shape
V = np. empty( ( m, 0 ) , float )
for i in range ( n) :
u_i = get( A, i)
for j in range ( i) :
v_j = get( V, j)
u_i -= orthogonalize( u_i, v_j)
v_i = u_i
V = append( V, v_i)
Q = np. empty( ( m, 0 ) , float )
for i in range ( n) :
v_i = get( V, i)
q_i = v_i / normalize( v_i)
Q = append( Q, q_i)
R = Q. T @ A
return Q, R
Các hàm phụ trợ
def inner_product ( v1: np. array, v2: np. array) - > float :
product = 0
for i in range ( len ( v1) ) :
product += v1[ i] * v2[ i]
return product
def normalize ( v: np. array) - > float :
norm = 0
for elem in v:
norm += ( elem ** 2 )
return norm ** 0.5
def orthogonalize ( u: np. array, v: np. array) - > np. array:
return ( inner_product( u, v) * v) / ( normalize( v) ** 2 )
def get ( A: np. array, i: int ) - > np. array:
a_i = A. copy( ) [ : , i]
return a_i
def append ( A: np. array, v: np. array) - > np. array:
return np. insert( A, len ( A[ 0 ] ) , v, axis= 1 )
Kiểm thử
A = np. array( [
[ 1 , 2 , 3 ] ,
[ 3 , 2 , 1 ] ,
[ 2 , 1 , 2 ] ,
] , dtype= float )
Q, R = qr_decomposition( A)
print ( "Q:" , Q)
print ( "R:" , R)
print ( "Q @ R:" , Q @ R)
> Q: [[ 0.26726124 0.94345635 0.19611614]
[ 0.80178373 -0.10482848 -0.58834841]
[ 0.53452248 -0.31448545 0.78446454]]
> R: [[3.74165739e+00 2.67261242e+00 2.67261242e+00]
[0.00000000e+00 1.36277029e+00 2.09656967e+00]
[0.00000000e+00 1.11022302e-16 1.56892908e+00]]
> Q @ R: [[1. 2. 3.]
[3. 2. 1.]
[2. 1. 2.]]
Sử dụng thư viện numpy
import numpy as np
A = np. array( [
[ 1 , 2 , 3 ] ,
[ 3 , 2 , 1 ] ,
[ 2 , 1 , 2 ] ,
] , dtype= float )
Q, R = np. linalg. qr( A)
print ( "Q:" , Q)
print ( "R:" , R)
print ( "Q @ R:" , Q @ R)
> Q: [[-0.26726124 0.94345635 0.19611614]
[-0.80178373 -0.10482848 -0.58834841]
[-0.53452248 -0.31448545 0.78446454]]
> R: [[-3.74165739 -2.67261242 -2.67261242]
[ 0. 1.36277029 2.09656967]
[ 0. 0. 1.56892908]]
> Q @ R: [[1. 2. 3.]
[3. 2. 1.]
[2. 1. 2.]]
