lanczos算法及C++实现（二）实对称阵奇异值分解的QR算法

本文属作者原创，转载请注明出处：

http://www.cnblogs.com/qxred/p/qralgorithm.html

首先推荐两个参考文献

https://www.math.kth.se/na/SF2524/matber15/qrmethod.pdf

http://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter4.pdf

1. 基本的QR算法

我们先讨论一般对阵矩阵的QR算法，再讨论对称三对角阵的QR算法

给定一个实对称阵X,假设其特征值分解为X=PSP‘，其中P对正交阵，S是对角阵。求P,S的QR算法如下，其中为正交阵，为上三角阵：


for k=1,2, ...
   （QR分解）
  
endfor

将的非对角线元素置0即得S

用matlab代码表示，即

function [P,S]=qreigen
A=rand(10);
A=A+A‘      %Symmetry real matrix, otherwise P‘AP is a block upper triangle matrix
m=size(A,1);
P=eye(m);
for i=1:1000
        [Q,R]=qr(A); %QR factorization
        A=R*Q;
        P=P*Q;
end
S=diag(diag(R));
end

这个算法很简单，每次迭代做一个QR分解和一个矩阵乘法。最基本的，我们用施密特正交化进行QR分解，这在上一篇lanczos算法及C++实现（一）框架及简单实现中使用。附录A讨论了QR方法和幂迭代方法的关系。

值得注意的是，如果矩阵X不是对称阵，这个方法往往不能收敛到一个上三角阵，而是一个块状上三角阵。（我一开始不知道，还以为程序错误，调了半天）有兴趣的同学可以试一下如下代码：

function [P,A]=qreigen
A=rand(10); %not symmetry
m=size(A,1);
P=eye(m);
for i=1:1000
        [Q,R]=qr(A); %QR factorization
        A=R*Q;
        P=P*Q;
end
% A is block upper triangle matrix, not upper triangle
A
end

运行结果如下：

A =

    4.7120   -0.6211    0.4352   -0.0450   -0.1297    0.5700   -0.3769    0.4790   -0.0754    0.2798
         0    0.0273    0.8587   -0.2188   -0.0095   -0.1913   -0.3785   -0.2142   -0.0565   -0.1236
         0   -0.9150    0.3761   -0.0914    0.0217   -0.1527    0.0317    0.4235    0.4612   -0.3284
         0   -0.0000    0.0000   -0.7400   -0.0686    0.2954    0.1561   -0.1755    0.1879   -0.2237
         0   -0.0000    0.0000   -0.0000    0.4096   -0.5238   -0.0330    0.1636    0.3090   -0.4849
         0    0.0000    0.0000   -0.0000    0.2002    0.7429    0.2064   -0.4820   -0.2602   -0.2005
         0         0         0    0.0000   -0.0000   -0.0000   -0.3812   -0.1792   -0.0188   -0.4469
         0         0         0    0.0000   -0.0000    0.0000    0.4117   -0.1513    0.1134   -0.2247
         0         0         0    0.0000   -0.0000    0.0000    0.0000    0.0000   -0.3465    0.0887
         0         0         0         0         0         0   -0.0000   -0.0000   -0.0000    0.1833

2. 快速QR算法

施密特正交法每次需要O(n^3)的计算，很慢。所以通常不用这个方法。一般来说，我们采用如下QR算法，该算法分成两步：

1、利用Householder变换 (如果你对Householder变换不了解，请见附录B)，将X正交变换为上海森伯格阵(Upper Hessenburg Matrix)，复杂度O(n^3)。

2、在每次迭代中，利用Givens变换 (如果你对Givens变换不了解，请见附录C)，求得上海森伯格阵的QR分解，复杂度O(n^2)。

这样只有在迭代之前需要O(n^3)的预处理，而每次迭代只要平方复杂度，因此比施密特正交法快很多。

2.1 矩阵的上海森伯格化 Hessenburg Reduction

在这一步，对一个方阵X(不要求对称)，我们要求一个正交变换P，使得PXP‘为上海森伯格矩阵，即低于对角线下一行的所有元素为0：

如何求得这样的P呢？基本思想如下，将矩阵X分成4块：

然后用Householder变换，将左下的一列变为只有一个非零元素：


其中是一个正交矩阵，


，

是一个除第i个分量为1以外，其他分量为0的向量。P1的具体求法，参见附录B

这样我们有

类似的，我们可以构造等共n-1个矩阵，使得X正交相似于一个上海森伯格矩阵。

值得一提的是，如果X是一个对称阵，那么根据对称性，X正交相似于一个三对角阵。

2.2 上海森伯矩阵的QR分解

首先，如果H是一个上海森伯格阵，其QR分解为H=QR，那么RQ也是一个上海森伯格阵。证明省略。

这就是说，如果我们有一个对海森伯格阵的快速QR算法，那么在每次迭代中我们可以反复使用该算法。

我们采用Givens变换来进行QR分解。对于n阶上海森伯格阵H,从上到下对每相邻的两行使用Given变换，每次变换都消去对角线下方的一个元素。共n-1次变换。由于Givens变换是正交变换，那么这n-1次变换的乘积也是正交变换，变换后，
变成了上三角阵。其中是第i次Givens行变换的矩阵表示。
根据QR算法，我们需要求RQ。因此在对H进行n-1次行变换后，再从左到右对每相邻的两列使用对应的Given变换的逆变换，同样一共也是n-1次变换：
，这样我们就完成了一次迭代。

