Eigen求最小二乘解

71 阅读 0 评论 47 点赞

我是靠谱客的博主背后胡萝卜，最近开发中收集的这篇文章主要介绍Eigen求最小二乘解，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

一个超定方程组，例如 $A x = b$ ，没有解。在这种情况下，在 $A x - b$ 尽可能小的意义上，搜索最接近解的向量 $x$ 是有意义的。该 $x$ 称为最小二乘解（如果使用欧几里得范数）。

此处使用的三种求解方法为：SVD分解，QR分解，正规矩阵(normal equation)。在这三种方程中，SVD分解通常准确率最高，但是速度最慢；正规矩阵速度最快，但是准确率最低；QR分解位于两者之间。

求解问题

  Eigen::MatrixXd A = Eigen::MatrixXd::Random(3, 2);
  std::cout << "Here is the matrix A:n" << A << std::endl;
  Eigen::VectorXd b = Eigen::VectorXd::Random(3);
  std::cout << "Here is the right hand side b:n" << b.transpose() << std::endl;

Here is the matrix A:
 0.680375   0.59688
-0.211234  0.823295
 0.566198 -0.604897
Here is the right hand side b:
-0.329554  0.536459 -0.444451

SVD分解

可以直接使用BDCSVD类中的solve()方法求解线性方程。单纯计算奇异值是不足够的，还需要计算奇异向量，但是稀疏SVD分解足够计算最小二乘解。

  // using the SVD decomposition
  std::cout << "The solution using the SVD decomposition is:n"
            << A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV)
                   .solve(b)
                   .transpose()
            << std::endl;

The solution using the SVD decomposition is:
-0.669626  0.314253

QR分解

Eigen中QR分解类有三个：HouseholderQR（没有枢转，快速但不稳定），ColPivHouseholderQR（列枢转，因此有一些慢但是更准确），FullPivHouseholderQR（完全旋转，最慢，但是最稳定）。

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

The solution using the QR decomposition is:
-0.669626  0.314253

正规矩阵

求解 $A x = b$ 的最小二乘解等于求解正规矩阵 $A^TAx = A^Tb$ 。

  // using normal equations
  std::cout << "The solution using normal equations is:n"
            << (A.transpose() * A).ldlt().solve(A.transpose() * b).transpose()
            << std::endl;

The solution using normal equations is:
-0.669626  0.314253

如果矩阵 $A$ 是病态的，那么这不是一个好方法，因为 $A^TA$ 是 $A$ 的平方。这意为着使用正规方程所造成的损失是使用其他方法所造成损失的两倍。

代码

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

int main(int argc, char **argv) {
  Eigen::MatrixXd A = Eigen::MatrixXd::Random(3, 2);
  std::cout << "Here is the matrix A:n" << A << std::endl;
  Eigen::VectorXd b = Eigen::VectorXd::Random(3);
  std::cout << "Here is the right hand side b:n" << b.transpose() << std::endl;

  // using the SVD decomposition
  std::cout << "The solution using the SVD decomposition is:n"
            << A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV)
                   .solve(b)
                   .transpose()
            << std::endl;

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

  // using normal equations
  std::cout << "The solution using normal equations is:n"
            << (A.transpose() * A).ldlt().solve(A.transpose() * b).transpose()
            << std::endl;

  return 0;
}