搜索
您的当前位置:首页正文

基于householder变换的QR分解

来源:步旅网

QR分解的定义

m和n为任意正整数,给出 A ∈ C m × n A\in C^{m\times n} ACm×n,任意矩阵都可以不需要满秩等条件,则 A A A可分解为 A = Q R A=QR A=QR,其中 Q ∈ C m × m Q\in C^{m\times m} QCm×m为一正交阵, R ∈ C m × n R\in C^{m \times n} RCm×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分解仍然存在,但可能不是唯一的。
    这是因为非列满秩矩阵的列向量组线性相关,导致正交化过程中存在多个选择。

householder变换

考虑一个位于 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,xRn]
该超平面是由无数垂直与 ω \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)=n1

对于任意模长为1的法向量 ω ∈ R n \omega \in{R^n} ωRn,有反射矩阵 H = I − 2 ω ω T H=I-2\omega\omega^T H=I2ωω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 ZRn,即 H Z HZ HZ,将得到一个与 Z Z Z关于 H H H超平面S对称的向量 Z ′ Z' Z,如上图所示,这种反射变换即householder变换。

基于householder变换的QR分解

核心推论:
任意两个范数相同的向量 z 1 , z 2 z_1,z_2 z1,z2 ∥ z 1 ∥ 2 = ∥ z 2 ∥ 2 , \|z_1\|_2=\|z_2\|_2, z12=z22,存在一个超平面 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=v2,于是迭代此操作将矩阵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。

因篇幅问题不能全部显示,请点此查看更多更全内容

Top