m和n为任意正整数,给出
A
∈
C
m
×
n
A\in C^{m\times n}
A∈Cm×n,任意矩阵都可以不需要满秩等条件,则
A
A
A可分解为
A
=
Q
R
A=QR
A=QR,其中
Q
∈
C
m
×
m
Q\in C^{m\times m}
Q∈Cm×m为一正交阵,
R
∈
C
m
×
n
R\in C^{m \times n}
R∈Cm×n为一上三角阵。
存在性:对于任何矩阵存在。
唯一性:对于实非奇异n阶矩阵或列满秩矩阵,QR分解是唯一的;对于非列满秩矩阵,QR分解可能不是唯一的。
注:
存在性
施密特正交化法:
对于实非奇异n阶矩阵A的各个列向量,通过施密特正交化可以得到n个标准正交向量,这些向量构成正交矩阵Q的列。
同时,可以得到一个上三角矩阵R,使得A=QR。
这证明了实非奇异n阶矩阵A可以进行QR分解。
Gram-Schmidt正交化过程:
对于列满秩的矩阵A,可以通过Gram-Schmidt正交化过程得到正交矩阵Q和上三角矩阵R,使得A=QR。
这个过程描述了A的列向量组进行正交化的矩阵形式,从而证明了QR分解的存在性。
Householder变换和Givens变换:
这两种变换都是正交变换,可以用来构造正交矩阵Q。
通过一系列Householder变换或Givens变换,可以将矩阵A转换为上三角矩阵R的形式,同时得到正交矩阵Q。
这证明了QR分解可以通过这两种变换来实现。
唯一性:
列满秩情况:
当矩阵A列满秩时,QR分解是唯一的。即存在一个具有正对角线元素的上三角矩阵R和一个单位列正交矩阵Q,使得A=QR。
唯一性的证明可以通过比较两个QR分解的形式,并利用正交矩阵和上三角矩阵的性质来推导。
非列满秩情况:
当矩阵A非列满秩时,虽然QR分解仍然存在,但可能不是唯一的。
这是因为非列满秩矩阵的列向量组线性相关,导致正交化过程中存在多个选择。
考虑一个位于
R
n
R^n
Rn空间的超平面,以向量
ω
\omega
ω为法向量,该超平面可表示为:
S
=
[
x
∣
ω
T
x
=
0
,
∀
x
∈
R
n
]
S=[x|\omega^Tx=0,\forall x\in R^n]
S=[x∣ωTx=0,∀x∈Rn]
该超平面是由无数垂直与
ω
\omega
ω的向量组成的,并且是
R
n
R^n
Rn的子空间,其维度/秩为n-1,因为,对于
ω
T
x
=
0
\omega^Tx=0
ωTx=0,有下列齐次线性方程组成立:
[
x
1
x
2
.
.
x
n
]
ω
=
0
\begin{bmatrix} x_1 \\x_2 \\ . \\.\\x_n\end{bmatrix}\omega=0
x1x2..xn
ω=0
对于该
A
X
=
0
AX=0
AX=0,有唯一解,因此
r
a
n
k
(
A
)
=
n
−
1
rank(A)=n-1
rank(A)=n−1。
对于任意模长为1的法向量 ω ∈ R n \omega \in{R^n} ω∈Rn,有反射矩阵 H = I − 2 ω ω T H=I-2\omega\omega^T H=I−2ωωT,显然,该反射矩阵为一对称正交阵,即满足: H = H T H=H^T H=HT, H H = I HH=I HH=I
将该反射矩阵 H H H作用与任意一个向量 Z ∈ R n Z \in R^n Z∈Rn,即 H Z HZ HZ,将得到一个与 Z Z Z关于 H H H超平面S对称的向量 Z ′ Z' Z′,如上图所示,这种反射变换即householder变换。
核心推论:
任意两个范数相同的向量
z
1
,
z
2
z_1,z_2
z1,z2,
∥
z
1
∥
2
=
∥
z
2
∥
2
,
\|z_1\|_2=\|z_2\|_2,
∥z1∥2=∥z2∥2,存在一个超平面
S
S
S,使得
z
1
z_1
z1和
z
2
z_2
z2关于超平面
S
S
S对称,因此,有householder变换矩阵
H
H
H,使得
z
1
=
H
z
2
z_1=Hz_2
z1=Hz2.
基于此推论,我们可以找到一个householder变换矩阵,使得矩阵A的最左边的向量 v v v变为仅第一行元素非0而其他元素为0的向量 [ x 0 . . 0 ] \begin{bmatrix} x \\0 \\ . \\.\\0\end{bmatrix} x0..0 , x = ∥ v ∥ 2 , x=\|v\|_2, x=∥v∥2,于是迭代此操作将矩阵A变换为一个上三角矩阵。算法思路简单,下面给出C++实现:
/**
* @brief: QR分解
* @details 将任意m.n维矩阵 A 分解为 Q(m.m) x R(m.n)
* @param {MatrixXd const&} m_in 输入矩阵
* @param {MatrixXd&} Q Q正交矩阵
* @param {MatrixXd&} R 上三角矩阵
* @return {*}
*/
void HouseholderQR(Eigen::MatrixXd const& m_in, Eigen::MatrixXd& Q, Eigen::MatrixXd& R) {
int row = m_in.rows();
int col = m_in.cols();
R = m_in;
Eigen::MatrixXd P = Eigen::MatrixXd::Identity(row, row);
for (int i = 0; i < col; i++) {
if (i == row) break;
Eigen::VectorXd v = R.block(i, i, row - i, 1);
double roi = v.norm();
Eigen::VectorXd base = Eigen::MatrixXd::Zero(row - i, 1);
base(0) = 1;
Eigen::VectorXd omega_v = v - roi * base;
omega_v.normalize();
R.block(i, i, row - i, col - i) = R.block(i, i, row - i, col - i) -
2 * omega_v * (omega_v.transpose() * R.block(i, i, row - i, col - i));
P.block(i, 0, row - i, row) = P.block(i, 0, row - i, row) -
2 * omega_v * (omega_v.transpose() * P.block(i, 0, row - i, row));
}
Q = P.transpose();
}
总的时间复杂度
O
(
n
3
)
O(n^3)
O(n3)。
注意,代码中
2 * omega_v * (omega_v.transpose() * R.block(i, i, row - i, col - i))
不要写成
2 * omega_v * omega_v.transpose() * R.block(i, i, row - i, col - i)
因为,前一种写法先执行(omega_v.transpose() * R.block(i, i, row - i, col - i)),为向量和矩阵相乘,时间复杂度
O
(
N
2
)
O(N^2)
O(N2),然后再将omega_v和括号里运算的结果进行相乘,为列向量乘行向量,时间复杂度为
O
(
N
2
)
O(N^2)
O(N2),因此叠加后总的时间复杂度依然为
O
(
N
2
)
O(N^2)
O(N2)。
而第二种写法先执行omega_v * omega_v.transpose() 变成了一个矩阵,接下来就是矩阵与矩阵相乘,时间复杂度为
O
(
N
3
)
O(N^3)
O(N3),因此我们将运算拆分成矩阵与向量的乘法,而不是矩阵与矩阵的乘法。
1、超定 m > n
Eigen::Matrix<double, 5, 4> m_in;
m_in << 3, 4, 2, 7,
7,4,9,8,
1,6,1,5,
3,4,6,2,
7,7,9,9;
结果:
Eigen::Matrix<double, 4, 4> m_in;
m_in << 3, 4, 2, 7,
7,4,9,8,
1,6,1,5,
3,4,6,2;
结果:
3、欠定 m < n
Eigen::Matrix<double, 3, 4> m_in;
m_in << 3, 4, 2, 7,
7,4,9,8,
1,6,1,5;
基于householder变换的QR分解对于任意维度的矩阵都适用,对于任意m行n列的矩阵,QR分解将其分解为
m行m列的正交矩阵Q以及m行n列的上三角矩阵R。
因篇幅问题不能全部显示,请点此查看更多更全内容