这个过程可以用下图表示：

QR分解过程,，其中红框代表当前需要考虑的H中的元素，这些元素将被更新。注意每次Givens变换将一个对角线下的元素变成0.

而R*Q的过程如下：

其中

将这个过程中所有产生的G乘起来就得到了H的特征向量。我们可以看到，每次迭代只需要O(n^2)的计算量，比施密特正交化快多了。

2.3 三对角阵的特征值分解

在lanczos方法中，我们得到的是一个对称三对角阵，是上海森伯格阵的特例。因此可以用上面所述方法计算特征值。值得注意的是，这里可以利用三对角阵的特性加速计算，在求第一个行变换G‘1H的时候，可以只考虑H左上的2*3的小阵，而不是整个前两行。而接下来的每次行变换也只需考虑H的2*3的子阵。在列变换中，每次只需考虑H的3*2的子阵即可。这样一次QR迭代的复杂度只有O(n)了。

2.4 通过平移(shifts)和deflation加速上海森伯格矩阵的特征值分解

由于QR算法的收敛速度取决于特征值之间的差异，差别越大，QR收敛越快。这和幂迭代算法是一致的。那么怎么加快特征值之间的差异呢？通常的做法就是平移。

假设H有两个特征值，那么QR迭代的速度为
。假如满足，那么矩阵和H有相同的特征向量，其特征值为，而求H^*的特征向量会很快，其收敛速度为

基于这样的观察，在求海森伯格矩阵的特征值时，每次迭代我们进行如下QR分解

（QR分解）
  
显然，这样的QR分解不改变H的特征向量，但是如果选取合适的
收敛速度要快很多。特别地，如果恰好取到某个的特征值时，那么一次迭代即可就可确定该特征值。见文献https://www.math.kth.se/na/SF2524/matber15/qrmethod.pdf中的lemma 2.3.2。事实上我们没那么幸运，但可以用每次迭代中，H对角线上的元素可以看作是特征值的估计。所以每次计算完RQ后，取对角线上固定某个位置的元素（比如始终取第1个）作为。但如果矩阵存在一对特征值互为相反数，那么这种方法不能收敛。比如矩阵 ，因此更好的做法是Wilkinson shift，取H对角线上的2*2的方阵的特征值作为。此外对复矩阵还有双位移方法(Double shift)，这里不作详细的讨论，有兴趣的读者可以参考http://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter4.pdf

那么这样只能加速计算某个特征值，比如始终取对角线上的最后一个则加快了求解最小特征值的速度。如何加速求解其他特征值的速度呢？这里就用到了deflation。

抱歉我不知道怎么翻译这个词deflation，因此就用英文了。

Deflation是指在迭代过程中，如果对角线上某个位置的下一行几乎为0的时候，说明该位置上的值很接近于特征值了，此时将改位置所对应的行和列从矩阵中刨去，对剩下的矩阵进行计算。

一般的步骤如下，取对角线最右下的元素做deflation，直到该元素左边的元素几乎为0时，这时最小特征值已经求出，将最后一行一列去掉，对剩下的矩阵继续计算。如此反复，deflate n-1次之后，所有特征值全部求出。

2.5 QR算法的C++实现

最后，附上我的QR算法的C++实现。这个算法实现了这篇文章中的所有内容，包括实对称阵到上海森伯格阵的转化（hessenburgReduction函数），上海森伯格阵特征值分解的QR算法(QRHessenberg函数)，三对角阵的特征值分解的QR算法(QRTridiagonal函数)，每个QR算法中都用了shift和deflation加速，采用的是Wilkinson shift。在Lanczos算法中，只有QRTridiagonal函数被调用。当然如果不用lanczos方法三对角化，也可以通过hessenburgReduction函数和QRHessenberg函数来进行QR迭代求解特征值分解。具体代码参见附录 D，有三个文件svds.h, svds.cpp, main.cpp。编译时用命令g++ main.cpp svds.cpp

附录

附录 A QR算法等价于并行幂迭代算法(Simultaneous Iteration)

QR算法实际上等价于并行幂迭代算法 (Simultaneous Iteration)。在幂迭代中，如果我们用一组线性独立的向量而不是一个向量与X反复相乘，则最终得到的将是一组线性独立的向量，对其正交化后，可以 得到X的特征向量。但实际上由于幂迭代的性质，这些向量快速的收敛到主特征向量，而其他的特征向量被快速丢弃。为了保留较小特征向量的成分，每次与X相乘 后，都做一次正交化，那么这样最终这些向量将收敛到X的不同的特征向量。如果一开始选取的一组向量是单位阵，那么这种并行的幂迭代算法，实际上就是QR算 法。为什么这么说呢？我们可以看到，第一迭代，这组向量相乘X，得到X，然后正交化。而正交化实际上就是QR分解：，接着第二次迭代，我们用X乘并正交化：，如此我们得到




...

比较一下QR算法




...

