QR分解

QR 分解的定义

将矩阵分解为一下正交阵 $Q$ 与上三角阵 $R$ 乘积的过程，称为QR分解。

且当 $R$ 的对角元均为正时，分解是唯一的。

QR 分解的方法

Schmidt 正交化

可以使用 Schmidt 正交化，根据 Schmidt 正交化过程，

$\alpha_i=\beta_i+\sum_{j=1}^{i-1} \frac{(\beta_{j},\alpha_i)}{(\beta_j,\beta_j)} \beta_j$

于是可以根据 $\beta$ 构造出 $Q$，根据前面的系数构造出 $R$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

template<class Num>
std::pair<Matrix<Num>,Matrix<Num> > QR(Matrix<Num>A){
    assert(A.row==A.col);
    Matrix<Num>R(A.row,A.col);
    auto alpha=breakAsVector(A,'C');
    int s=alpha.size();
    std::vector<Matrix<Num> >beta;
    for (int i=1;i<=s;++i){
        Matrix<Num> beta_i=alpha[i-1];
        for (int j=1;j<=i-1;++j){
            Num lambda=(alpha[i-1]&beta[j-1])/(beta[j-1]&beta[j-1]);
            beta_i=beta_i-lambda*beta[j-1];
            R[j][i]=lambda;
        }
        R[i][i]=1;
        beta_i.Message="\\beta_"+std::to_string(i);
        beta.push_back(beta_i);
    }
    for (int i=1;i<=s;++i){
        Num l=length(beta[i-1]);
        beta[i-1]=(Num(1)/Num(l))*beta[i-1];
        R.times('R',i,l);
    }
    Matrix<Num>Q=addH(beta);
    return std::make_pair(Q.message("Q"),R.message("R"));
}

HouseHolder 变换

$H=E-2\omega\omega^T$

其中 $||\omega||=1$。

$Hx$ 代表将向量 $x$ 根据以 $\omega$ 为法向量的超平面进行对称形成的新向量。

$H$ 变换是一种特殊的镜面正交变换，而 $Q$ 对应普通的正交变换，我们要把 $Q$ 表示为若干 $H$ 变换的乘积。

问题：对 $x=(x_1,\cdots,x_n)^T \not=0$，如何取 $\omega$ 和常数 $k$，使得 $Hx=ke_1$，$e_1=(1,0,\cdots,0)^T$

先确定常数 $k$，注意到 $H$ 变换不改变向量的模长，则取 $k=-sign(x_1) ||x||_2$ 即可。

这里取 $-sign(x_1)$ 是因为 $H$ 变换会改变向量的朝向，于是需要取相反的符号。

令 $u=x-ke_1$，则 $\omega=\frac{u}{||u||}$。

证明：几何方法可以证明

如图，研究 $x,e_1$ 公平面，$AB$ 代表 $x$，$AC$ 代表 $e_1$，需要找到对应法向量的方向。我们首先求点 $D$，使得三角形 $ABD$ 是等腰三角形，则 $AD=AB-AC\times |AB|$，观察到 $BE // AD$，而 $BE$ 就是我们要求的法向量，因此，得证。

注意到，可以进一步简化 $H$ 矩阵的表达式。

$H=E-2\omega\omega^T=E-\beta uu^T$

$\beta=2(||u||^2)^{-1}=2(||x-ke_1||^2)^{-1}=2(2(||x||^2+sign(x_1)x_1||x||))^{-1}\\=(||x||(||x||+|x_1|))^{-1}$

代码实现：

1
2
3
4
5
6
7
8

Matrix<Num>x=subMatrix(A,1,n,1,1);
Num lx=length(x);
Num beta=(Num)(1)/(Num)(lx*lx+lx*myabs(x[1][1]));
Num sign=(Num)(x[1][1]>=0?1:-1);
Num k=-sign*lx;
x[1][1]-=k;
Matrix<Num>H=Matrix<Num>(n,n,1)-beta*x*x.transpose();
A=H*A;

有了上述结论，如何通过 HouseHolder 变换将任意矩阵 $A$ 化为上三角矩阵？

将矩阵 $A$ 看做列向量的组合 $(x,x_1,\cdots,x_{n-1})$，我们只关心将 $x$ 化为 $ke_1$，于是，可以通过以上的变换将 $A$ 化为 $A'=(ke_1,Ax_1,\cdots,Ax_{n-1})$。

