对于直线拟合y=kx+b,构造方程组kx+b-y=0-》Ax=0,求解A^T*A的最小特征值对应的特征向量,即使问题的答案k,b,1,下面用C++ Eigen3来实现
#include <iostream>
#include <fstream>
#include <vector>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
//cout << "Start Program" << endl;
ifstream fread("data.txt");
if (fread.fail())
{
cout << "Data File does not exit." << endl;
return 0;
}
double x = 0;
double y = 0;
vector<double> xs;
vector<double> ys;
while (!fread.eof())
{
fread >> x >> y;
xs.push_back(x);
ys.push_back(y);
}
int len = xs.size();
//cout << len << endl;
Eigen::MatrixXd A(len, 3);
for (int i = 0; i < len; i++)
{
A(i, 0) = xs[i];
A(i, 1) = 1;
A(i, 2) = -ys[i];
}
cout << "Matrix A" << endl;
cout << A << endl;
Eigen::SelfAdjointEigenSolver < Eigen::Matrix<double, 3, 3> > eigen_sovler(A.transpose()*A); // 特征值求解器
cout << "Eigen Values : " << endl;
cout << eigen_sovler.eigenvalues().transpose() << endl; // 输出特征值
cout << "Eigen Vector :" << endl;
cout << eigen_sovler.eigenvectors() << endl; // 输出特征向量
double v[3];
v[0] = eigen_sovler.eigenvectors()(0, 0); // 提取第一列特征向量
v[1] = eigen_sovler.eigenvectors()(1, 0);
v[2] = eigen_sovler.eigenvectors()(2, 0);
v[0] /= v[2];
v[1] /= v[2];
v[2] /= v[2];
cout <<"k:"<< v[0] << " " <<"b:"
<< v[1] << " " << v[2] << " " << endl;
return 0;
}