用数学归纳法可以证明QR算法中的

和并行幂迭代中的

相同。因为当k=1时，显然成立，当k=2时
并行幂迭代中

而QR算法中

由于

是正交阵，一个矩阵在进行行正交变换（左乘正交阵）后，其QR分解中的R不变。因此两者的
相同
以此类推。可以证明两个算法中所有的
相同

然后，我们可以得到在并行幂迭代中，

而QR算法中

于是，可以得到两个算法中
的联系：

附录 B Householder变换

又名Householder Reflections，顾名思义，是一种镜像变换。如图，假设给定一个点x和一个过原点的平面（虚线）,平面用一个垂直于它的单位向量v表示，那么对x的Householder变换就是求x关于平面v的镜像点x‘。现在我们来看怎么计算x‘

首先0指向x的向量在向量v上的投影为

，在平面v上的投影为。所以其镜像点x‘的坐标为

这里矩阵就是一个householder矩阵，是一个对称的正交矩阵，有两个特征值，1和-1

现在我们考虑如何利用householder变换将一个向量变x成只有一个非0元素的向量,如下图所示

关键就是求那个平面，即v向量。由图可得出x-x‘平行于v，因此我们有

其中

附录 C Givens变换

又名Givens旋转，顾名思义，是对一个向量进行旋转变换，不改变向量的长度，因此是一个正交变换，其逆变换就是反方向旋转回去。

如图所示，我们考虑一个二维向量，将它逆时针旋转角度后，求所得到的坐标。

这里我们可以将向量重新表示为，其中是向量v与x坐标轴的夹角。

这样变换后的坐标为

其中为Givens变换矩阵，是一个正交阵，其逆阵为

附录 D: QR分解的C++代码

svds.h

#ifndef SVDS_H
#define SVDS_H
#include <vector>
using namespace std;
class Cell{
public:
    unsigned int row;
    unsigned int col;
    double value;
    Cell():row(0),col(0),value(0){};
    Cell(int r, int c, double v):row(r),col(c),value(v){};
};
class SparseMatrix{
public:
    unsigned int rows;
    unsigned int cols;
    vector<Cell> cells;

    int cellID;

    //read data sequentially, especially when loading large matrix from the disk.
    void moveFirst(){
        cellID=0;
    }
    bool hasNext(){
        return cellID<cells.size();
    }
    Cell next(){
        return cells[cellID++];
    }
};

void svds(SparseMatrix &A, int r, vector<vector<double> > &U, vector<double> &s, vector<vector<double> > &V);
//decompose A = U * diag(s) *V‘
//where A is a m*n matrix
//U is a m*r column orthogonal matrix
//s is a length r vector
//V is a n*r column orthogonal matrix

void print(vector<vector<double> > &X);
void hessenbergReduction(vector<vector<double> > &A, vector<vector<double> > &U);
void QRHessenbergBasic(vector<vector<double> > &A, vector<vector<double> > &Q);
void QRbasic(vector<vector<double> > &T, vector<vector<double> > &W);
void transpose(vector<vector<double> > &A, vector<vector<double> > &T);
void multiply(vector<vector<double> > &A, vector<vector<double> > &B, vector<vector<double> > &C);
void QRHessenberg(vector<vector<double> > &A, vector<vector<double> > &Q);
void QRTridiagonal(vector<vector<double> > &A, vector<vector<double> > &Q);
#endif

svds.cpp

#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <cmath>
#include <climits>
#include <cstdlib>
#include <algorithm>
#include<iomanip>
#include "svds.h"
using namespace std;
const double EPS=1e-15;
void print(vector<vector<double> > &X){
    cout.precision(6);
    cout.setf(ios::fixed);
    for(int i=0;i<X.size();i++){
        for(int j=0;j<X[i].size();j++)
            cout<<X[i][j]<<‘ ‘;
        cout<<endl;
    }
    cout<<endl;
}

void transpose(vector<vector<double> > &A, vector<vector<double> > &T){
    if(A.empty()||A[0].empty()) return;
    int m=A.size();
    int n=A[0].size();
    T.clear();
    T.resize(n,vector<double>(m,0));
    for(int i=0;i<m;i++)
        for(int j=0;j<n;j++)
            T[j][i]=A[i][j];
}
void transpose(vector<vector<double> > &A){
    int m=A.size();
    for(int i=0;i<m;i++)
        for(int j=i+1;j<m;j++)
            swap(A[i][j],A[j][i]);
}

void randUnitVector(int n, vector<double> &v){
    srand(time(NULL));
    v.clear();v.resize(n);
    while(true){
        double r=0;
        for(int i=0;i<n;i++){
            v[i]=rand()*1.0/RAND_MAX-0.5;
            r+=v[i]*v[i];
        }
        r=sqrt(r);
        if(r>EPS){
            for(int i=0;i<n;i++)
                v[i]/=r;
            break;
        }
    }
}
void multiply(vector<vector<double> > &A, vector<vector<double> > &B, vector<vector<double> > &C){
    //C=A*B
    C.clear();
    if(A.empty() || A[0].empty() || B.empty() || B[0].empty()) return ;
    C.resize(A.size(),vector<double>(B[0].size(),0));
    for(int i=0;i<A.size();i++){
        for(int j=0;j<B[0].size();j++){
            C[i][j]=0;
            for(int k=0;k<A[0].size();k++){
                C[i][j]+=A[i][k]*B[k][j];
            }
        }
    }
}

