Eigen：进阶使用

最新推荐文章于 2024-11-25 09:00:00 发布

非晚非晚

最新推荐文章于 2024-11-25 09:00:00 发布

阅读量1.1k

点赞数

分类专栏：数学理论与工具文章标签：数学开发工具线性代数矩阵

本文链接：https://blog.youkuaiyun.com/QLeelq/article/details/111660907

版权

数学理论与工具专栏收录该内容

10 篇文章

订阅专栏

本文继续完善和拓展Eigen的功能使用，介绍了初始化操作，如逗号初始化、逗号连接初始化；特殊矩阵和数组操作，如全零操作等；还涉及范数计算、布尔约定、元素地址访问、行列访问，最后探讨了别名问题的出现、解决方法及不会出现别名问题的运算。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

我们上一篇文章对Eigen的基础使用进行了介绍，下面我们继续对Eigen的功能使用进行完善和拓展。

一. 初始化

我们上一篇文章对初始化进行了简单介绍，比如可以使用逗号操作进行初始化。下面分别介绍几种初始化操作。

1.1 逗号初始化

逗号初始化对矩阵、向量、数组都可以使用。

Matrix3f m;
m << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;
std::cout << m

初始化结果：

1 2 3
4 5 6
7 8 9

1.2 逗号连接初始化

连接初始化对矩阵、向量进行操作，不过需要注意的是需要提前定义好连接的大小。

向量连接

RowVectorXd vec1(3);
vec1 << 1, 2, 3;
std::cout << "vec1 = " << vec1 << std::endl;
 
RowVectorXd vec2(4);
vec2 << 1, 4, 9, 16;
std::cout << "vec2 = " << vec2 << std::endl;
 
RowVectorXd joined(7);//提前定义大小
joined << vec1, vec2;
std::cout << "joined = " << joined << std::endl;

输出结果：

vec1 = 1 2 3
vec2 =  1  4  9 16
joined =  1  2  3  1  4  9 16

矩阵连接

MatrixXf matA(2, 2);
matA << 1, 2, 3, 4;
MatrixXf matB(4, 4);
matB << matA, matA/10, matA/10, matA;
std::cout << matB << std::endl;

输出结果：

 1   2 0.1 0.2
  3   4 0.3 0.4
0.1 0.2   1   2
0.3 0.4   3   4

连接块初始化

Matrix3f m;
m.row(0) << 1, 2, 3;
m.block(1,0,2,2) << 4, 5, 7, 8;
m.col(2).tail(2) << 6, 9;                   
std::cout << m;

输出结果：

1 2 3
4 5 6
7 8 9

二. 特殊矩阵和数组

Eigen可以定义了特殊矩阵、数组操作，例如全零矩阵、固定数组矩阵。

2.1 全零操作

可以对矩阵和数组类进行操作。

zero()：只能用于固定大小对象
zero(n)：对一维进行操作
zero(m,n)：对2维进行操作

std::cout << "A fixed-size array:\n";
Array33f a1 = Array33f::Zero();
std::cout << a1 << "\n\n";
 
std::cout << "A one-dimensional dynamic-size array:\n";
ArrayXf a2 = ArrayXf::Zero(3);//3列为0
std::cout << a2 << "\n\n";
 
std::cout << "A two-dimensional dynamic-size array:\n";
ArrayXXf a3 = ArrayXXf::Zero(3, 4);
std::cout << a3 << "\n";

输出结果：

A fixed-size array:
0 0 0
0 0 0
0 0 0

A one-dimensional dynamic-size array:
0
0
0

A two-dimensional dynamic-size array:
0 0 0 0
0 0 0 0
0 0 0 0

2.2 constant、Random、Identity、LinSpaced操作

Constant(value)：矩阵或数组的元素全为value
Constant(rows,cols,value):指定矩阵元素的全value，用于动态数组（如下所示）
Random()：指定矩阵或数组的元素为随机数
Identity():指定矩阵的对角线的元素为1。（对数组不可操作）
LinSpaced(size, low, high):对向量或者一维数组进行操作，从low到high等分为size份。

