视觉SLAM十四讲读书笔记（5）P40-P52

Eigen是一个C++模板库，提供了线性代数运算的基本功能，包括矩阵、向量、数值分解、求解线性方程组等等。它的设计目标是高性能、易用性和通用性。Eigen提供了类似于MATLAB和NumPy的API，可以简单易懂地进行各种矩阵和向量计算。它还提供了多种算法来处理各种矩阵运算，比如求解特征值、SVD分解、QR分解等。Eigen是一个开源库，可以自由地使用和分发。它使用头文件实现，因此很容易将其集成到

文章共8,283字 · 阅读需要大约28分钟

一键AI生成摘要，助你高效阅读

问答

嘉子的秃头日记

2042人浏览 · 2023-04-12 22:01:22

嘉子的秃头日记 · 2023-04-12 22:01:22 发布

Q：什么是eigen

Eigen是一个C++模板库，提供了线性代数运算的基本功能，包括矩阵、向量、数值分解、求解线性方程组等等。它的设计目标是高性能、易用性和通用性。

Eigen提供了类似于MATLAB和NumPy的API，可以简单易懂地进行各种矩阵和向量计算。它还提供了多种算法来处理各种矩阵运算，比如求解特征值、SVD分解、QR分解等。

Eigen是一个开源库，可以自由地使用和分发。它使用头文件实现，因此很容易将其集成到您的项目中。它支持各种平台和编译器，包括Windows、Linux、Mac OS和多种编译器，如GCC、Clang、MSVC等。

由于Eigen的设计优秀，实现高效，因此在很多领域得到广泛应用，包括计算机图形学、机器人学、信号处理等等。

Q：正交矩阵的定义和性质

正交矩阵是指矩阵的列向量（或行向量）两两正交且模长为1的矩阵。具体地，设Q是一个n阶矩阵，则Q为正交矩阵当且仅当满足以下两个条件：

Q^TQ=QQ^T=I_n，其中I_n为n阶的单位矩阵。
Q的每一列（或每一行）都是单位向量，即|q_i|=1，其中q_i为Q的第i列（或第i行）。

正交矩阵的性质如下：

正交矩阵的逆矩阵是其转置矩阵，即Q^{-1}=Q^T。
正交矩阵的行列式的绝对值为1，即|det(Q)|=1。
正交矩阵的列向量（或行向量）组成一个标准正交基，即任何向量都可以由这组向量线性表示。
正交矩阵保持向量的长度和角度不变，即对于任意向量x，有|Qx|=|x|，且angle(Qx,Qy)=angle(x,y)。
正交矩阵的转置矩阵仍然是正交矩阵。
两个正交矩阵的乘积仍然是正交矩阵。

正交矩阵在许多数学和应用领域中都有重要的应用，例如在线性代数中用于求解特征值和特征向量，以及在计算机图形学中用于变换和旋转等操作。

Q：什么是群

群（Group）是代数学中的一个基本概念，是一种带有代数结构的集合。群的定义包含了一个二元运算和一些基本公理，它们决定了群的性质。

具体地说，设G是一个集合， $\circ$ 是一个在G上的二元运算，若满足以下四个基本公理：

封闭性：对于任意的 $a,b \in G$ ，有 $a \circ b \in G$ 。
结合律：对于任意的 $a,b,c \in G$ ，有 $(a \circ b) \circ c = a \circ (b \circ c)$ 。
单位元素：存在一个元素 $e \in G$ ，称为单位元素，满足对于任意的 $a \in G$ ，有 $a \circ e = e \circ a = a$ 。
逆元素：对于任意的 $a \in G$ ，存在一个元素 $a^{-1} \in G$ ，称为 $a$ 的逆元素，满足 $a \circ a^{-1} = a^{-1} \circ a = e$ 。

则称 $(G, \circ)$ 为一个群。如果群 $(G, \circ)$ 中的运算满足交换律，即对于任意的 $a,b \in G$ ，有 $a \circ b = b \circ a$ ，则称 $(G, \circ)$ 为交换群或阿贝尔群。

Q：什么是特殊正交群

特殊正交群（Special Orthogonal Group）通常记作SO(n)，是由所有n阶行列式为1的正交矩阵所组成的群。其中的元素都是n阶的矩阵，满足矩阵的逆等于它的转置矩阵。在二维中，特殊正交群是一个圆周，因为旋转矩阵可以由一个角度参数表示；在三维中，它是一个球面，因为旋转矩阵可以由三个角度参数表示。SO(n)群在物理学中经常用来描述空间的对称性，例如在旋转对称的体系中。