void multiply(const vector<vector<double> > &X, const vector<double> & v, vector<double> & res){
    res.clear();
    if(X.empty() || v.empty()) return;
    int m=X[0].size();
    res.resize(m,0);
    for(int i=0;i<m;i++){
        for(int j=0;X[i].size();j++){
            res[i]+=X[i][j]*v[j];
        }
    }
}
double dotProduct(const vector<double> & a, const vector<double> & b){
    double res=0;
    for(int i=0;i<a.size();i++)
        res+=a[i]*b[i];
    return res;
}

void rightMultiply(SparseMatrix &A, const vector<double> & v, vector<double> & res){
    //res= A*A‘*v
    int m=A.rows;
    int n=A.cols;
    res.clear();
    res.resize(m,0);
    vector<double> w(n,0);
    A.moveFirst();
    while(A.hasNext()){
        Cell c=A.next();
        w[c.col]+=c.value*v[c.row];
    }
    A.moveFirst();
    while(A.hasNext()){
        Cell c=A.next();
        res[c.row]+=c.value*w[c.col];
    }
}
void leftMultiply(SparseMatrix &A, const vector<double> & v, vector<double> & res){
    //res= A‘*A*v
    int m=A.rows;
    int n=A.cols;
    res.clear();
    res.resize(n,0);
    vector<double> w(m,0);
    A.moveFirst();
    while(A.hasNext()){
        Cell c=A.next();
        w[c.row]+=c.value*v[c.col];
    }
    A.moveFirst();
    while(A.hasNext()){
        Cell c=A.next();
        res[c.col]+=c.value*w[c.row];
    }
}

void rightMultiply(const vector<vector<double> > & B,SparseMatrix &A, vector<vector<double> > & C){
    //C= B‘*A
    int m=B[0].size();
    int k=B.size();
    int n=A.cols;
    for(int i=0;i<C.size();i++) fill(C[i].begin(),C[i].end(),0);
    A.moveFirst();
    while(A.hasNext()){
        Cell c=A.next();
        for(int i=0;i<m;i++){
            C[c.col][i]+=c.value*B[c.row][i];
        }
    }
}
void leftMultiply(const vector<vector<double> > & B,SparseMatrix &A, vector<vector<double> > & C){
    //C <- A B
    int r=B[0].size();
    int n=B.size();
    int m=A.rows;
    C.clear();
    C.resize(m,vector<double>(r,0));
    A.moveFirst();
    while(A.hasNext()){
        Cell c=A.next();
        for(int i=0;i<r;i++){
            C[c.row][i]+=c.value*B[c.col][i];
        }
    }
}

double norm(const vector<double> &v){
    double r=0;
    for(int i=0;i<v.size();i++)
        r+=v[i]*v[i];
    return sqrt(r);
}

double normalize(vector<double> &v){
    double r=0;
    for(int i=0;i<v.size();i++)
        r+=v[i]*v[i];
    r=sqrt(r);
    if(r>EPS){
        for(int i=0;i<v.size();i++)
            v[i]/=r;
    }
    return r;
}

void multiply(vector<double> &v, double d){
    for(int i=0;i<v.size();i++)
        v[i]*=d;
}
void hessenbergReduction(vector<vector<double> > &A, vector<vector<double> > &U){
    //A: m*m matrix
    //Reduction A to Hessenberg form: H=U‘AU (A=UHU‘), H overwrite A to save memory
    int m=A.size();
    vector<double> v(m);
    U.clear();
    U.resize(m,vector<double>(m,0));
    for(int i=0;i<m;i++)
        U[i][i]=1;
    for(int i=1;i<m;i++){
        fill(v.begin(),v.end(),0);
        for(int j=i;j<m;j++)
            v[j]=A[j][i-1];
        v[i]-=norm(v);
        normalize(v);    //P=I-2*v*v‘
        //A=PAP
        //1. PA
        for(int j=i-1;j<m;j++){
            double d=0;
            for(int k=i;k<m;k++)
                d+=A[k][j]*v[k];
            for(int k=i;k<m;k++)
                A[k][j]-=d*2*v[k];
        }
        //2. AP
        for(int j=0;j<m;j++){//row j
            double d=0;
            for(int k=i;k<m;k++)
                d+=A[j][k]*v[k];
            for(int k=i;k<m;k++)
                A[j][k]-=d*2*v[k];
        }

        //U=U*P
        for(int j=0;j<m;j++){
            double d=0;
            for(int k=i;k<m;k++)
                d+=U[j][k]*v[k];
            for(int k=i;k<m;k++)
                U[j][k]-=d*2*v[k];
        }
    }
}

void givensRotation(double a, double b, double &c, double &s){
    if(fabs(b)<EPS){
        c=1;s=0;
    }else{
        if(fabs(b)>fabs(a)){
            double r=a/b;
            s=1/sqrt(1+r*r);
            c=s*r;
        }else{
            double r=b/a;
            c=1/sqrt(1+r*r);
            s=c*r;
        }
    }
}