代码举例：

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

int main()
{
    Matrix3d a = Matrix3d::Constant(3);
    cout<<"Constant(3):\n"<<a<<endl;
    Matrix3d b = Matrix3d::Random();
    cout<<"Random():\n"<<b<<endl;
    Matrix3d m = Matrix3d::Identity();
    cout<<"Identity():\n"<<m<<endl;

    ArrayXXf table(10, 4);
    table.col(0) = ArrayXf::LinSpaced(10, 0, 90);
    table.col(1) = M_PI / 180 * table.col(0);
    table.col(2) = table.col(1).sin();
    table.col(3) = table.col(1).cos();
    std::cout << "  Degrees   Radians      Sine    Cosine\n";
    std::cout << table << std::endl;

    Matrix3d c;
    c = MatrixXd::Constant(3,3,4);
    cout<<"Constant(3,3,4)\n"<<c<<endl;
}

输出结果：

Constant(3):
3 3 3
3 3 3
3 3 3
Random():
 0.680375   0.59688 -0.329554
-0.211234  0.823295  0.536459
 0.566198 -0.604897 -0.444451
Identity():
1 0 0
0 1 0
0 0 1
  Degrees   Radians      Sine    Cosine
           0            0            0            1
          10     0.174533     0.173648     0.984808
          20     0.349066      0.34202     0.939693
          30     0.523599          0.5     0.866025
          40     0.698132     0.642788     0.766044
          50     0.872665     0.766044     0.642788
          60       1.0472     0.866025          0.5
          70      1.22173     0.939693      0.34202
          80      1.39626     0.984808     0.173648
          90       1.5708            1 -4.37114e-08
Constant(3,3,4)
4 4 4
4 4 4
4 4 4

Eigen还提供相应的另外一个版本， setZero(), MatrixBase::setIdentity() and DenseBase::setLinSpaced().

代码举例：

const int size = 6;
MatrixXd mat1(size, size);
mat1.topLeftCorner(size/2, size/2)     = MatrixXd::Zero(size/2, size/2);
mat1.topRightCorner(size/2, size/2)    = MatrixXd::Identity(size/2, size/2);
mat1.bottomLeftCorner(size/2, size/2)  = MatrixXd::Identity(size/2, size/2);
mat1.bottomRightCorner(size/2, size/2) = MatrixXd::Zero(size/2, size/2);
std::cout << mat1 << std::endl << std::endl;
 
MatrixXd mat2(size, size);
mat2.topLeftCorner(size/2, size/2).setZero();
mat2.topRightCorner(size/2, size/2).setIdentity();
mat2.bottomLeftCorner(size/2, size/2).setIdentity();
mat2.bottomRightCorner(size/2, size/2).setZero();
std::cout << mat2 << std::endl << std::endl;
 
MatrixXd mat3(size, size);
mat3 << MatrixXd::Zero(size/2, size/2), MatrixXd::Identity(size/2, size/2),
        MatrixXd::Identity(size/2, size/2), MatrixXd::Zero(size/2, size/2);
std::cout << mat3 << std::endl;

输出结果：

0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0

0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0

0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0

三. 范数计算

3.1 范数的概念

先来介绍一些范数的概念

向量的范数

（1）1-范数：向量中所有元素的绝对值之和：
在这里插入图片描述
（2）2-范数：向量中每个元素的平方和再开根号，即欧式距离：

（3） $\infty$ -范数:即所有向量元素中的最大值。

（4） $-\infty$ -范数:即所有向量元素中的最小值。

（5）p-范数:即所有向量元素的p次方之和的1/p次幂。
在这里插入图片描述

矩阵的范数

假设矩阵大小为m行n列

（1）1-范数：也称列和范数。即矩阵列向量中绝对值之和的最大值。
在这里插入图片描述
（2）2-范数：又称谱范数，计算方法是 $A^{T}A$ 矩阵的最大特征值的平方根。

（3） $\infty$ -范数：又称行和范数。即矩阵行向量中绝对值之和的最大值。

