SVD 最小二乘法解 亲测ok!_svd.matrixu()

线性最小二乘问题

m个方程求解n个未知数，有三种情况：

1. m=n且A为非奇异，则有唯一解，x=A.inverse()*b
2. m>n，约束的个数大于未知数的个数，称为超定问题（overdetermined）
3. m<n，负定/欠定问题（underdetermined）

通常我们遇到的都是超定问题，此时Ax=b的解是不存在的，从而转向解最小二乘问题：

J(x)为凸函数，一阶导数为0，得到：

,称之为正规方程

一般解：

奇异值分解与线性最小二乘问题

设

列满秩，A的奇异值分解:（公式1）

其中U和V为半酉阵，分别满足

其中，这里的几个符号的大小：U:mm   :nn   V:nn

注意：为什么这里是nn，是因为我上面写的公式1中UGV，G=，G的尺寸是mn 

为 U 的 前 n 列矩阵，即

则：

等号当且仅

时成立，所以：

这就是我们千辛万苦要求的线性最小二乘问题的解！！！

用eigen库计算的例子：

    //Ax=b
    MatrixXd A= MatrixXd::Zero(15,8);
    Eigen::JacobiSVD<MatrixXd> svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV);
 
    Eigen::Matrix<double, 15, 15> U = svd.matrixU();
    Eigen::Matrix<double, 8, 8> V = svd.matrixV();
    Eigen::Matrix<double, 8, 8> Gama = svd.singularValues().asDiagonal();
        
    Eigen::Matrix<double, 15, 1> b;
    Eigen::Matrix<double, 8, 1> x = V * Gama.inverse()*(U.block<15,8>(0,0).transpose())*b;

 svd解法：修改 A 和 data（代表XYZ）

         修改未知量的个数；

//解决问题： Ax=XYZ
{
	int numberOfpara = 3;
    std::vector<double> data = { 1.7,1.6,1.6,1.5,1.2,1.3,2.3,2,1.7,1.2};
    typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> MatrixXd;
    MatrixXd A;
    A.resize(data.size(), numberOfpara);
    A.fill(0);
    /*A << 1.2, 0.1, 1.2, 2.3, 1.4,
        1.12, 1.03, 2, 2.06, 1.4,
        1.32, 1.24, 2.1, 2.05, 1.3,
        1.4, 0.77, 1.5, 2.02, 1.3,
        1.3, 1.37, 2, 2, 1.2,
        1.16, 0.8, 1.7, 2, 1.3,
        1.9, 1.23, 1.8, 1.99, 1.5,
        1.65, 1.27, 1.9, 1.99, 1.4,
        1.09, 0.89, 2, 1.99, 1.3,
        1.5, 0.75, 1.1, 1.99, 1.2;*/ 
    cout << "Here is the matrix m:" << endl << A << endl;
 
    int row=A.rows();   
    int col = A.cols();
 
    typedef Eigen::Matrix<double, Eigen::Dynamic, 1> VectorXd;
    Eigen::Map<VectorXd> XYZ(data.data(), row, 1);
    cout << "Here is the matrix XYZ:" << endl << XYZ << endl;
 
   // svd解法
    Eigen::JacobiSVD<MatrixXd> svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV);
 
    MatrixXd U, V, Gama, x;
    U.resize(row, row);
    V.resize(col, col);
    Gama.resize(col, col);
    x.resize(col, 1);
 
     U = svd.matrixU();
     V = svd.matrixV();
     Gama = svd.singularValues().asDiagonal();
 
    x = V * Gama.inverse()*(U.block(0, 0, row, col).transpose())*XYZ;    
    cout << "Here is the x:" << endl << x << endl;
 
    //验证结果
    cout << "Here is the eps:" << endl << A * x - XYZ;
    //x转化为vector
    std::vector<double> out(&x(0, 0), x.data() + x.cols()*x.rows());