void QRTridiagonal(vector<vector<double> > &A, vector<vector<double> > &Q){
    //input: A m*m symmetry tridiagonal matrix
    //output: Upper triangular R, and orthogonal Q, such that A=QRQ‘
    int n=A.size();
    Q.clear();
    Q.resize(n,vector<double>(n,0));
    vector<double> cs(n-1,0);
    vector<double> ss(n-1,0);
    for(int i=0;i<n;i++) Q[i][i]=1;
    for(int m=n;m>=2;m--){
        while(1){
            fill(cs.begin(),cs.end(),0);
            fill(ss.begin(),ss.end(),0);
            double delta=(A[m-2][m-2]-A[m-1][m-1])/2;
            double sign=1;
            if(delta<0) sign=-1;
            //Wilkinson shift
            double shift=A[m-1][m-1]-sign*A[m-1][m-2]*A[m-2][m-1]/(fabs(delta)+sqrt(delta*delta+A[m-1][m-2]*A[m-2][m-1]));
            for(int i=0;i<m;i++)
                A[i][i]-=shift;
            for(int i=0;i<m-1;i++){
                double a=A[i][i];
                double b=A[i+1][i];
                givensRotation(a,b,cs[i],ss[i]);
                for(int j=i;j<=i+2 && j<m;j++){
                    a=A[i][j];
                    b=A[i+1][j];
                    A[i][j]=cs[i]*a+ss[i]*b;
                    A[i+1][j]=-ss[i]*a+cs[i]*b;
                }
            }

            for(int j=1;j<m;j++){// cols    j-1, j                           c -s
                for(int i=max(0,j-2);i<=j;i++){// rows  0 ... j        [a ,b]     s  c
                    double a=A[i][j-1];
                    double b=A[i][j];
                    A[i][j-1]=cs[j-1]*a+ss[j-1]*b;
                    A[i][j]=-ss[j-1]*a+cs[j-1]*b;
                }
            }
            for(int i=0;i<m;i++)
                A[i][i]+=shift;
            //Q=Q*G1...Gm
            for(int j=1;j<m;j++){
                for(int i=0;i<n;i++){
                    double a=Q[i][j-1];
                    double b=Q[i][j];
                    Q[i][j-1]=cs[j-1]*a+ss[j-1]*b;
                    Q[i][j]=-ss[j-1]*a+cs[j-1]*b;
                }
            }
            if(fabs(A[m-1][m-2])<1e-10)
                break;
        }
    }
}

void QRHessenberg(vector<vector<double> > &A, vector<vector<double> > &Q){
    //input: A m*m square Hessenberg Matrix
    //output: Upper triangular R, and orthogonal Q, such that A=QRQ‘
    int n=A.size();
    Q.clear();
    Q.resize(n,vector<double>(n,0));
    vector<double> cs(n-1,0);
    vector<double> ss(n-1,0);
    for(int i=0;i<n;i++) Q[i][i]=1;
    for(int m=n;m>=2;m--){
        while(1){
            fill(cs.begin(),cs.end(),0);
            fill(ss.begin(),ss.end(),0);
            double delta=(A[m-2][m-2]-A[m-1][m-1])/2;
            double sign=1;
            if(delta<0) sign=-1;
            //Wilkinson shift
            double shift=A[m-1][m-1]-sign*A[m-1][m-2]*A[m-2][m-1]/(fabs(delta)+sqrt(delta*delta+A[m-1][m-2]*A[m-2][m-1]));
            for(int i=0;i<m;i++)
                A[i][i]-=shift;
            for(int i=0;i<m-1;i++){
                double a=A[i][i];
                double b=A[i+1][i];
                givensRotation(a,b,cs[i],ss[i]);
                for(int j=i;j<m;j++){
                    a=A[i][j];
                    b=A[i+1][j];
                    A[i][j]=cs[i]*a+ss[i]*b;
                    A[i+1][j]=-ss[i]*a+cs[i]*b;
                }
            }

            for(int j=1;j<m;j++){// cols    j-1, j                           c -s
                for(int i=0;i<=j;i++){// rows  0 ... j        [a ,b]     s  c
                    double a=A[i][j-1];
                    double b=A[i][j];
                    A[i][j-1]=cs[j-1]*a+ss[j-1]*b;
                    A[i][j]=-ss[j-1]*a+cs[j-1]*b;
                }
            }
            for(int i=0;i<m;i++)
                A[i][i]+=shift;
            //Q=Q*G1...Gm
            for(int j=1;j<m;j++){
                for(int i=0;i<n;i++){
                    double a=Q[i][j-1];
                    double b=Q[i][j];
                    Q[i][j-1]=cs[j-1]*a+ss[j-1]*b;
                    Q[i][j]=-ss[j-1]*a+cs[j-1]*b;
                }
            }
            if(fabs(A[m-1][m-2])<1e-10)
                break;
        }
    }
}
void QRHessenbergBasic(vector<vector<double> > &A, vector<vector<double> > &Q){
    //input: A m*m square Hessenberg Matrix
    //output: Upper triangular R such taht A=QRQ‘ for orthogonal matrix Q
    int m=A.size();
    vector<double> cs(m-1,0);
    vector<double> ss(m-1,0);
    double diff=1;
    Q.clear();
    Q.resize(m,vector<double>(m,0));
    for(int i=0;i<m;i++) Q[i][i]=1;
    while(1){
        for(int i=0;i<m-1;i++){
            double a=A[i][i];
            double b=A[i+1][i];
            givensRotation(a,b,cs[i],ss[i]);
            for(int j=i;j<m;j++){
                a=A[i][j];
                b=A[i+1][j];
                A[i][j]=cs[i]*a+ss[i]*b;
                A[i+1][j]=-ss[i]*a+cs[i]*b;
            }
        }

        for(int j=1;j<m;j++){// cols    j-1, j                c -s
            for(int i=0;i<=j;i++){// rows  0 ... j        [a ,b]    s  c
                double a=A[i][j-1];
                double b=A[i][j];
                A[i][j-1]=cs[j-1]*a+ss[j-1]*b;
                A[i][j]=-ss[j-1]*a+cs[j-1]*b;
            }
        }
        //Q=Q*G1...Gm
        for(int j=1;j<m;j++){
            for(int i=0;i<m;i++){
                double a=Q[i][j-1];
                double b=Q[i][j];
                Q[i][j-1]=cs[j-1]*a+ss[j-1]*b;
                Q[i][j]=-ss[j-1]*a+cs[j-1]*b;
            }
        }
        diff=0;
        for(int i=0;i+1<m;i++){
            diff+=A[i+1][i]*A[i+1][i];
        }
        diff=sqrt(diff);
        if(diff<1e-10)
            break;
    }
}