（4）F-范数：Frobenius范数，计算方法就是矩阵中各个元素的绝对值的平方和在开方。
在这里插入图片描述

3.2范数的计算

squaredNorm() ：向量中所有元素的平方和。
norm()：2-范数，即squaredNorm()开根号。如果是矩阵，仅仅用在F-范数。
lpNorm<p>()：用于计算多维范数

#include <Eigen/Dense>
#include <iostream>
using namespace std;
using namespace Eigen;
int main()
{
  VectorXf v(2);
  MatrixXf m(2,2), n(2,2);

  v << -1,
       2;

  m << 1,-2,
       -3,4;

  cout << "v.squaredNorm() = " << v.squaredNorm() << endl;
  cout << "v.norm() = " << v.norm() << endl;
  cout << "v.lpNorm<1>() = " << v.lpNorm<1>() << endl;
  cout << "v.lpNorm<Infinity>() = " << v.lpNorm<Infinity>() << endl;

  cout << endl;
  cout << "m.squaredNorm() = " << m.squaredNorm() << endl;
  cout << "m.norm() = " << m.norm() << endl;
  cout << "m.lpNorm<1>() = " << m.lpNorm<1>() << endl;
  cout << "m.lpNorm<Infinity>() = " << m.lpNorm<Infinity>() << endl;
}

输出结果：

v.squaredNorm() = 5
v.norm() = 2.23607
v.lpNorm<1>() = 3
v.lpNorm<Infinity>() = 2

m.squaredNorm() = 30
m.norm() = 5.47723
m.lpNorm<1>() = 10
m.lpNorm<Infinity>() = 4

四. 布尔约定

all()：所有元素均满足条件，返回值true。
any()：任何元素满足条件，返回true。
count()：满足条件的元素的个数。

#include <Eigen/Dense>
#include <iostream>
using namespace std;
using namespace Eigen;
int main()
{
  ArrayXXf a(2,2);
  a << 1,2,
       3,4;
 
  cout << "(a > 0).all()   = " << (a > 0).all() << endl;
  cout << "(a > 0).any()   = " << (a > 0).any() << endl;
  cout << "(a > 0).count() = " << (a > 0).count() << endl;
  cout << endl;
  cout << "(a > 2).all()   = " << (a > 2).all() << endl;
  cout << "(a > 2).any()   = " << (a > 2).any() << endl;
  cout << "(a > 2).count() = " << (a > 2).count() << endl;
}

输出结果：

(a > 0).all()   = 1
(a > 0).any()   = 1
(a > 0).count() = 4

(a > 2).all()   = 0
(a > 2).any()   = 1
(a > 2).count() = 2

五. 访问元素的地址

使用index进行访问。

#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()
{
  Eigen::MatrixXf m(2,2);
  
  m << 1, 2,
       3, 4;
 
  //get location of maximum
  MatrixXf::Index maxRow, maxCol;
  float max = m.maxCoeff(&maxRow, &maxCol);
 
  //get location of minimum
  MatrixXf::Index minRow, minCol;
  float min = m.minCoeff(&minRow, &minCol);
 
  cout << "Max: " << max <<  ", at: " <<
     maxRow << "," << maxCol << endl;
  cout << "Min: " << min << ", at: " <<
     minRow << "," << minCol << endl;
}

输出结果：

Max: 4, at: 1,1
Min: 1, at: 0,0

六. 对行列访问

之前我们讲过matrix.row(i),和matrix.col(j)进行单独的一行或者一列进行操作。 rowwise() ，colwise()则可以用于每一或者每一列进行访问。

与最大最小操作一起使用

#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
int main()
{
  Eigen::MatrixXf mat(2,4);
  mat << 1, 2, 6, 9,
         		3, 1, 7, 2;
  
  std::cout << "Column's maximum: " << std::endl
   << mat.colwise().maxCoeff() << std::endl;
}

输出结果：

Column's maximum: 
3 2 7 9

与sum等操作一起使用

