Eigen学习教程(一)
Eigen学习教程(一)
1.Eigen介绍
Eigen是一个用于线性代数的C++模板库：矩阵、向量、数值求解器和相关算法。他可以很方便的植入到C++的程序，然后在矩阵计算中有着很广泛的应用。
本教程按照官方教程进行组织，主要内容如下：
- 稠密矩阵与数组操作
- 稠密线性问题与分解
- 稀疏线性代数
- 几何学
2. 稠密矩阵与数组操作
2.1 矩阵类
Eigen中的matrix和vector底层采用Matrix<>模板类表示，因为向量是特殊的矩阵，例如行向量是一行多列，列向量是一列多行。
Matrix<>模板类定义：
Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>
typedef Matrix<double, Dynamic, Dynamic> MatrixXd;
typedef Matrix<int, Dynamic, 1> VectorXi;
Eigen中的MatrixXt基本上是动态内存分配和初始化分开
2.1.1 矩阵类的初始化
固定大小的向量提供了初始化操作
Vector2d a(5.0, 6.0);
Vector3d b(5.0, 6.0, 7.0);
Vector4d c(5.0, 6.0, 7.0, 8.0);
获取元素：
所有元素：m(index)的方式可以获取数据，但是和存储方式有关，
eigen默认采用列优先存储的方法, 但是可以通过设置row-major来实现行优先存储, 注意，存储顺序与逗号初始化是无关系的
[]运算符只能获取向量元素，因为C++[]只支持一个参数；
MatrixXd m(2, 2);
以上程序定义了矩阵大小，但未初始化矩阵,该方式在堆上分配内存。
m(0, 0) = 3;
m(1, 0) = 2.5;
m(0, 1) = -1;
m(1, 1) = m(1, 0) + m(0, 1);
std::cout << m << std::endl;
上述形式通过逐个赋值的形式初始化了定义的m矩阵。
矩阵和矢量系数可以使用所谓的逗号初始化语法方便地设置。现在，只需要知道以下示例即可:
Matrix3f m;
m << 1, 2, 3,
4, 5, 6,
7, 8, 9;
std::cout << m;
2.1.2 矩阵大小的返回与操作
矩阵当前大小可以通过rows()，cols()和size()获取。
这些方法分别返回行数，列数和系数数。调整动态大小矩阵的大小是通过resize()方法完成的。
如果实际矩阵大小不变，则resize()方法为空操作。
否则具有破坏性：系数的值可能会更改。如果你想调整大小，它不改变系数，使用conservativeResize()。
MatrixXd m(2, 5);
m.resize(4, 3);
std::cout << "The matrix m is of size "
<< m.rows() << "x" << m.cols() << std::endl;
std::cout << "It has " << m.size() << " coefficients" << std::endl;
VectorXd v(2);
v.resize(5);
std::cout << "The vector v is of size " << v.size() << std::endl;
std::cout << "As a matrix, v is of size "
<< v.rows() << "x" << v.cols() << std::endl;
为了实现API统一性，所有这些方法仍可用于固定大小的矩阵。当然，您实际上不能调整固定大小的矩阵的大小。尝试将固定大小更改为实际不同的值将触发断言失败。
但是以下代码是合法的:
Matrix4d m;
m.resize(4, 4); // no operation
std::cout << "The matrix m is of size "
<< m.rows() << "x" << m.cols() << std::endl;
固定大小的矩阵or向量在栈上分配内存，因为它在编译时期就可以确定大小
Matrix4f mymatrix ≈ float mymatrix[16];
MatrixX表示运行时才确定矩阵的大小，因为它在堆上分配内存
MatrixXf mymatrix(rows,columns) ≈ float *mymatrix = new float[rows*columns];
示例代码如下:
MatrixXd m = MatrixXd::Random(3, 3); // random values between -1 and 1
m = (m + MatrixXd::Constant(3, 3, 1.2)) * 50; // MatrixXd::Constant(3, 3, 1.2) represents a 3-by-3 matrix expression having all coefficients equal to 1.2
cout << "m =" << endl
<< m << endl;
VectorXd v(3); // 尚未初始化
v << 1, 2, 3; // uses the so-called comma-initializer
cout << "m * v =" << endl
<< m * v << endl;
编译时可确定矩阵尺寸，在栈上分配内存
Matrix3d m = Matrix3d::Random(); //a fixed size
m = (m + Matrix3d::Constant(1.2)) * 50;
cout << "m =" << endl
<< m << endl;
Vector3d v(1, 2, 3);
cout << "m * v =" << endl
<< m * v << endl;
operator= 将矩阵复制到另一个矩阵中的操作。Eigen自动调整左侧矩阵的大小，使其与右侧大小的矩阵大小匹配。例如：
MatrixXf a(2, 2);
std::cout << "a is of size " << a.rows() << "x" << a.cols() << std::endl;
MatrixXf b(3, 3);
a = b;
std::cout << "a is now of size " << a.rows() << "x" << a.cols() << std::endl;
2.1.3 固定尺寸与动态尺寸
什么时候应该使用固定尺寸(例如Matrix4f)，什么时候应该使用动态尺寸(例如MatrixXf)?
简单的答案是: 在可能的地方使用固定尺寸来显示非常小的尺寸，在需要的地方使用动态尺寸来显示较大的尺寸。
- 对于小尺寸，尤其是对于小于（大约）16的尺寸，使用固定尺寸对性能有极大的好处，因为它使Eigen避免了动态内存分配并展开了循环。
- 在内部，固定大小的本征矩阵只是一个简单的数组，即
Matrix4f mymatrix;真的等于只是在做float[16];
因此这确实具有零运行时间成本。
相比之下，动态大小矩阵的数组始终分配在堆上，因此MatrixXf mymatrix（行，列）;等于做 float * mymatrix = new [行*列];，除此之外，MatrixXf对象还将其行数和列数存储为成员变量。
当然，使用固定大小的限制是，只有当您在编译时知道大小时，才有可能这样做。
同样，对于足够大的尺寸(例如，对于大于(大约)32的尺寸)，使用固定尺寸的性能优势变得可以忽略不计。
更糟糕的是，尝试使用函数内部的固定大小创建非常大的矩阵可能会导致堆栈溢出，因为Eigen会尝试自动将数组分配为局部变量，而这通常是在堆栈上完成的。
最后，视情况而定，当使用动态尺寸时，Eigen还可尝试进行矢量化(使用SIMD指令)，请参见参考矢量化。
2.1.4 可选模板参数
我们在页面开始时提到Matrix类采用六个模板参数，但到目前为止，我们仅讨论了前三个。其余三个参数是可选的。这是模板参数的完整列表：
Matrix<typename Scalar,
int RowsAtCompileTime,
int ColsAtCompileTime,
int Options = 0,
int MaxRowsAtCompileTime = RowsAtCompileTime,
int MaxColsAtCompileTime = ColsAtCompileTime>
- Options是位字段。
- 在这里，我们只讨论一点:
RowMajor。它指定这种类型的矩阵使用行优先存储顺序。默认情况下，存储顺序为“按列的顺序存储”。
请参阅有关存储顺序的页面。例如，此类型表示行优先存储的3x3矩阵:
Matrix<float，3、3，RowMajor>
MaxRowsAtCompileTime和MaxColsAtCompileTime在您希望指定时很有用，即使在编译时不知道矩阵的确切大小，在编译时也知道固定的上限。
您可能要这样做的最大原因是避免动态内存分配。例如，以下矩阵类型使用12个浮点数的普通数组，而不分配动态内存：
Matrix<float，Dynamic，Dynamic，0、3、4>
Eigen定义了以下Matrix typedef：
MatrixNt for Matrix<type, N, N>. For example, MatrixXi for Matrix<int, Dynamic, Dynamic>.
VectorNt for Matrix<type, N, 1>. For example, Vector2f for Matrix<float, 2, 1>.
RowVectorNt for Matrix<type, 1, N>. For example, RowVector3d for Matrix<double, 1, 3>.
- N可以是任何一个2，3，4，或X（意思Dynamic）
- t可以是i（表示int），f（表示float），d（表示double），cf（表示complex ）或cd（表示complex ）的任何一种
- 仅针对这五个类型定义typedef的事实并不意味着它们是唯一受支持的标量类型。例如，支持所有标准整数类型，请参阅标量类型
2.2 矩阵和矢量运算
本部分介绍了矩阵在线性代数中的基本运算，相信学过线性代数的同学理解起来十分容易。
需要注意的是：
- 对于a = a.transpose()操作，应该使用a.transposeInPlace(), 见TranspositionAndConjugation()转置
- 向量叉积只适用于大小为3的向量，点积适用于任意向量, 见DotProductAndCrossProduct()
- 理解矩阵Reduction操作的意义
2.2.1 矩阵的加减法
binary operator + as in a+b
binary operator - as in a-b
unary operator - as in -a
compound operator += as in a+=b
compound operator -= as in a-=b
示例代码:
Matrix2d a;
a << 1, 2,
3, 4;
MatrixXd b(2, 2);
b << 2, 3,
1, 4;
std::cout << "a + b =\n"
<< a + b << std::endl;
std::cout << "a - b =\n"
<< a - b << std::endl;
std::cout << "Doing a += b;" << std::endl;
a += b;
std::cout << "Now a =\n"
<< a << std::endl;
Vector3d v(1, 2, 3);
Vector3d w(1, 0, 0);
std::cout << "-v + w - v =\n"
<< -v + w - v << std::endl;
//Output is:
// a + b =
// 3 5
// 4 8
// a - b =
// -1 -1
// 2 0
// Doing a += b;
// Now a =
// 3 5
// 4 8
// -v + w - v =
// -1
// -4
// -6
2.2.2 矩阵的乘除法
binary operator * as in matrixscalar
binary operator * as in scalarmatrix
binary operator / as in matrix/scalar
compound operator = as in matrix=scalar
compound operator /= as in matrix/=scalar
示例代码:
Matrix2d a;
a << 1, 2,
3, 4;
Vector3d v(1, 2, 3);
std::cout << "a * 2.5 =\n"
<< a * 2.5 << std::endl;
std::cout << "0.1 * v =\n"
<< 0.1 * v << std::endl;
std::cout << "Doing v *= 2;" << std::endl;
v *= 2;
std::cout << "Now v =\n"
<< v << std::endl;
//Output is:
// a * 2.5 =
// 2.5 5
// 7.5 10
// 0.1 * v =
// 0.1
// 0.2
// 0.3
// Doing v *= 2;
// Now v =
// 2
// 4
// 6
在Eigen中，诸如算术运算符（例如）operator+自己并不执行任何计算,它们仅返回描述要执行的计算的“表达式对象”。当计算整个表达式时，实际的计算将在稍后进行，通常在中operator=。尽管这听起来很沉重，但是任何现代的优化编译器都可以优化该抽象，从而获得完美优化的代码。例如，当您这样做时：
VectorXf a(50), b(50), c(50), d(50);
...
a = 3*b + 4*c + 5*d;
Eigen将其编译为一个for循环，因此数组仅被遍历一次。简化(例如忽略SIMD优化)，此循环如下所示：
for（int i = 0; i <50; ++ i）
a [i] = 3 * b [i] + 4 * c [i] + 5 * d [i]；
因此，您不必担心Eigen使用相对较大的算术表达式：它只会为Eigen提供更多优化机会。
2.2.3 矩阵的转置矩阵、共轭矩阵、共轭转置矩阵
矩阵或向量的转置 a T a ^ T aT，共轭 a ˉ \bar{a} aˉ和伴随(即共轭转置)分别通过成员函数transpose()，conjugate()和adjoint()获得。
MatrixXcf a = MatrixXcf::Random(2, 2); //MatrixXcf 为复数矩阵
cout << "Here is the matrix a\n"
<< a << endl;
cout << "Here is the matrix a^T\n"
<< a.transpose() << endl;
cout << "Here is the conjugate of a\n"
<< a.conjugate() << endl;
cout << "Here is the matrix a^*\n"
<< a.adjoint() << endl;
//Output is:
// Here is the matrix a
// (-0.211,0.68) (-0.605,0.823)
// (0.597,0.566) (0.536,-0.33)
// Here is the matrix a^T
// (-0.211,0.68) (0.597,0.566)
// (-0.605,0.823) (0.536,-0.33)
// Here is the conjugate of a
// (-0.211,-0.68) (-0.605,-0.823)
// (0.597,-0.566) (0.536,0.33)
// Here is the matrix a^*
// (-0.211,-0.68) (0.597,-0.566)
// (-0.605,-0.823) (0.536,0.33)
注意：对于一个矩阵自身的转置，应该使用.transposeInPlace()方法
Matrix2i a;
a << 1, 2, 3, 4;
cout << "Here is the matrix a:\n"
<< a << endl;
// ～～～～～～不应该这样～～～～～～～～
// a = a.transpose(); // !!! do NOT do this !!!
// cout << "and the result of the aliasing effect:\n"
// << a << endl;
// 应该这样～～～～
a.transposeInPlace();
cout << "and after being transposed:\n"
<< a << endl;
//Output is:
// Here is the initial matrix a:
// 1 2 3
// 4 5 6
// and after being transposed:
// 1 4
// 2 5
// 3 6
2.2.4 矩阵矩阵和矩阵向量乘法
Matrix2d mat;
mat << 1, 2,
3, 4;
Vector2d u(-1, 1), v(2, 0);
std::cout << "Here is mat*mat:\n"
<< mat * mat << std::endl;
std::cout << "Here is mat*u:\n"
<< mat * u << std::endl;
std::cout << "Here is u^T*mat:\n"
<< u.transpose() * mat << std::endl;
std::cout << "Here is u^T*v:\n"
<< u.transpose() * v << std::endl;
std::cout << "Here is u*v^T:\n"
<< u * v.transpose() << std::endl;
std::cout << "Let's multiply mat by itself" << std::endl;
mat = mat * mat;
std::cout << "Now mat is mat:\n"
<< mat << std::endl;
//Output is:
// Here is mat*mat:
// 7 10
// 15 22
// Here is mat*u:
// 1
// 1
// Here is u^T*mat:
// 2 2
// Here is u^T*v:
// -2
// Here is u*v^T:
// -2 -0
// 2 0
// Let's multiply mat by itself
// Now mat is mat:
// 7 10
// 15 22
2.2.5 点积与叉积
// 叉积只适用于大小为3的向量，点积适用于任意向量
Vector3d v(1, 2, 3);
Vector3d w(0, 1, 2);
cout << "Dot product: " << v.dot(w) << endl;
double dp = v.adjoint() * w; // automatic conversion of the inner product to a scalar
cout << "Dot product via a matrix product: " << dp << endl;
cout << "Cross product:\n"
<< v.cross(w) << endl;
// Output is:
// Dot product: 8
// Dot product via a matrix product: 8
// Cross product:
// 1
// -2
// 1
2.2.6 基本算术归约运算
// Eigen also provides some reduction operations to reduce a given matrix or vector to a single value
// such as the sum (computed by sum()), product (prod()), or the maximum (maxCoeff()) and minimum (minCoeff()) of all its coefficients.
Eigen::Matrix2d mat;
mat << 1, 2,
3, 4;
//元素和，元素乘积，元素均值，最小系数，最大系数，踪
cout << "Here is mat.sum(): " << mat.sum() << endl;
cout << "Here is mat.prod(): " << mat.prod() << endl;
cout << "Here is mat.mean(): " << mat.mean() << endl;
cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl;
cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl;
cout << "Here is mat.trace(): " << mat.trace() << endl;
// 可以返回元素位置
Matrix3f m = Matrix3f::Random();
std::ptrdiff_t i, j; // std::ptrdiff_t 是二个指针相减结果的有符号整数类型
float minOfM = m.minCoeff(&i, &j);
cout << "Here is the matrix m:\n"
<< m << endl;
cout << "Its minimum coefficient (" << minOfM
<< ") is at position (" << i << "," << j << ")\n\n";
RowVector4i v = RowVector4i::Random();
int maxOfV = v.maxCoeff(&i);
cout << "Here is the vector v: " << v << endl;
cout << "Its maximum coefficient (" << maxOfV
<< ") is at position " << i << endl;
// Output is:
// Here is mat.sum(): 10
// Here is mat.prod(): 24
// Here is mat.mean(): 2.5
// Here is mat.minCoeff(): 1
// Here is mat.maxCoeff(): 4
// Here is mat.trace(): 5
// Here is the matrix m:
// -0.444451 0.257742 0.904459
// 0.10794 -0.270431 0.83239
// -0.0452059 0.0268018 0.271423
// Its minimum coefficient (-0.444451) is at position (0,0)
2.3 数组和系数运算
与用于线性代数的Matrix类相反，Array类提供了通用数组。此外，Array类提供了一种执行逐系数运算的简便方法，该运算可能没有线性代数含义，例如将常数添加到数组中的每个系数或按系数乘两个数组。
Array<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>
Eigen还为某些常见情况提供了typedef，其方式类似于Matrix typedef，但有一些细微的差异,因为单词“ array”既用于一维数组，也用于二维数组。
- 一维数组采用ArrayNt形式的typedef代表，其中N和t是大小和标量类型;
- 二维数组采用ArrayNNt形式的typedef。
Type Typedef
Array<float,Dynamic,1> ArrayXf
Array<float,3,1> Array3f
Array<double,Dynamic,Dynamic> ArrayXXd
Array<double,3,3>
2.3.1 在数组中访问值
ArrayXXf m(2, 2); // 二维动态float类型数组
// assign some values coefficient by coefficient
m(0, 0) = 1.0;
m(0, 1) = 2.0;
m(1, 0) = 3.0;
m(1, 1) = m(0, 1) + m(1, 0);
// print values to standard output
cout << m << endl
<< endl;
// using the comma-initializer is also allowed
m << 1.0, 2.0,
3.0, 4.0;
// print values to standard output
cout << m << endl;
// Output is:
// 1 2
// 3 5
// 1 2
// 3 4
2.3.2 数组的加减法
两个数组的加减法与矩阵相同。 如果两个数组的大小相同，并且该加法或减法是按系数进行的，则此操作有效。数组还支持以下形式的表达式，该表达式array + scalar将标量添加到数组中的每个系数。这提供了不能直接用于Matrix对象的功能。
ArrayXXf a(3, 3);
ArrayXXf b(3, 3);
a << 1, 2, 3,
4, 5, 6,
7, 8, 9;
b << 1, 2, 3,
1, 2, 3,
1, 2, 3;
// Adding two arrays
cout << "a + b = " << endl
<< a + b << endl
<< endl;
// Subtracting a scalar from an array
cout << "a - 2 = " << endl
<< a - 2 << endl;
// Output is:
// a + b =
// 2 4 6
// 5 7 9
// 8 10 12
// a - 2 =
// -1 0 1
// 2 3 4
// 5 6 7
2.3.3 数组乘法
您可以将一个数组乘以标量，这与矩阵的工作方式相同。数组与矩阵根本不同的地方是将两个矩阵相乘。矩阵将乘法解释为矩阵乘积，而数组将乘法解释为按系数乘积。因此，当且仅当两个数组具有相同的维数时，它们才能相乘。
ArrayXXf a(2, 2);
ArrayXXf b(2, 2);
a << 1, 2,
3, 4;
b << 5, 6,
7, 8;
cout << "a * b = " << endl
<< a * b << endl;
// Output is:
// a * b =
// 5 12
// 21 32
2.3.4 其他系数操作
数组类定义除上述加法，减法和乘法运算符其他系数为单位的运算。例如，.abs()方法获取每个系数的绝对值，而.sqrt（）计算系数的平方根。如果您有两个大小相同的数组，则可以调用.min()来构造其系数是两个给定数组中对应系数的最小值的数组。注意:.array()和.matrix() 转换为相同的维数，但是不同的对象具有不同的方法。
ArrayXf a = ArrayXf::Random(5);
a *= 2;
cout << "a =" << endl
<< a << endl;
cout << "a.abs() =" << endl
<< a.abs() << endl;
cout << "a.abs().sqrt() =" << endl
<< a.abs().sqrt() << endl;
cout << "a.min(a.abs().sqrt()) =" << endl
<< a.min(a.abs().sqrt()) << endl;
// Output is:
// a =
// 1.36
// -0.422
// 1.13
// 1.19
// 1.65
// a.abs() =
// 1.36
// 0.422
// 1.13
// 1.19
// 1.65
// a.abs().sqrt() =
// 1.17
// 0.65
// 1.06
// 1.09
// 1.28
// a.min(a.abs().sqrt()) =
// 1.17
// -0.422
// 1.06
// 1.09
// 1.28
2.3.5 数组和矩阵表达式之间的转换
什么时候应该使用Matrix类的对象，什么时候应该使用Array类的对象？
您不能对数组应用矩阵运算，也不能对矩阵应用数组运算。因此，如果您需要进行线性代数运算（例如矩阵乘法），则应使用矩阵。如果需要进行系数运算，则应使用数组。但是，有时并不是那么简单，但是您需要同时使用Matrix和Array操作。在这种情况下，您需要将矩阵转换为数组或反向转换。无论选择将对象声明为数组还是矩阵，都可以访问所有操作。
矩阵表达式具有.array()方法，可以将它们“转换”为数组表达式，因此可以轻松地应用按系数进行运算。相反，数组表达式具有.matrix()方法。与所有Eigen表达式抽象一样，这没有任何运行时开销（只要您让编译器进行优化。.array()和.matrix()可被用作右值和作为左值。
Eigen禁止在表达式中混合矩阵和数组。
例如，您不能直接矩阵和数组相加。运算符+的操作数要么都是矩阵，要么都是数组。但是，使用.array（）和.matrix（）可以轻松地将其转换为另一种。
～～～～注意～～～～～
此规则的例外是赋值运算符=：允许将矩阵表达式分配给数组变量，或将数组表达式分配给矩阵变量。
下面的示例演示如何通过使用.array()方法对Matrix对象使用数组操作。
例如，语句需要两个矩阵和，他们两个转换为阵列，用来将它们相乘系数明智并将结果指定给矩阵变量（这是合法的，因为本征允许分配数组表达式到矩阵的变量）。result = m.array() * n.array()mnresult
实际上，这种使用情况非常普遍，以至于Eigen为矩阵提供了const .cwiseProduct()方法来计算系数乘积。
MatrixXf m(2, 2);
MatrixXf n(2, 2);
MatrixXf result(2, 2);
m << 1, 2,
3, 4;
n << 5, 6,
7, 8;
result = m * n;
cout << "-- Matrix m*n: --" << endl
<< result << endl
<< endl;
result = m.array() * n.array();
cout << "-- Array m*n: --" << endl
<< result << endl
<< endl;
result = m.cwiseProduct(n);
cout << "-- With cwiseProduct: --" << endl
<< result << endl
<< endl;
result = m.array() + 4;
cout << "-- Array m + 4: --" << endl
<< result << endl
<< endl;
// Output is:
// -- Matrix m*n: --
// 19 22
// 43 50
// -- Array m*n: --
// 5 12
// 21 32
// -- With cwiseProduct: --
// 5 12
// 21 32
// -- Array m + 4: --
// 5 6
// 7 8
2.4 块动作
块是矩阵或数组的某一部分。块表达式既可以用作右值，也可以用作左值。
Block operation
Block of size (p,q), starting at (i,j)
起始位置+行列数matrix.block(i,j,p,q);
模板参数为行列数，函数参数为起始位置
matrix.block<p,q>(i,j);
2.4.1 块操作
两种版本均可用于固定大小和动态大小的矩阵和数组。这两个表达式在语义上是等效的。唯一的区别是，如果块大小较小，则固定大小版本通常会为您提供更快的代码，但是需要在编译时知道此大小。
Eigen::MatrixXf m(4, 4);
m << 1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 14, 15, 16;
cout << "Block in the middle" << endl;
cout << m.block<2, 2>(1, 1) << endl
<< endl;
for (int i = 1; i <= 3; ++i)
{
cout << "Block of size " << i << "x" << i << endl;
cout << m.block(0, 0, i, i) << endl
<< endl;
}
// Output is:
// Block in the middle
// 6 7
// 10 11
// Block of size 1x1
// 1
// Block of size 2x2
// 1 2
// 5 6
// Block of size 3x3
// 1 2 3
// 5 6 7
// 9 10 11
2.4.2 block作为左值
Array22f m;
m << 1, 2,
3, 4;
Array44f a = Array44f::Constant(0.6);
cout << "Here is the array a:" << endl
<< a << endl
<< endl;
a.block<2, 2>(1, 1) = m; //可作为左值
cout << "Here is now a with m copied into its central 2x2 block:" << endl
<< a << endl
<< endl;
a.block(0, 0, 2, 3) = a.block(2, 1, 2, 3);
cout << "Here is now a with bottom-right 2x3 block copied into top-left 2x2 block:" << endl
<< a << endl
<< endl;
// Output is:
// ere is the array a:
// 0.6 0.6 0.6 0.6
// 0.6 0.6 0.6 0.6
// 0.6 0.6 0.6 0.6
// 0.6 0.6 0.6 0.6
// Here is now a with m copied into its central 2x2 block:
// 0.6 0.6 0.6 0.6
// 0.6 1 2 0.6
// 0.6 3 4 0.6
// 0.6 0.6 0.6 0.6
// Here is now a with bottom-right 2x3 block copied into top-left 2x2 block:
// 3 4 0.6 0.6
// 0.6 0.6 0.6 0.6
// 0.6 3 4 0.6
// 0.6 0.6 0.6 0.6
2.4.3 单行单列
单独的列和行是块的特殊情况。Eigen提供了可以轻松解决它们的方法：.col（）和.row（）。
Eigen::MatrixXf m(3, 3);
m << 1, 2, 3,
4, 5, 6,
7, 8, 9;
cout << "Here is the matrix m:" << endl
<< m << endl;
cout << "2nd Row: " << m.row(1) << endl;
m.col(2) += 3 * m.col(0);
cout << "After adding 3 times the first column into the third column, the matrix m is:\n";
cout << m << endl;
// Output is:
// Here is the matrix m:
// 1 2 3
// 4 5 6
// 7 8 9
// 2nd Row: 4 5 6
// After adding 3 times the first column into the third column, the matrix m is:
// 1 2 6
// 4 5 18
// 7 8 30
2.4.4 一系列矩阵的边角块操作
// Block operation Version constructing a dynamic-size block expression Version constructing a fixed-size block expression
// Top-left p by q block * matrix.topLeftCorner(p,q); matrix.topLeftCorner<p,q>();
// Bottom-left p by q block * matrix.bottomLeftCorner(p,q); matrix.bottomLeftCorner<p,q>();
// Top-right p by q block * matrix.topRightCorner(p,q); matrix.topRightCorner<p,q>();
// Bottom-right p by q block * matrix.bottomRightCorner(p,q); matrix.bottomRightCorner<p,q>();
// Block containing the first q rows * matrix.topRows(q); matrix.topRows<q>();
// Block containing the last q rows * matrix.bottomRows(q); matrix.bottomRows<q>();
// Block containing the first p columns * matrix.leftCols(p); matrix.leftCols<p>();
// Block containing the last q columns * matrix.rightCols(q); matrix.rightCols<q>();
示例代码:
Eigen::Matrix4f m;
m << 1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 14, 15, 16;
cout << "m.leftCols(2) =" << endl
<< m.leftCols(2) << endl
<< endl;
cout << "m.bottomRows<2>() =" << endl
<< m.bottomRows<2>() << endl
<< endl;
m.topLeftCorner(1, 3) = m.bottomRightCorner(3, 1).transpose();
cout << "After assignment, m = " << endl
<< m << endl;
// Output is:
// m.leftCols(2) =
// 1 2
// 5 6
// 9 10
// 13 14
// m.bottomRows<2>() =
// 9 10 11 12
// 13 14 15 16
// After assignment, m =
// 8 12 16 4
// 5 6 7 8
// 9 10 11 12
// 13 14 15 16
2.4.5 向量的块运算
同样，向量也有一系列类似于矩阵的边角块操作
Block operation Version constructing adynamic-size block expression Version constructing afixed-size block expression
Block containing the first n elements * vector.head(n); vector.head<n>();
Block containing the last n elements * vector.tail(n); vector.tail<n>();
Block containing n elements, starting at position i * vector.segment(i,n); vector.segment<n>(i);
示例代码:
Eigen::ArrayXf v(6);
v << 1, 2, 3, 4, 5, 6;
cout << "v.head(3) =" << endl
<< v.head(3) << endl
<< endl;
cout << "v.tail<3>() = " << endl
<< v.tail<3>() << endl
<< endl;
v.segment(1, 4) *= 2;
cout << "after 'v.segment(1,4) *= 2', v =" << endl
<< v << endl;
// Output is:
// v.head(3) =
// 1
// 2
// 3
// v.tail<3>() =
// 4
// 5
// 6
// after 'v.segment(1,4) *= 2', v =
// 1
// 4
// 6
// 8
// 10
// 6
2.5 高级初始化方法
本节讨论了几种初始化矩阵的高级方法。
它提供了有关之前介绍的逗号初始化程序的更多详细信息。
它还说明了如何获得特殊矩阵，例如单位矩阵和零矩阵。
2.5.1 逗号初始化值
Eigen提供了一种逗号初始化器语法，该语法使用户可以轻松设置矩阵，向量或数组的所有系数，简单地列出系数，开始在左上角，并从左至右，从顶部向底部移动。 需要预先指定对象的大小。如果列出的系数太少或太多，就会报错。
Matrix3f mat;
mat << 1, 2, 3,
4, 5, 6,
7, 8, 9;
std::cout << mat << endl;
此外，初始化列表的元素本身可以是向量或矩阵。通常的用途是将向量或矩阵连接在一起。例如，这是将两个行向量连接在一起的方法。
注意：
请记住，必须先设置大小，然后才能使用逗号初始化程序。
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;
// 类似分块矩阵的初始化方式
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;
// 对矩阵的某一块赋值
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;
2.5.2 特殊矩阵和数组
模板类Matrix<>和Array<>有静态方法，可以帮助初始化。
有三种变体:
- 第一个变体不带参数，只能用于固定大小的对象。如果要将动态尺寸对象初始化为零，则需要指定尺寸；
- 第二个变体需要一个参数，并且可以用于一维动态尺寸对象；
- 第三个变体需要两个参数，并且可以用于二维对象。
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);
std::cout << a2 << "\n\n";
std::cout << "A two-dimensional dynamic-size array:\n";
ArrayXXf a3 = ArrayXXf::Zero(3, 4);
std::cout << a3 << "\n";
同样，静态方法Constant(value)会将所有系数设置为value。如果需要指定对象的大小，则附加参数放在value参数之前，如MatrixXd::Constant(rows, cols, value)
- Random()用随机系数填充矩阵或数组;
- Identity()获得单位矩阵, 此方法仅适用于Matrix，不适用于Array，因为“单位矩阵”是线性代数概念;
- 该方法LinSpaced（尺寸，低，高）是仅可用于载体和一维数组; 它产生一个指定大小的向量，其系数在low和之间平均间隔high;
- 方法LinSpaced()以下示例说明了该示例，该示例打印一张表格，其中包含以度为单位的角度，以弧度为单位的相应角度以及它们的正弦和余弦值。
示例代码:
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;
Eigen定义了诸如setZero()，MatrixBase :: setIdentity()DenseBase::setLinSpaced()之类的实用程序函数来方便地执行此操作。即可以采用对象的成员函数设置初值。
下面的示例对比了三种构造矩阵的法J=[O I ; I O]
使用静态方法和operator=
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;
// 使用.setXxx()方法
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;
2.5.3 用作临时对象
如上所示，可以在声明时或在赋值运算符的右侧使用静态方法Zero()和Constant()来初始化变量。你可以将这些方法视为返回矩阵或数组。实际上，它们返回所谓的表达式对象，这些表达式对象在需要时求值到矩阵或数组，因此该语法不会产生任何开销。这些表达式也可以用作临时对象。
MatrixXd m = MatrixXd::Random(3, 3);
m = (m + MatrixXd::Constant(3, 3, 1.2)) * 50;
cout << "m =" << endl
<< m << endl;
VectorXd v(3);
v << 1, 2, 3;
cout << "m * v =" << endl
<< m * v << endl;
逗号初始值设定项也可用于构造临时对象。
下面的示例构造了一个大小为2×3的随机矩阵,然后将左边的矩阵乘以[0 1; 1 0]。
MatrixXf mat = MatrixXf::Random(2, 3);
std::cout << mat << std::endl
<< std::endl;
完成临时子矩阵的逗号初始化后，这里需要finished()方法来获取实际的矩阵对象。finished 类似于endl，让它立即生成。
mat = (MatrixXf(2, 2) << 0, 1, 1, 0).finished() * mat;
std::cout << mat << std::endl;
本文主要是作者转载在github的学习教程，链接如下 https://github.com/qixianyu-buaa/EigenChineseDocument，如果涉及到版权问题，请留言作者删除，谢谢！