void QRFactorization(vector<vector<double> > &A, vector<vector<double> > &Q, vector<vector<double> > &R){
    //A=Q*R
    //A: m*n matrix
    //Q: m*m matrix
    //R: m*n matrix
    int m=A.size();
    int n=A[0].size();
    for(int i=0;i<Q.size();i++)
        fill(Q[i].begin(),Q[i].end(),0);
    for(int i=0;i<R.size();i++)
        fill(R[i].begin(),R[i].end(),0);

    vector<double> q(m);
    vector<double> r(n);
    for(int i=0;i<n;i++){
        for(int j=0;j<m;j++)
            q[j]=A[j][i];
        fill(r.begin(),r.end(),0);
        for(int j=0;j<i && j<m;j++){
            r[j]=0;
            for(int k=0;k<q.size();k++)
                r[j]+=Q[k][j]*q[k];
            for(int k=0;k<q.size();k++)
                q[k]-=Q[k][j]*r[j];
        }
        if(i<m){
            r[i]=normalize(q);
            for(int j=0;j<m;j++)
                Q[j][i]=q[j];
        }
        for(int j=0;j<m;j++){
            R[j][i]=r[j];
        }
    }
}

void lanczos(SparseMatrix &A, vector<vector<double> > &P, vector<double> &alpha, vector<double> &beta, unsigned int rank){
    //P‘*A*A‘*P = T = diag(alpha) + diag(beta,1) + diag(beta, -1)
    //P=[p1,p2, ... , pk]
    rank=min(A.cols,min(A.rows,rank));
    vector<double> p;
    unsigned int m=A.rows;
    unsigned int n=A.cols;
    vector<double> prevP(m,0);
    randUnitVector(m,p);
    P.clear();
    P.resize(m,vector<double>(rank,0));
    vector<double> v;
    alpha.clear();alpha.resize(rank);
    beta.clear();beta.resize(rank);
    beta[0]=0;
    for(int i=0;i<rank;i++){
        for(int j=0;j<p.size();j++){
            P[j][i]=p[j];
        }
        rightMultiply(A, p, v);
        alpha[i]=dotProduct(p,v);
        if(i+1<rank){
            for(int j=0;j<m;j++)
                v[j]=v[j]-beta[i]*prevP[j]-alpha[i]*p[j];
            beta[i+1]=norm(v);
            prevP=p;
            for(int j=0;j<m;j++)
                p[j]=v[j]/beta[i+1];
        }
    }
}

void lanczosT(SparseMatrix &A, vector<vector<double> > &P, vector<double> &alpha, vector<double> &beta, unsigned int rank){
    //P‘*A‘*A*P = T = diag(alpha) + diag(beta,1) + diag(beta, -1)
    //P=[p1,p2, ... , pk]
    rank=min(A.cols,min(A.rows,rank));
    vector<double> p;
    unsigned int m=A.rows;
    unsigned int n=A.cols;
    vector<double> prevP(n,0);
    randUnitVector(n,p);
    P.clear();
    P.resize(n,vector<double>(rank,0));
    vector<double> v;
    alpha.clear();alpha.resize(rank);
    beta.clear();beta.resize(rank);
    beta[0]=0;
    for(int i=0;i<rank;i++){
        for(int j=0;j<p.size();j++){
            P[j][i]=p[j];
        }
        leftMultiply(A, p, v);
        alpha[i]=dotProduct(p,v);
        if(i+1<rank){
            for(int j=0;j<n;j++)
                v[j]=v[j]-beta[i]*prevP[j]-alpha[i]*p[j];
            beta[i+1]=norm(v);
            prevP=p;
            for(int j=0;j<n;j++)
                p[j]=v[j]/beta[i+1];
        }
    }
}