#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
int main()
{
  MatrixXf mat(2,4);
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;
  
  MatrixXf::Index   maxIndex;
  float maxNorm = mat.colwise().sum().maxCoeff(&maxIndex);//注意维度，这里只有一个索引
  
  std::cout << "Maximum sum at position " << maxIndex << std::endl;
 
  std::cout << "The corresponding vector is: " << std::endl;
  std::cout << mat.col( maxIndex ) << std::endl;
  std::cout << "And its sum is is: " << maxNorm << std::endl;

Maximum sum at position 2
The corresponding vector is: 
6
7
And its sum is is: 13

与广播一起使用

#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
int main()
{
  Eigen::MatrixXf mat(2,4);
  Eigen::VectorXf v(4);
  
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;
         
  v << 0,1,2,3;
       
  //add v to each row of m
  mat.rowwise() += v.transpose();//广播：逐项相加
  
  std::cout << "Broadcasting result: " << std::endl;
  std::cout << mat << std::endl;
}

输出结果：

Broadcasting result: 
 1  3  8 12
 3  2  9  5

混合使用

#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()
{
  Eigen::MatrixXf m(2,4);
  Eigen::VectorXf v(2);
  
  m << 1, 23, 6, 9,
       3, 11, 7, 2;
       
  v << 2,
       3;
 
  MatrixXf::Index index;
  // find nearest neighbour
  (m.colwise() - v).colwise().squaredNorm().minCoeff(&index);
 
  cout << "Nearest neighbour is column " << index << ":" << endl;
  cout << m.col(index) << endl;
}

输出结果：

Nearest neighbour is column 0:
1
3

七. 别名问题

7.1 别名问题的出现

别名(alias[ˈeɪliəs])表示：同一个矩阵（数组或者向量）同时出现在赋值语句左边和右边。有些别名是没问题的，有一些别名问题会出现问题。比如：

mat = 2*mat;//合理
mat = mat.transpose();//不合理；

下面展示一个典型的别名问题：

MatrixXi mat(3,3); 
mat << 1, 2, 3,   4, 5, 6,   7, 8, 9;
cout << "Here is the matrix mat:\n" << mat << endl;
 
// This assignment shows the aliasing problem
mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);
cout << "After the assignment, mat = \n" << mat << endl;

输出结果：

Here is the matrix mat:
1 2 3
4 5 6
7 8 9
After the assignment, mat = 
1 2 3
4 1 2
7 4 1

赋值语句中把左上角的矩阵赋值给右下角。但是有一个问题，mat(1,1)同时出现在左上角和右下角。这个执行过程为：
mat(1,1) = mat(0,0);
mat(1,2) = mat(0,1);
mat(2,1) = mat(1,0);
mat(2,2) = mat(1,1);
所以导致执行第一步的时候mat(1,1)=1，继续执行第四步就不再是自己想要的结果了。

7.2 别名问题的解决

既然知道别名问题是怎么出现的，那我们就知道怎么解决了：先把右边的矩阵当做一个临时矩阵，然后再赋值给左值。eval()就是实现这一功能的函数。

例如上一问题的解决方案：

MatrixXi mat(3,3); 
mat << 1, 2, 3,   4, 5, 6,   7, 8, 9;
cout << "Here is the matrix mat:\n" << mat << endl;
 
// The eval() solves the aliasing problem
mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2).eval();//实现自己想要的功能。
cout << "After the assignment, mat = \n" << mat << endl;

输出：

Here is the matrix mat:
1 2 3
4 5 6
7 8 9
After the assignment, mat = 
1 2 3
4 1 2
7 4 5

转置问题有两种解决方案：

a = a.transpose().eval();
a.transposeInPlace();//Eigen提供的更好的解决方案

提供了In-Place的函数：

Original function	In-place function
MatrixBase::adjoint()	MatrixBase::adjointInPlace()
DenseBase::reverse()	DenseBase::reverseInPlace()
LDLT::solve()	LDLT::solveInPlace()
LLT::solve()	LLT::solveInPlace()
TriangularView::solve()	TriangularView::solveInPlace()
DenseBase::transpose()	DenseBase::transposeInPlace()