C++通过Eigen库实现最小二乘法的三种方法

C++通过Eigen库实现最小二乘法的三种方法

1、最小二乘法的数学原理

http://image.mamicode.com/info/201804/20180420183835504092.png

2、矩阵伪逆的C++实现

Maltab中通过pinv函数可以实现矩阵伪逆的求解,那如何在C++中实现矩阵的计算呢,比如Ax=B,这种表达式是在最小二乘法中常见的表达式,如何通过最小二乘法求解出x呢?通过调用Eigen库即可实现最小二乘法。

方法一 :通过SVD分解实现

实现案例:

#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;

}

方法二 :通过QR分解实现

实现案例:

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 is equivalent to solving the normal equation ATAx = ATb

实现案例

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

VectorXf b = VectorXf::Random(3);

cout << "The solution using normal equations is:\n"

     << (A.transpose() * A).ldlt().solve(A.transpose() * b) << endl;

 

当矩阵A为病态矩阵时,通过常规表达式求解时效果不好。

 

总结

SVD分解方法最准确,但是运算速度最慢;常规求解方法运算速度快,但准确度低;QR分解在两者之间。

  • 3
    点赞
  • 23
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
使用Eigen实现最小二乘法拟合圆的步骤如下: 1. 定义数据点数组和目标变量数组。 ```c++ #include <Eigen/Dense> using namespace Eigen; int main() { // 定义数据点和目标变量 MatrixXd data(2, n); // 2维数据点,共n个点 VectorXd b(n); // 目标变量 // 将数据点和目标变量赋值 // ... } ``` 2. 构造设计矩阵。 ```c++ // 构造设计矩阵 MatrixXd A(n, 3); for (int i = 0; i < n; i++) { double x = data(0, i); double y = data(1, i); A(i, 0) = x * x; A(i, 1) = y * y; A(i, 2) = x * y; } ``` 3. 使用最小二乘法求解系数。 ```c++ // 求解系数 VectorXd x = A.bdcSvd(ComputeThinU | ComputeThinV).solve(b); ``` 4. 根据系数计算圆的参数。 ```c++ // 计算圆的参数 double a = x(0); double b = x(1); double c = x(2); double xc = b / (2 * c); double yc = a / (2 * c); double r = sqrt((a * b - c * c) / c); ``` 完整代码如下: ```c++ #include <Eigen/Dense> #include <iostream> using namespace Eigen; int main() { // 定义数据点和目标变量 int n = 10; MatrixXd data(2, n); VectorXd b(n); data << 0.5, 1.2, 2.0, 3.6, 4.1, 4.9, 5.5, 6.3, 7.0, 8.4, 0.9, 2.2, 3.6, 5.2, 4.8, 6.2, 7.4, 8.1, 9.0, 10.0; b << 1.2, 2.1, 3.7, 4.7, 5.5, 6.3, 7.4, 8.1, 9.0, 10.1; // 构造设计矩阵 MatrixXd A(n, 3); for (int i = 0; i < n; i++) { double x = data(0, i); double y = data(1, i); A(i, 0) = x * x; A(i, 1) = y * y; A(i, 2) = x * y; } // 求解系数 VectorXd x = A.bdcSvd(ComputeThinU | ComputeThinV).solve(b); // 计算圆的参数 double a = x(0); double b = x(1); double c = x(2); double xc = b / (2 * c); double yc = a / (2 * c); double r = sqrt((a * b - c * c) / c); // 输出结果 std::cout << "Center: (" << xc << ", " << yc << ")\n"; std::cout << "Radius: " << r << std::endl; return 0; } ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值