void QRbasic(vector<vector<double> > &T, vector<vector<double> > &W){
    int l=T.size();
    W.clear();
    W.resize(l,vector<double>(l,0));
    vector<vector<double> > Q(l,vector<double>(l,0));
    vector<vector<double> > R(l,vector<double>(l,0));
    vector<vector<double> > nextW(l,vector<double>(l,0));
    for(int i=0;i<l;i++)
        W[i][i]=1;
    while(true){
        //T=Q*R
        QRFactorization(T,Q,R);
        //T=R*Q
        multiply(R,Q,T);
    print(T);
        //W=W*Q;
        multiply(W,Q,nextW);
        double diff=0;
        for(int i=0;i<l;i++)
            for(int j=0;j<l;j++)
                diff+=(W[i][j]-nextW[i][j]) * (W[i][j]-nextW[i][j]);
    cout<<diff<<endl;
        W=nextW;
        if(diff<EPS*EPS)
            break;
    }
}
void svds(SparseMatrix &A, int r, vector<vector<double> > &U, vector<double> &s, vector<vector<double> > &V){
    //A=U*diag(s)*V‘
    //A:m*n matrix sparse matrix
    //U:m*r matrix, U[i]=i th left singular vector
    //s:r vector
    //V:n*r matrix, V[i]=i th right singular vector
    int m=A.rows;
    int n=A.cols;

    //lanczos: A*A‘=P*T*P‘
    if(m<=n){
        int l=m;
        vector<vector<double> > P(m,vector<double>(l,0));
        vector<double> alpha(l,0);
        vector<double> beta(l,0);
        lanczos(A,P,alpha,beta,l);
        vector<vector<double> > T(l,vector<double>(l,0));
        for(int i=0;i<l;i++){
        T[i][i]=alpha[i];
        if(i)
            T[i-1][i]=T[i][i-1]=beta[i];
        }
        vector<vector<double> > W;
        QRTridiagonal(T,W);
        U.clear();U.resize(m,vector<double>(l));
        multiply(P,W,U);
        for(int i=0;i<U.size();i++)
        U[i].resize(r);
        V.clear();V.resize(n,vector<double>(r));
        rightMultiply(U,A,V);
        s.clear();s.resize(r,0);
        for(int i=0;i<r;i++){
        for(int j=0;j<n;j++)
            s[i]+=V[j][i]*V[j][i];
        s[i]=sqrt(s[i]);
        if(s[i]>EPS){
            for(int j=0;j<n;j++)
            V[j][i]/=s[i];
        }
        }
    }else{
        int l=n;
        vector<vector<double> > P(n,vector<double>(l,0));
        vector<double> alpha(l,0);
        vector<double> beta(l,0);
        lanczosT(A,P,alpha,beta,l);
        vector<vector<double> > T(l,vector<double>(l,0));
        for(int i=0;i<l;i++){
        T[i][i]=alpha[i];
        if(i)
            T[i-1][i]=T[i][i-1]=beta[i];
        }
        vector<vector<double> > W;
        QRTridiagonal(T,W);
        V.clear();V.resize(n,vector<double>(l,0));
        U.clear();U.resize(m,vector<double>(r,0));
        multiply(P,W,V);
        for(int i=0;i<V.size();i++)
        V[i].resize(r);
        leftMultiply(V,A,U);
        s.clear();s.resize(r,0);
        for(int i=0;i<r;i++){
            for(int j=0;j<m;j++)
                s[i]+=U[j][i]*U[j][i];
            s[i]=sqrt(s[i]);
            if(s[i]>EPS){
                for(int j=0;j<m;j++)
                U[j][i]/=s[i];
            }
        }
    }
}

main.cpp

#include "svds.h"
#include <iostream>
#include <cstdlib>
int main(){

    cout.precision(6);
    cout.setf(ios::fixed);
    srand(time(NULL));
    SparseMatrix A;
    A.rows=12;
    A.cols=15;
    int rank=4;
    for(int i=0;i<A.rows;i++){
        for(int j=0;j<A.cols;j++){
            A.cells.push_back(Cell(i,j,rand()*1.0/RAND_MAX-0.5));
        }
    }

    vector<vector<double> > U,V;
    vector<double> s;
    svds(A,rank,U,s,V);

    cout<<"[U,S,V]=svds(A,r)"<<endl;
    cout<<"r = "<<rank<<endl;
    cout<<"A = "<<endl;
    A.moveFirst();
    for(int i=0;i<A.rows;i++){
        for(int j=0;j<A.cols;j++)
            cout<<A.next().value<<‘ ‘;
        cout<<endl;
    }
    cout<<endl<<endl;

    cout<<endl<<endl;

    cout<<"U = "<<endl;
    print(U);

    cout<<"s = "<<endl;
    for(int i=0;i<s.size();i++){
        for(int j=0;j<s.size();j++){
            if(i==j)
                cout<<s[i]<<‘ ‘;
            else
                cout<<0.0<<‘ ‘;
        }
        cout<<endl;
    }
    cout<<endl;
    cout<<endl;

    cout<<"V = "<<endl;
    print(V);

    return 0;
}
时间： 2025-01-04 04:27:12

