矩阵求逆 QR、Cholesky 分解

【参考】

[1] 视觉slam14讲

[2] 维基百科

[3] [数值计算] QR分解

[4] QR分解

[5] Cholesky分解

目录

一、QR分解

二、Cholesky分解

三、c++实现

---

，根据m，n的相对大小，可以分为三种情况

 

 

一、QR分解

QR分解法 是Francis于1961年发表的用于求解所有特征值的算法，把矩阵分解成一个正交矩阵与一个上三角矩阵的积。QR分解经常用来解线性最小二乘法问题。QR分解也是特定特征值算法即QR算法的基础。

1.1 定义

一个矩阵可以被分解成A=QR，其中：

对于一个over-determined线性最小二乘问题，其目标函数是

，

  

其中，

由于等式右边第二项是常数，能最小化只有第一项，当时，最小值为

 

1.2 三种计算方法

• Gram–Schmidt Orthogonalization（格拉姆-施密特正交化方法）
• Householder Triangularization（Householder变换）
• Givens Rotations（吉文斯旋转）

 

二、Cholesky分解

线性代数中，Cholesky分解（英语：Cholesky decomposition or Cholesky factorization）是指将一个正定的埃尔米特矩阵分解成一个下三角矩阵与其共轭转置之乘积。这种分解方式在提高代数运算效率、蒙特卡罗方法等场合中十分有用。实数矩阵的Cholesky分解由安德烈·路易·科列斯基最先发明。实际应用中，Cholesky分解在求解线性方程组中的效率约两倍于LU分解。

 

 

(未完待续)

三、c++实现

#include <iostream>
#include <ctime>

#include <Eigen/Core>
#include <Eigen/Dense>

using namespace Eigen;

#define MATRIX_SIZE 50

using namespace std;

int main(int argc, char **argv){
    //定义变量
    Matrix<float,2,3> matrix_23;
    
    Vector3d v_3d;//与下面一样是列向量 double
    Matrix<float,3,1> vd_3d;//lie
    
    Matrix3d matrix_33 = Matrix3d::Zero();
    Matrix<double,Dynamic,Dynamic> matrix_dynamic;
    MatrixXd matrix_x;
    
    //初始化变量
    matrix_23 << 1,2,3,4,5,6;
    v_3d << 3, 2, 1;
    vd_3d << 4, 5, 6;
    
    cout<<"matrix_23:\n"<<matrix_23<<endl;
    cout<<"v_3d:\n"<<v_3d<<endl;
    cout<<"vd_3d:\n"<<vd_3d<<endl;
    for(int i=0;i<2;i++)
        for(int j=0;j<3;j++)
            cout<<matrix_23(i,j)<<"\t";
    cout<<endl;
    
    //一个float，一个double，不能直接相乘，应该显式转换
    Matrix<double, 2, 1> result = matrix_23.cast<double>() * v_3d;
    cout << "[1,2,3;4,5,6]*[3,2,1]=" << result.transpose() << endl;
    //一个float，一个float
    Matrix<float, 2, 1> result2 = matrix_23 * vd_3d;
    cout << "[1,2,3;4,5,6]*[4,5,6]: " << result2.transpose() << endl;
    
    matrix_33 = Matrix3d::Random();      // 随机数矩阵
    cout << "random matrix: \n" << matrix_33 << endl;
    cout << "transpose: \n" << matrix_33.transpose() << endl;      // 转置
    cout << "sum: " << matrix_33.sum() << endl;            // 各元素和
    cout << "trace: " << matrix_33.trace() << endl;          // 迹
    cout << "times 10: \n" << 10 * matrix_33 << endl;               // 数乘
    cout << "inverse: \n" << matrix_33.inverse() << endl;        // 逆
    cout << "det: " << matrix_33.determinant() << endl;    // 行列式
    
    
    // 特征值
    // 实对称矩阵可以保证对角化成功
    SelfAdjointEigenSolver<Matrix3d> eigen_solver(matrix_33.transpose() * matrix_33);
    cout << "Eigen values = \n" << eigen_solver.eigenvalues() << endl;
    cout << "Eigen vectors = \n" << eigen_solver.eigenvectors() << endl;
    
    
    Matrix<double, MATRIX_SIZE, MATRIX_SIZE> matrix_NN
      = MatrixXd::Random(MATRIX_SIZE, MATRIX_SIZE);
    matrix_NN = matrix_NN * matrix_NN.transpose();  // 保证半正定
    Matrix<double, MATRIX_SIZE, 1> v_Nd = MatrixXd::Random(MATRIX_SIZE, 1);

    clock_t time_stt = clock(); // 计时
    // 直接求逆
    Matrix<double, MATRIX_SIZE, 1> x = matrix_NN.inverse() * v_Nd;
    cout << "time of normal inverse is "
        << 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
    cout << "x = " << x.transpose() << endl;

    // 通常用矩阵分解来求，例如QR分解，速度会快很多
    time_stt = clock();
    x = matrix_NN.colPivHouseholderQr().solve(v_Nd);
    cout << "time of Qr decomposition is "
        << 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
    cout << "x = " << x.transpose() << endl;

    // 对于正定矩阵，还可以用cholesky分解来解方程
    time_stt = clock();
    x = matrix_NN.ldlt().solve(v_Nd);
    cout << "time of ldlt decomposition is "
        << 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
    cout << "x = " << x.transpose() << endl;
    return 0;
}

random matrix: 
 0.680375   0.59688 -0.329554
-0.211234  0.823295  0.536459
 0.566198 -0.604897 -0.444451

time of normal inverse is 0.299ms

time of Qr decomposition is 0.174ms

time of cholesky decomposition is 0.077ms