Q：什么叫齐次坐标

在计算几何中，齐次坐标是一种向量表示方法，它在计算和描述点、直线、平面等几何对象的位置和运动时非常有用。

齐次坐标是通过在原有的笛卡尔坐标系中增加一个维度，将二维或三维的笛卡尔坐标扩展为三维或四维的向量，从而实现点、直线、平面等几何对象的无限远点的表示。

在齐次坐标中，如果一个点的坐标为 $(x,y)$ ，则其齐次坐标为 $(x', y', w') = (wx, wy, w)$ ，其中 $w$ 为一个非零常数。这个齐次坐标表示的实际上是一个与原来点相同的等价点，因为当 $w$ 不等于 0 时， $(x', y', w')$ 与 $(x'/w', y'/w')$ 表示的点相同。

齐次坐标中的“齐次”一词源自于拉丁语“homogeneus”，表示“相似的、相同的、同种的”。这是因为在齐次坐标中，点与向量之间的运算可以看作是“同种”对象之间的运算，即点与点之间的运算、向量与向量之间的运算。

Q：什么是特殊欧氏群

特殊欧氏群（Special Euclidean Group）是指在欧几里得空间中保持长度和角度不变的所有等距变换所组成的群，通常用SE(n)表示。这个群包括旋转和平移两种基本变换，即旋转矩阵和平移向量构成的变换矩阵，满足矩阵的行列式为1。在二维欧氏空间中，特殊欧氏群就是平面上所有旋转和平移的组合。在三维欧氏空间中，特殊欧氏群则包括旋转和平移的组合。这个群在计算机图形学和机器人学中应用广泛，例如用于描述和控制三维物体的运动和变换。

Q：什么是运算符重载

运算符重载是指在编程中对原有运算符进行重新定义，以实现自定义类型的运算操作。例如，对于内置类型的运算符如“+”、“-”、“*”、“/”等，它们可以被重载，使其能够适用于自定义类型，实现类似加减乘除等运算。

通过运算符重载，可以让自定义类型的代码更加直观、易读和易用，提高代码的可读性和可维护性。重载的运算符可以像内置运算符一样使用，使代码看起来更像自然语言的形式。运算符重载在 C++、Java、Python、C#等编程语言中都得到了广泛应用。

Q：什么是g2o

g2o是一种用于求解非线性最小二乘问题的C++库。它广泛应用于机器人、计算机视觉、SLAM（Simultaneous Localization and Mapping）等领域中。g2o库提供了一种灵活的方式来描述问题的模型，可以使用它来优化包括位姿、相机参数、三维重建、图像匹配等方面的参数。g2o库使用了稀疏矩阵和自动微分等技术，具有高效的计算速度和内存利用率。由于其优秀的性能和广泛应用，g2o已经成为了SLAM领域中最受欢迎的求解器之一。

g2o的官网链接是：GitHub - RainerKuemmerle/g2o: g2o: A General Framework for Graph Optimization

Q：什么是sophus

Sophus是一个用于计算机视觉和机器人领域的C++库，实现了李代数和李群相关的计算和优化，以及常见的位姿表示（例如四元数、欧拉角、转子）。它主要被用于求解非线性优化问题，例如在视觉SLAM（Simultaneous Localization and Mapping）和机器人导航中。Sophus的优点之一是可以处理李代数和李群的数学特性，如非线性性和流形结构，从而提高了算法的数值稳定性和精度。

sophus的官网链接为：GitHub - strasdat/Sophus: C++ implementation of Lie Groups using Eigen.

cmake_minimum_required(VERSION 2.8)
project(useEigen)

set(CMAKE_BUILD_TYPE "Release")
set(CMAKE_CXX_FLAGS "-O3")

# 添加Eigen头文件
include_directories("/usr/include/eigen3")
add_executable(eigenMatrix eigenMatrix.cpp)

这段代码是一个使用Eigen库的C++项目的CMakeLists.txt文件。该文件指定了项目名称为"useEigen"，设置了C++编译器的优化级别为O3，指定了Eigen头文件的路径，并将其添加到项目的include路径中。同时，它还定义了一个名为"eigenMatrix"的可执行文件，并将名为"eigenMatrix.cpp"的源文件与之关联。

#include <iostream>

using namespace std;

#include <ctime>
// Eigen 核心部分
#include <Eigen/Core>
// 稠密矩阵的代数运算（逆，特征值等）
#include <Eigen/Dense>

