Eigen密集矩阵求解 2 - 求解最小二乘系统

简介

本篇介绍如何使用Eigen求解线性最小二乘系统。

一个系统可能无精确的解，比如Ax=b的线性方程式，不存在解。这时，找到一个最接近的解x，使得偏差Ax-b尽可能地小，能够满足误差要求error-margin。那这个x就称为最小二乘解。

这里讨论3个方法： SVD分解法，QR分解法，和规范等式。这中间，SVD分解法精度最高，但效率最差；规范式最快但精度最小；而QR分解法居中。

奇异值分解(SVD)法（Singular value decomposition）

使用Eigen中的BDCSVD类的的solve()方法，就能直接解出线性二乘系统了。但对奇异值计算，这样并不足够，你也会需要计算奇异向量。

示例如下，其使用了Matrix的bdcsvd()方法来创建一个BDCSVD类的实例：

#include <iostream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

int main()
{
   MatrixXf A = MatrixXf::Random(3, 2);
   cout << "Here is the matrix A:\n" << A << endl;

   VectorXf b = VectorXf::Random(3);
   cout << "Here is the right hand side b:\n" << b << endl;

   cout << "The least-squares solution is:\n"
        << A.bdcSvd(ComputeThinU | ComputeThinV).solve(b) << endl;
}

执行：

$ g++ -I /usr/local/include/eigen3  matrix_svd1.cpp -o matrix_svd1
$ 
$ ./matrix_svd1 
Here is the matrix A:
 -0.999984 -0.0826997
 -0.736924  0.0655345
  0.511211  -0.562082
Here is the right hand side b:
-0.905911
 0.357729
 0.358593
The least-squares solution is:
  0.46358
0.0429898

QR分解法

在Eigen中，QR分解类的solve()方法用于计算最小二乘解。Eigen内提供了3中QR分解类：

• HouseholderQR: 无需行列转换pivoting，速度快，但不稳定。
• ColPivHouseholderQR: 需要列转换，稍慢，但精度高。
• FullPivHouseholderQR: 完全的行列转换，所以最慢，但最稳定。

参考下面简单示例，使用A.colPivHouseholderQr().solve(b)：

MatrixXf A = MatrixXf::Random(3, 2);
VectorXf b = VectorXf::Random(3);

cout << "The solution using the QR decomposition is:\n"
     << A.colPivHouseholderQr().solve(b) << endl;

使用规范等式

有变换等式转换： $Ax=b$ ⇒ $A^{T}Ax=A^{T}b$， 这里对矩阵A有要求，$A^{T}A$结果为方阵，如果矩阵的条件比较病态，则计算可能会急剧变大，考虑MatrixXf(1000,2)，需要计算一个1000X1000的矩阵了。

下面示例：

//matrix_decom_norm.cpp
#include <iostream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

int main()
{
    MatrixXf A = MatrixXf::Random(3, 2);
    cout<<"A: "<<endl<<A<<endl;

    VectorXf b = VectorXf::Random(3);
    cout<<"b: "<<endl<<b<<endl;
    
    cout << "The solution using normal equations is:\n"
        << (A.transpose() * A).ldlt().solve(A.transpose() * b) << endl;
}

执行：

$ g++ -I /usr/local/include/eigen3  matrix_decom_norm.cpp -o matrix_decom_norm
$ ./matrix_decom_norm 
A: 
 -0.999984 -0.0826997
 -0.736924  0.0655345
  0.511211  -0.562082
b: 
-0.905911
 0.357729
 0.358593
The solution using normal equations is:
 0.463581
0.0429899