C++ 实现householder变换

C++ 实现householder变换

此次写householder是为了实现矩阵的QR分解。
householder有个重要的性质就是，他可以将一个向量，除了第一个值以外，通过一个H矩阵把其余的数都给消去：Hx=y。

代码引入了Eigen3库，对矩阵做基础操作，新建，分块，加减乘除等等。

以下是代码：
第一份代码，返回结果是变换以后的向量。那个fanshu（vec）函数就是范数的意思，直接用vec.norm()也行，我当时脑子抽筋自己写了一个，懒得改了。
第二份代码是返回H矩阵。

//householder变换出只有一个元素的向量
Eigen::VectorXd houseHolder(Eigen::VectorXd vec){

    //找出v向量，并求出w向量
    double firstNum=sign(vec(0))*fanShu(vec);
    Eigen::VectorXd newVec=Eigen::VectorXd::Zero(vec.size());
    newVec[0]=firstNum;
    Eigen::VectorXd v=newVec+vec;

    //计算Hx=x-2(vTx/vTv)v
    double vTx=v.transpose()*vec;
    double vTv=v.transpose()*v;
    double divideRes=vTx/vTv;

    //得出最终只有第一个数的列向量
    Eigen::VectorXd Hx=vec-2*divideRes*v;
    return Hx;
}

//householder变换得到 H=I-wwT
Eigen::MatrixXd houseHolder_H(Eigen::VectorXd vec){

    //找出v向量，并求出w向量
    double firstNum=sign(vec(0))*fanShu(vec);
    Eigen::VectorXd newVec=Eigen::VectorXd::Zero(vec.size());
    newVec[0]=firstNum;
    Eigen::VectorXd v=newVec+vec;

    //计算H=I-2vvT/(|v|^2)
    Eigen::MatrixXd vvT=v*v.transpose();
    double scale=2/v.squaredNorm();
    Eigen::MatrixXd I=Eigen::MatrixXd::Identity(vvT.rows(),vvT.cols());
    Eigen::MatrixXd H=I-scale*vvT;
    return H;
}

鄙人很菜，这是第一次用C++写东西，namespace啥的也没管，目前也没太搞明白。
看了两眼C++ premier的书放弃了，太厚了。直接上手，不懂再查。
能用就完了。