lanczos算法及C++实现（二）实对称阵奇异值分解的QR算法的相关文章

实对称阵可对角化的几种证明

实对称阵是一类常见的矩阵, 它与实二次型和实内积空间上的自伴随算子有着密切的联系. 任一实对称阵 $A$ 均正交相似于对角阵, 即存在正交阵 $P$, 使得 $P'AP=\mathrm{diag\,}\{\lambda_1,\lambda_2,\cdots,\lambda_n\}$, 这是实对称阵的一条重要性质, 通常在内积空间理论的框架中加以证明. 然而, 实对称阵可对角化这一性质可以在引入矩阵可对角化的定义和判定准则后直接加以证明, 也可以利用 Jordan 标准型理论加以证明. 下面给出实

【Math】证明：实对称阵属于不同特征值的的特征向量是正交的

证明:实对称阵属于不同特征值的的特征向量是正交的. 设Ap=mp,Aq=nq,其中A是实对称矩阵,m,n为其不同的特征值,p,q分别为其对应得特征向量. 则 p1(Aq)=p1(nq)=np1q (p1A)q=(p1A1)q=(AP)1q=(mp)1q=mp1q 因为 p1(Aq)= (p1A)q 上两式作差得: (m-n)p1q=0 由于m不等于n, 所以p1q=0 即(p,q)=0,从而p,q正交.说明:p1表示p的转置,A1表示A的转置,(Ap)1表示Ap的转置

lanczos算法及C++实现（三）实对称三对角阵特征值分解的分治算法

本文属作者原创,转载请注明出处 http://www.cnblogs.com/qxred/p/dcalgorithm.html 本系列目录: lanczos算法及C++实现(一)框架及简单实现 lanczos算法及C++实现(二)实对称阵奇异值分解的QR算法 lanczos算法及C++实现(三)实对称三对角阵特征值分解的分治算法 0. 参考文献 https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm A. Mel

关于正交阵与实对称正交阵的专题讨论

正交阵 $\bf命题:$设$A,B$为$n$阶实正交阵,且$\det \left( {A + B} \right) = \det \left( A \right) -\det \left( B \right)$,证明:$\det \left( A \right) = \det \left( B \right)$ 1 $\bf命题:$设$A,B$为$n$阶正交阵,则$n-r(A+B)$为偶数当且仅当$\det \left( A \right) = \det \left( B \right)$ 1

求实对称阵的 特征值 和 特征向量(转)

/// <summary> /// 求实对称阵的 特征值 和 特征向量 /// </summary> /// <param name="data">实对称阵</param> /// <param name="num">维数</param> /// <param name="eigenvalue">引用参数 特征值 回传</param> /// <

[算法系列之二十六]字符串匹配之KMP算法

一 简介 KMP算法是一种改进的字符串匹配算法,由D.E.Knuth与V.R.Pratt和J.H.Morris同时发现,因此人们称它为克努特-莫里斯-普拉特操作(简称KMP算法).KMP算法的关键是利用匹配失败后的信息,尽量减少模式串与主串的匹配次数以达到快速匹配的目的. 二 基于部分匹配表的KMP算法 举例来说,有一个字符串"BBC ABCDAB ABCDABCDABDE",我想知道,里面是否包含搜索串"ABCDABD"? 步骤1:字符串"BBC ABC

Java加密技术（二）——对称加密算法DES&AES

接下来我们介绍对称加密算法,最常用的莫过于DES数据加密算法. DES DES-Data Encryption Standard,即数据加密算法.是IBM公司于1975年研究成功并公开发表的.DES算法的入口参数有三个:Key.Data.Mode.其中Key为8个字节共64位,是DES算法的工作密钥;Data也为8个字节64位,是要被加密或被解密的数据;Mode为DES的工作方式,有两种:加密或解密. DES算法把64位的明文输入块变为64位的密文输出块,它所使用的密钥也是64位. 通过java

堆排序：什么是堆?什么是最大堆?二叉堆是什么？堆排序算法是怎么样的？PHP如何实现堆排序？

本文标签:  堆排序 php php算法 堆排序算法 二叉堆 数据结构 REST   服务器 什么是堆 这里的堆(二叉堆),指得不是堆栈的那个堆,而是一种数据结构. 堆可以视为一棵完全的二叉树,完全二叉树的一个"优秀"的性质是,除了最底层之外,每一层都是满的,这使得堆可以利用数组来表示,每一个结点对应数组中的一个元素. 数组与堆之间的关系 二叉堆一般分为两种:最大堆和最小堆. 什么是最大堆 堆中每个父节点的元素值都大于等于其孩子结点(如果存在),这样的堆就是一个最大堆 因此,最大堆中的

Linux-0.11源代码阅读二 实模式到保护模式

bootsect部分已经执行完成,程序也跳转到setup部分: start: ! ok, the read went well so we get current cursor position and save it for ! posterity. mov ax,#INITSEG ! this is done in bootsect already, but... mov ds,ax mov ah,#0x03 ! read cursor pos xor bh,bh int 0x10 ! sa