基于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;

结果:
在这里插入图片描述2、方阵 m = n

    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。

### 使用Givens变换实现矩阵QR分解 #### Givens变换简介 Givens变换是一种用于将向量中的特定元素置零的方法。通过一系列这样的旋转操作,可以使原始矩阵逐渐变为上三角形式。每次旋转变换仅影响两个选定行之间的元素[^2]。 #### QR分解过程概述 给定一个 \( m \times n \) 的实数矩阵 A (假设 \(m \geq n\)),可以通过构建一组平面旋转矩阵 \( G_i \),使得: \[ Q = G_k ... G_2 G_1, R = G_k^T...G_2^TG_1^TA \] 最终得到的结果满足 \(A=QR\) ,其中\(Q\) 是正交矩阵而 \(R\) 则为上三角矩阵。 #### 实现步骤说明 对于每一个要处理的列 j (从左至右),依次选取该列下方非零位置 i 处的元素作为目标消元对象;构造对应的 Givens 旋转矩阵并应用到整个矩阵上直到所有指定位置都被清零为止。 #### Python代码示例 下面是一个简单的Python函数来展示如何利用Givens变换执行QR分解: ```python import numpy as np def givens_rotation(a, b): """Compute the cosine and sine of a Givens rotation.""" r = np.hypot(a, b) c = a / r s = -b / r return c, s def qr_decomposition_givens(A): m, n = A.shape Q = np.eye(m) R = A.copy() for j in range(n): for i in range(j + 1, m): if R[i,j] != 0: c, s = givens_rotation(R[j,j], R[i,j]) # Apply Givens rotation to matrix R temp_row_j = c * R[j,:] + s * R[i,:] temp_row_i = -s * R[j,:] + c * R[i,:] R[j,:] = temp_row_j R[i,:] = temp_row_i # Update corresponding rows in Q temp_q_j = c * Q[:,j] + s * Q[:,i] temp_q_i = -s * Q[:,j] + c * Q[:,i] Q[:,j] = temp_q_j Q[:,i] = temp_q_i return Q.T, R # 测试例子 if __name__ == "__main__": A = np.array([[12., -51, 4], [6, 167, -68], [-4, 24, -41]]) Q, R = qr_decomposition_givens(A) print("Matrix Q:\n", Q.round(6)) print("\nMatrix R:\n", R.round(6)) ``` 此程序实现了基本的基于Givens变换QR分解算法,并打印出了测试矩阵 `A` 对应的正交矩阵 `Q` 和上三角矩阵 `R` 。注意这里为了简化起见忽略了对输入数据的一些边界情况检查以及性能优化措施。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值