第一列已经满足上三角的形式，于是问题转化为将 $A'_{2\sim n,2\sim n}$ 的子矩阵进行 QR 分解。

1
2
3

Matrix<Num>H=Matrix<Num>(n,n,1)-beta*x*x.transpose();
A=H*A;
auto p=QR_HouseHolder(subMatrix(A,2,n,2,n));

得到了子矩阵的 $QR$ 分解后，得到 $Q',R'$，容易发现，

$Q=H^T\begin{bmatrix}1&O\\O & Q'^T\end{bmatrix}$

$R=\begin{bmatrix}x_1 &...\\O&R'\end{bmatrix}$

如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

Matrix<Num>_Q(n,n);
_Q[1][1]=1;
for (int i=1;i<=n-1;++i){
    for (int j=1;j<=n-1;++j){
        _Q[i+1][j+1]=p.first[i][j];
    }
}
Matrix<Num>R(n,n);
for (int i=1;i<=n;++i) R[1][i]=A[1][i];
for (int i=1;i<=n-1;++i){
    for (int j=1;j<=n-1;++j){
        R[i+1][j+1]=p.second[i][j];
    }
}
return std::make_pair(H.transpose()*_Q.transpose(),R);

Givens 变换

和 HouseHolder 变换差不多，这里不再赘述。

QR 分解的应用：求解特征值

幂迭代法

根据特征值的定义 $Ax=\lambda x$，这样不断迭代，向量就会趋于特征向量，当然初始值不同，趋于的特征向量是不同的。计算特征值，对于单位向量 $x$，我们可以计算 $x,Ax$ 的内积。

幂迭代法的假设：

1. $A$ 可对角化
2. $|\lambda_1|>|\lambda_2|$，当两数值非常接近的时候收敛会非常慢

反迭代算法

反迭代法和之后的带偏移的QR迭代法思路相近。

Rayleigh 商迭代

使得位移 $σ$ 尽可能地接近所求的特征值。

期望能直接给出一个理想位移是不太现实的，比较现实的方法就是动态调整，使得位移逐渐靠近某个特征值。

Rayleigh 商迭代的 $\sigma$ 选取，和带偏移的 QR 迭代法中偏移量的选取思路相近。

QR 分解法

QR 迭代的收敛性证明比较复杂，这里略去。

带位移的 QR 迭代

这里简要证明 $A_{k+1}$ 和 $A_k$ 是相似的。

$A_{k+1}=R_kQ_k+\sigma_kI=Q_k^T(A_k-\sigma_kI)Q_k+\sigma_kI=Q_k^TA_kQ_k$

$\sigma_k$ 该怎么选取呢，类似 Rayleigh 商迭代法，我们选取我们希望收敛的位置的值，将其作为 $\sigma_k$，这里，我们希望 $A_{n,n}$ 收敛到特征值，于是这里 $\sigma_k$ 选择 $A_{n,n}$，迭代结束条件怎么选择，需要特征多项式不改变，这里，我们要求 $A_{n,1\sim n-1}$ 全部为 $0$ 即可。

这里有递归的代码实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

template<class Num>
std::vector<Num> EigenVals(Matrix<Num>A){
    assert(A.row==A.col);
    int n=A.row;
    if (n==1) return std::vector<Num>{A[1][1]};
    Matrix eye=Matrix<Num>(n,n,1);
    while(true){
        Num sigma=A[n][n];
        auto p=QR(A-sigma*eye);
        Matrix Q=p.first,R=p.second;
        A=R*Q+sigma*eye;
        bool flag=true;
        for (int i=1;i<=n-1;++i){
            if (!equals(A[n][i],(Num)(0))){
                flag=false;
                break;
            }
        }
        if (flag) break;
    }
    Num eig_val=A[n][n];
    std::vector<Num> eig=EigenVals(subMatrix(A,1,n-1,1,n-1));
    eig.push_back(eig_val);
    return eig;
}

进一步的时间复杂度优化

可以使用隐式迭代的方法

1. 首先通过相似变化将 A 转化成一个 上 Hessenberg 矩阵
2. 对这个 Hessenberg 矩阵实施 隐式 QR 迭代

等有时间再写。