这段代码包含了四个头文件的引用：

iostream：C++标准库中的输入输出流库，用于输入输出操作。
ctime：C++标准库中的时间库，用于获取当前系统时间等操作。
Eigen/Core：Eigen库的核心部分，包含了矩阵和向量等数学基础操作。
Eigen/Dense：Eigen库的稠密矩阵代数运算库，包含了矩阵分解、求逆、求特征值等高级数学操作。

using namespace Eigen;

#define MATRIX_SIZE 50

这段代码定义了一个常数 MATRIX_SIZE 并使用 using namespace Eigen; 将 Eigen 命名空间中的类型和函数引入到当前命名空间中。其中，Eigen 是一个 C++ 模板库，用于线性代数运算和矩阵计算。接下来的代码中将使用 Eigen 中的矩阵类型和运算。

int main(int argc, char **argv) {
  // Eigen 中所有向量和矩阵都是Eigen::Matrix，它是一个模板类。它的前三个参数为：数据类型，行，列
  // 声明一个2*3的float矩阵
  Matrix<float, 2, 3> matrix_23;

这行代码声明了一个类型为 Matrix<float, 2, 3> 的对象 matrix_23，该对象是一个 2x3 的矩阵，元素类型为 float。Matrix 是 Eigen 中的一个模板类，第一个模板参数指定矩阵元素的类型，第二个和第三个模板参数分别指定矩阵的行数和列数。

  // 同时，Eigen 通过 typedef 提供了许多内置类型，不过底层仍是Eigen::Matrix
  // 例如 Vector3d 实质上是 Eigen::Matrix<double, 3, 1>，即三维向量
  Vector3d v_3d;
  // 这是一样的
  Matrix<float, 3, 1> vd_3d;

这段代码中，我们声明了一个类型为 Vector3d 的三维向量 v_3d，和一个类型为 Matrix<float, 3, 1> 的三维向量 vd_3d，它们的本质是相同的，都是三维向量。由于 Eigen 通过 typedef 提供了许多内置类型，因此 Vector3d 可以看作是 Eigen::Matrix<double, 3, 1> 的别名，而 vd_3d 实质上就是一个列向量，它的行数为 3，列数为 1。

  // Matrix3d 实质上是 Eigen::Matrix<double, 3, 3>
  Matrix3d matrix_33 = Matrix3d::Zero(); //初始化为零
  // 如果不确定矩阵大小，可以使用动态大小的矩阵
  Matrix<double, Dynamic, Dynamic> matrix_dynamic;
  // 更简单的
  MatrixXd matrix_x;
  // 这种类型还有很多，我们不一一列举

这段代码使用了Eigen库定义的不同矩阵类型。例如，Matrix3d实际上是Eigen::Matrix<double, 3, 3>，用于表示3x3的矩阵。Matrix3d matrix_33 = Matrix3d::Zero()创建了一个3x3的零矩阵。如果不确定矩阵大小，可以使用动态大小的矩阵，例如Matrix<double, Dynamic, Dynamic>或更简单地使用MatrixXd。Eigen库支持许多其他的矩阵和向量类型，可以根据需要进行选择。

在Eigen中，Dynamic是一个特殊的值，表示一个矩阵或向量的维度是在运行时确定的。这意味着我们可以使用变量来指定矩阵或向量的大小，而不是在编译时指定固定的大小。

使用动态大小的矩阵的优点是它们非常灵活，可以在运行时根据需要动态分配内存。这使得它们非常适合需要动态调整大小的情况，如读取文件中的数据，或在运行时从用户获取输入。

// 下面是对Eigen阵的操作
  // 输入数据（初始化）
  matrix_23 << 1, 2, 3, 4, 5, 6;
  // 输出
  cout << "matrix 2x3 from 1 to 6: \n" << matrix_23 << endl;

  // 用()访问矩阵中的元素
  cout << "print matrix 2x3: " << endl;
  for (int i = 0; i < 2; i++) {
    for (int j = 0; j < 3; j++) cout << matrix_23(i, j) << "\t";
    cout << endl;
  }

以上代码是对 Eigen 矩阵的操作，具体内容如下：

matrix_23 << 1, 2, 3, 4, 5, 6;：将 1 到 6 的数据依次填入一个 2 行 3 列的矩阵 matrix_23 中，也就是：

1 2 3
4 5 6

cout << "matrix 2x3 from 1 to 6: \n" << matrix_23 << endl;：输出这个矩阵的内容，输出结果为：

matrix 2x3 from 1 to 6: 
1 2 3
4 5 6

for 循环中的语句则是分别输出矩阵中每一个元素的值，结果与上面的输出是一样的，即：

print matrix 2x3: 
1 2 3	
4 5 6

  // 矩阵和向量相乘（实际上仍是矩阵和矩阵）
  v_3d << 3, 2, 1;
  vd_3d << 4, 5, 6;

  // 但是在Eigen里你不能混合两种不同类型的矩阵，像这样是错的
  // Matrix<double, 2, 1> result_wrong_type = matrix_23 * v_3d;
  // 应该显式转换
  Matrix<double, 2, 1> result = matrix_23.cast<double>() * v_3d;
  cout << "[1,2,3;4,5,6]*[3,2,1]=" << result.transpose() << endl;

  Matrix<float, 2, 1> result2 = matrix_23 * vd_3d;
  cout << "[1,2,3;4,5,6]*[4,5,6]: " << result2.transpose() << endl;

这段代码是演示如何对矩阵和向量进行乘法运算。在Eigen中，矩阵和向量相乘仍然是矩阵和矩阵相乘。首先，定义了两个向量v_3d和vd_3d，然后使用cast()方法将matrix_23矩阵中的元素转换为double类型，与v_3d向量相乘，得到result向量。然后，matrix_23矩阵与vd_3d向量相乘，得到result2向量。最后，使用transpose()方法将向量转置，以便于输出。需要注意的是，在Eigen中，不同类型的矩阵和向量不能混合相乘，需要进行类型转换。

矩阵 matrix_23 的数据类型为 Matrix<float, 2, 3>，向量 v_3d 的数据类型为 Vector3d，这两个类型是不同的。在 Eigen 中，对于矩阵乘法操作，两个矩阵的数据类型必须相同，否则编译器将会报错。

在这种情况下，需要将 matrix_23 的数据类型显式地转换为 Matrix<double, 2, 3>，这样才能与向量 v_3d 相乘。具体做法是使用 cast 函数将 matrix_23 中的每个元素转换为 double 类型。因此正确的写法是：Matrix<double, 2, 1> result = matrix_23.cast<double>() * v_3d;

matrix_23.cast<double>()可以将matrix_23矩阵中的元素类型转换为double类型，因为矩阵matrix_23是Matrix<float, 2, 3>类型的，其中元素类型是float，而向量v_3d是Vector3d类型的，其中元素类型是double，两者不同类型的矩阵和向量不能直接相乘。

因此，需要使用cast方法将矩阵matrix_23中的元素类型转换为double，得到一个Matrix<double, 2, 3>类型的矩阵，才能和Vector3d类型的向量v_3d相乘。

在代码中，v_3d 是 Vector3d 类型的向量，元素类型是 double，而 vd_3d 是 Matrix<float, 3, 1> 类型的向量，元素类型是 float。

  // 同样你不能搞错矩阵的维度
  // 试着取消下面的注释，看看Eigen会报什么错
  // Eigen::Matrix<double, 2, 3> result_wrong_dimension = matrix_23.cast<double>() * v_3d;

matrix_23 是一个 Matrix<float, 2, 3> 类型的矩阵，v_3d 是一个 Vector3d 类型的向量。这两者可以相乘，因为它们的维度匹配，但是得到的结果需要是一个 Matrix<double, 2, 1> 类型的矩阵。

Matrix<double, 2, 1> result_wrong_type = matrix_23 * v_3d 的问题在于，左侧矩阵的维度是不匹配的。

  // 一些矩阵运算
  // 四则运算就不演示了，直接用+-*/即可。
  matrix_33 = Matrix3d::Random();      // 随机数矩阵
  cout << "random matrix: \n" << matrix_33 << endl;
  cout << "transpose: \n" << matrix_33.transpose() << endl;      // 转置
  cout << "sum: " << matrix_33.sum() << endl;            // 各元素和
  cout << "trace: " << matrix_33.trace() << endl;          // 迹
  cout << "times 10: \n" << 10 * matrix_33 << endl;               // 数乘
  cout << "inverse: \n" << matrix_33.inverse() << endl;        // 逆
  cout << "det: " << matrix_33.determinant() << endl;    // 行列式

这段代码展示了一些常用的矩阵运算，下面是一些解释：

matrix_33 = Matrix3d::Random();：使用Random()函数生成一个随机数矩阵matrix_33，它的元素在[-1, 1]之间。
matrix_33.transpose()：返回矩阵matrix_33的转置矩阵。
matrix_33.sum()：返回矩阵matrix_33中所有元素的和。
matrix_33.trace()：返回矩阵matrix_33的迹，即主对角线上元素的和。
10 * matrix_33：对矩阵matrix_33进行数乘，即矩阵中的每个元素乘以10。
matrix_33.inverse()：返回矩阵matrix_33的逆矩阵。如果矩阵不可逆，会抛出一个Eigen::Exception异常。
matrix_33.determinant()：返回矩阵matrix_33的行列式。如果矩阵不可逆，返回0。

这些运算在线性代数中经常用到，尤其是在求解方程组和特征值、特征向量等问题中。

下面是一组输出示例：

random matrix: 
-0.981616   0.396752    1.32985
-0.0557061 -0.772753 -0.0692563
 0.236567  -0.471076   0.975621
transpose: 
-0.981616 -0.0557061  0.236567
 0.396752 -0.772753  -0.471076
  1.32985 -0.0692563  0.975621
sum: -0.361194
trace: -0.779747
times 10: 
-9.81616   3.96752    13.2985
-0.557061 -7.72753  -0.692563
 2.36567  -4.71076    9.75621
inverse: 
  -0.46513    1.20915  -0.100626
   0.16754  -0.183321  -0.369492
  -0.05287   0.644857    0.09406
det: -0.0352479

  // 特征值
  // 实对称矩阵可以保证对角化成功
  SelfAdjointEigenSolver<Matrix3d> eigen_solver(matrix_33.transpose() * matrix_33);
  cout << "Eigen values = \n" << eigen_solver.eigenvalues() << endl;
  cout << "Eigen vectors = \n" << eigen_solver.eigenvectors() << endl;

这段代码计算了实对称矩阵 matrix_33.transpose() * matrix_33 的特征值和特征向量，输出了特征值和特征向量的矩阵形式。

输出的特征值矩阵是一个列向量，每个元素表示对应的特征值。输出的特征向量矩阵是一个矩阵，每一列表示对应的特征向量。特征值和特征向量按照特征值的大小从大到小排列。

C++中可以通过定义类来封装数据和方法，并通过创建对象来访问它们。语法格式是：类名对象名；其中，类名表示定义好的类，对象名表示你创建的对象，可以调用该类中的成员函数和成员变量。

SelfAdjointEigenSolver<Matrix3d> eigen_solver(matrix_33.transpose() * matrix_33);

这行代码使用了 SelfAdjointEigenSolver 模板类，它是 Eigen 中专门用于对实对称矩阵求特征值和特征向量的类。它的模板参数是矩阵的类型。在这里，参数类型是 Matrix3d，表示一个 $3\times 3$ 的双精度实数矩阵。

该行代码的实际操作是对矩阵 matrix_33.transpose() * matrix_33 求特征值和特征向量。由于该矩阵是实对称矩阵，因此可以通过 SelfAdjointEigenSolver 求解，得到的特征值和特征向量是实数类型的，分别存储在 eigen_solver.eigenvalues() 和 eigen_solver.eigenvectors() 中。

SelfAdjointEigenSolver<MatrixType> eigen_solver(A);

其中，MatrixType是矩阵的类型，A是需要求解特征值和特征向量的矩阵，eigen_solver是一个SelfAdjointEigenSolver类型的对象，用于存储求解的结果。

这个语句定义了一个名为eigen_solver的对象，类型为SelfAdjointEigenSolver<MatrixType>，其中MatrixType是一个模板参数，表示矩阵类型。这个对象的构造函数使用一个矩阵A作为参数，用于计算这个矩阵的特征值和特征向量。

  // 解方程
  // 我们求解 matrix_NN * x = v_Nd 这个方程
  // N的大小在前边的宏里定义，它由随机数生成
  // 直接求逆自然是最直接的，但是求逆运算量大

  Matrix<double, MATRIX_SIZE, MATRIX_SIZE> matrix_NN
      = MatrixXd::Random(MATRIX_SIZE, MATRIX_SIZE);
  matrix_NN = matrix_NN * matrix_NN.transpose();  // 保证半正定
  Matrix<double, MATRIX_SIZE, 1> v_Nd = MatrixXd::Random(MATRIX_SIZE, 1);

这段代码定义了一个随机的 MATRIX_SIZE × MATRIX_SIZE 的矩阵 matrix_NN 和一个 MATRIX_SIZE × 1 的随机向量 v_Nd。其中，matrix_NN 先与其转置相乘，使其成为半正定矩阵。这是为了确保其可逆，从而方便求解方程。

MatrixXd是Eigen库中的一个类，表示动态大小的矩阵，即行数和列数可以在运行时确定。其构造函数的形式为MatrixXd(rows,cols)，其中rows和cols表示矩阵的行数和列数。使用MatrixXd可以方便地定义任意大小的矩阵。

MatrixXd::Random(MATRIX_SIZE, MATRIX_SIZE)创建一个MATRIX_SIZE x MATRIX_SIZE的矩阵，其中每个元素都是一个[-1, 1]之间的随机数。

MATRIX_SIZE 的值为 50。

一个 $n\times n$ 实矩阵A被称为半正定矩阵，如果对于任意n维非零实向量x，都有 $x^\mathrm{T}Ax\geq 0$ 。其中， $\mathrm{T}$ 表示矩阵的转置。

半正定矩阵具有以下性质：

半正定矩阵的特征值均为非负实数。
半正定矩阵可以被分解成 $A=LL^\mathrm{T}$ 的形式，其中L是一个下三角矩阵。
半正定矩阵是一类非常重要的矩阵，它们在统计学、机器学习和优化等领域有着广泛的应用，例如用于协方差矩阵的估计、二次规划问题的求解等等。

设A是一个矩阵，则 $A\cdot A^\top$ 是一个对称矩阵，因为 $(A\cdot A^\top)^\top=A\cdot A^\top$ 。又因为对于任意非零向量x，有 $x^\top\cdot A\cdot A^\top\cdot x=|A^\top\cdot x|^2\geq 0$ ，所以 $A\cdot A^\top$ 是半正定矩阵。

  clock_t time_stt = clock(); // 计时
  // 直接求逆
  Matrix<double, MATRIX_SIZE, 1> x = matrix_NN.inverse() * v_Nd;
  cout << "time of normal inverse is "
       << 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
  cout << "x = " << x.transpose() << endl;

这段代码是对矩阵求逆并求解方程 Ax=b，其中 A是半正定矩阵，b是列向量。

矩阵求逆的方法是使用 inverse() 函数。然而，这种方法在计算量大的时候效率比较低。

第二行代码的作用是求解线性方程组 Ax=b，其中 A 是 nxn的矩阵，b 是nx1 的列向量，x 是 nx1的列向量，即求解 x满足 Ax=b。在这里，A 是半正定矩阵，使用求逆的方式来求解线性方程组。

1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC 是一个计算程序运行时间的表达式。在这个表达式中，clock()函数返回从程序开始运行到调用该函数时为止的CPU时钟周期数，单位为“时钟打点”。CLOCKS_PER_SEC表示CPU每秒运行的时钟打点数。因此，(clock() - time_stt) / (double) CLOCKS_PER_SEC的结果就是程序运行的CPU时间（秒），乘以1000即可转换为毫秒。

  // 通常用矩阵分解来求，例如QR分解，速度会快很多
  time_stt = clock();
  x = matrix_NN.colPivHouseholderQr().solve(v_Nd);
  cout << "time of Qr decomposition is "
       << 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
  cout << "x = " << x.transpose() << endl;

这部分代码使用QR分解的方法来求解方程组，使用了矩阵库中的 colPivHouseholderQr() 函数对系数矩阵进行QR分解，然后使用 solve() 函数对分解后的矩阵和向量进行求解。

具体来说，QR分解是将一个矩阵分解为一个正交矩阵和一个上三角矩阵的乘积，即A=QR，其中Q为正交矩阵，R为上三角矩阵。对于一个满秩矩阵A，我们可以通过QR分解来求解线性方程组Ax=b。首先将系数矩阵A进行QR分解，即A=QR，然后将方程组变为 $Rx=Q^Tb$ ，再通过回带求解即可得到x的解。

相比直接求逆的方法，QR分解方法具有更高的求解效率和更好的数值稳定性。

  // 对于正定矩阵，还可以用cholesky分解来解方程
  time_stt = clock();
  x = matrix_NN.ldlt().solve(v_Nd);
  cout << "time of ldlt decomposition is "
       << 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
  cout << "x = " << x.transpose() << endl;

  return 0;

Cholesky分解是一种用于解决线性方程组的数值算法，对于对称正定矩阵，Cholesky分解可以将其分解为一个下三角矩阵L$和其转置矩阵 $L^T$ 的乘积，即 $A = L \cdot L^T$ 。其中L的对角线元素都是正数，因此可以将A分解为 $A = L \cdot L^T$ 的形式，这个过程可以减少计算量，提高计算速度。