Eigen学习教程(二)
2.6 约简、访问与广播
2.6.1 约简
在Eigen中，约简是一个采用矩阵或数组并返回单个标量值的函数。最常用的归约方法之一是.sum()，它返回给定矩阵或数组内所有系数的总和。
示例代码如下:
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; // 主对角线上的各个元素的和
2.6.2 范数计算
2范数的平方可以通过squaredNorm()获得。它本身等于矢量的点积，并且等效于其系数的平方绝对值的总和。norm()方法返回squaredNorm()的平方根。
这些运算也可以在矩阵上运算。在这种情况下，n×p矩阵被视为大小(n * p)的向量，因此，例如norm()方法返回 Frobenius或 Hilbert-Schmidt范数。
我们避免谈论
l
2
l^2
l2矩阵的范数，因为那可能意味着不同的事情。
- 如果需要其他按系数分配的l^p范数，请使用
lpNorm <p>(); - 如果需要无穷范数，则模板参数p可以采用特殊值Infinity，这是系数绝对值的最大值。
示例代码:
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;
// 也可自己计算1范数和无穷范数
MatrixXf mat(2, 2);
mat << 1, -2,
-3, 4;
cout << "1-norm(mat) = " << mat.cwiseAbs().colwise().sum().maxCoeff()
<< " == " << mat.colwise().lpNorm<1>().maxCoeff() << endl;
cout << "infty-norm(mat) = " << mat.cwiseAbs().rowwise().sum().maxCoeff()
<< " == " << mat.rowwise().lpNorm<1>().maxCoeff() << endl;
2.6.3 布尔约简
矩阵或者数组皆有与，或，个数的方法
all() 返回true,如果数组或者矩阵中每个元素都符合条件;
any() 返回true,如果数组或者矩阵中有一个元素都符合条件;
count() 返回数组或者矩阵中满足条件的个数。
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;
2.6.4 访问
这是在矩阵和数组的所有元素中，想要获得一个系数在Matrix或Array中的位置时，访问者很有用。最简单的示例是maxCoeff(＆x，＆y)和minCoeff(＆x，＆y)，可用于查找Matrix或Array中最大或最小系数的位置。传递给访问者的参数是指向要存储行和列位置的变量的指针。这些变量应为Index类型。
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;
2.6.5 部分约简
这是在矩阵或数组的列向量和行向量中，记住element-wise是按元素的，那么colwise()或rowwise()表示按列或行的，部分归约是可以在Matrix或Array上按列或按行操作的归约，对每个列或行应用归约运算并返回具有相应值的列或行向量。部分缩减适用于colwise()或rowwise()。一个简单的示例是获取给定矩阵中每一列中元素的最大值，并将结果存储在行向量中: column-wise返回行向量，row-wise返回列向量。
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; // 对于矩阵mat的每一列，取最大系数值
// 也可以对行操作
std::cout << "Row's maximum: " << std::endl
<< mat.rowwise().maxCoeff() << std::endl; // 对于矩阵mat的每一行，取最大系数值
2.6.6 部分约简与其他操作相结合
MatrixXf mat(2, 4);
mat << 1, 2, 6, 9,
3, 1, 7, 2;
MatrixXf::Index maxIndex;
float maxNorm = mat.colwise().sum().maxCoeff(&maxIndex); // 对于矩阵的每一列中的元素求和，结果的最大系数在第2列
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;
// Output is :
//Maximum sum at position 2
// The corresponding vector is:
// 6
// 7
// And its sum is is: 13
2.6.7 广播
广播背后的概念类似于部分归约，区别在于广播构造了一个表达式，其中向量（列或行）通过在一个方向上复制而被解释为矩阵。一个简单的示例是将某个列向量添加到矩阵中的每一列。这可以通过以下方式完成:
Eigen::MatrixXf mat(2, 4);
Eigen::VectorXf v(2);
mat << 1, 2, 6, 9,
3, 1, 7, 2;
v << 0,
1;
//add v to each column of m
//mat.colwise() += v用两种等效的方式解释指令。
//它将向量添加v到矩阵的每一列。或者，可以将其解释为重复向量v四次以形成四乘二矩阵，然后将其加到mat
mat.colwise() += v;
std::cout << "Broadcasting result: " << std::endl;
std::cout << mat << std::endl;
// Output is:
// Broadcasting result:
// 1 2 6 9
// 4 2 8 3
在矩阵上，我们可以执行-=,+=,+,-操作，但是不能进行*=,/=,*,/操作，在数组上我们可执行*=,/=,*,/操作If you want multiply column 0 of a matrix mat with v(0), column 1 with v(1), and so on, then use mat = mat * v.asDiagonal().要逐列或逐行添加的向量必须为Vector类型，并且不能为Matrix。如果不满足，则会出现编译时错误。这也意味着，在使用Matrix进行操作时，广播操作只能应用于Vector类型的对象。这同样适用于数组类VectorXf是ArrayXf。与往常一样，您不应在同一表达式中混合使用数组和矩阵。同样，也可以对行执行此操作:
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
2.6.8 广播与其他操作结合
广播还可以与其他操作(例如矩阵或阵列操作)，归约和部分归约相结合。现在已经介绍了广播，约简和部分约简，我们可以深入研究一个更高级的示例，该示例v在矩阵的列中找到向量的最近邻m。欧几里德距离将在本示例中使用，计算具有部分归约名为squaredNorm()的平方欧几里德距离:
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;
2.7 Map类
本节介绍了如何使用“原始” C / C++数组。这在各种情况下都可能有用，特别是在将其他库中的向量和矩阵“导入”Eigen。有时您可能要在Eigen中使用预定义的数字数组(C++)作为矢量或矩阵(Eigen)。一种选择是复制数据，但最常见的情况是您可能希望将此内存。幸运的是，使用Map类非常容易。Map类实现C++中的数组内存和Eigen对象的交互:
Map< Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> >
请注意，在这种默认情况下，Map仅需要一个模板参数。要构造Map变量，您还需要其他两条信息：指向定义系数数组的内存区域的指针，以及所需的矩阵或矢量形状。(注意区分模板参数和函数形参)例如，要定义一个float在编译时确定大小的矩阵，您可以执行以下操作:
Map <MatrixXf> mf (pf，rows，columns);
其中pf是一个float *指向的存储器阵列。固定大小的整数只读向量可能会声明为Map <const Vector4i> mi (pi);，其中pi是int *。在这种情况下，不必将大小传递给构造函数，因为它已经由Matrix / Array类型指定。请注意，Map没有默认的构造函数。您必须传递一个指针来初始化对象。但是，您可以解决此要求(请参阅更改Map数组)。
Map足够灵活，可以容纳各种不同的数据表示形式。还有其他两个（可选）模板参数:
Map<typename MatrixType, int MapOptions, typename StrideType>
MapOptions指定指针是Aligned还是Unaligned。默认值为Unaligned。
StrideType允许您使用Stride类为内存阵列指定自定义布局。一个示例是指定数据数组以行优先格式进行组织MapConstruct()。
2.7.1 Map构造
int array[8];
for (int i = 0; i < 8; ++i)
array[i] = i;
cout << "Column-major:\n"
<< Map<Matrix<int, 2, 4>>(array) << endl;
cout << "Row-major:\n"
<< Map<Matrix<int, 2, 4, RowMajor>>(array) << endl;
cout << "Row-major using stride:\n"
<< Map<Matrix<int, 2, 4>, Unaligned, Stride<1, 4>>(array) << endl;
//Output is:
//Column-major:
// 0 2 4 6
// 1 3 5 7
// Row-major:
// 0 1 2 3
// 4 5 6 7
// Row-major using stride:
// 0 1 2 3
// 4 5 6 7
2.7.2 使用Map变量
typedef Matrix<float, 1, Dynamic> MatrixType;
typedef Map<MatrixType> MapType;
typedef Map<const MatrixType> MapTypeConst; // a read-only map
const int n_dims = 5;
MatrixType m1(n_dims), m2(n_dims);
m1.setRandom();
m2.setRandom();
float *p = &m2(0); // get the address storing the data for m2
MapType m2map(p, m2.size()); // m2map shares data with m2
MapTypeConst m2mapconst(p, m2.size()); // a read-only accessor for m2
cout << "m1: " << m1 << endl;
cout << "m2: " << m2 << endl;
cout << "Squared euclidean distance: " << (m1 - m2).squaredNorm() << endl;
cout << "Squared euclidean distance, using map: " << (m1 - m2map).squaredNorm() << endl;
m2map(3) = 7; // this will change m2, since they share the same array
cout << "Updated m2: " << m2 << endl;
cout << "m2 coefficient 2, constant accessor: " << m2mapconst(2) << endl;
/* m2mapconst(2) = 5; */ // this yields a compile-time error
//Output is:
// m1: 0.68 -0.211 0.566 0.597 0.823
// m2: -0.605 -0.33 0.536 -0.444 0.108
// Squared euclidean distance: 3.26
// Squared euclidean distance, using map: 3.26
// Updated m2: -0.605 -0.33 0.536 7 0.108
// m2 coefficient 2, constant accessor: 0.536
2.7.3 改变Map的数组
可以使用C ++placement new(位置new，在程序员给定的内存放置元素) 语法更改已声明的Map对象的数组。尽管有出现，但它不会调用内存分配器，因为语法指定了存储结果的位置(具体参考《C++primer》中的内容)。简单的说，位置new只是在指定位置写入内容，并不会像new一样，先在堆上分配内存，然后再依次初始化对象，这也是为什么叫位置new，因为它会按照我们指定的位置构造对象。
int data[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
Map<RowVectorXi> v(data, 4);
cout << "The mapped vector v is: " << v << "\n";
new (&v) Map<RowVectorXi>(data + 4, 5);
cout << "Now v is: " << v << "\n";
//The mapped vector v is: 1 2 3 4
//Now v is: 5 6 7 8 9
2.8 重塑与切片
2.8.1 重塑
整形操作在于修改矩阵的大小，同时保持相同的系数。除了修改输入矩阵本身（这对于编译时大小而言是不可能的）之外，该方法还包括使用Map类在存储上创建不同的视图。这是创建矩阵的一维线性视图的典型示例:
MatrixXf M1(3, 3); // Column-major storage
// 注意：逗号初始化是为了方便我们输入矩阵，但是底层存储是按照列主的顺序存储的
M1 << 1, 2, 3,
4, 5, 6,
7, 8, 9;
Map<RowVectorXf> v1(M1.data(), M1.size());
cout << "v1:" << endl
<< v1 << endl;
Matrix<float, Dynamic, Dynamic, RowMajor> M2(M1);
Map<RowVectorXf> v2(M2.data(), M2.size());
cout << "v2:" << endl
<< v2 << endl;
//Output is:
// v1:
// 1 4 7 2 5 8 3 6 9
// v2:
// 1 2 3 4 5 6 7 8 9
注意输入矩阵的存储顺序如何修改线性视图中系数的顺序。这是另一个将2x6矩阵重塑为6x2矩阵的示例:
MatrixXf M1(2, 6); // Column-major storage
M1 << 1, 2, 3, 4, 5, 6,
7, 8, 9, 10, 11, 12;
Map<MatrixXf> M2(M1.data(), 6, 2);
cout << "M2:" << endl
<< M2 << endl;
//Output is:
// M2:
// 1 4
// 7 10
// 2 5
// 8 11
// 3 6
// 9 12
2.8.2 切片
切片包括获取一组在矩阵内均匀间隔的行，列或元素。再次，Map类可以轻松模仿此功能。例如，可以跳过向量中的每个P元素:
RowVectorXf v = RowVectorXf::LinSpaced(20, 0, 19);
cout << "Input:" << endl
<< v << endl;
Map<RowVectorXf, 0, InnerStride<2>> v2(v.data(), v.size() / 2);
cout << "Even:" << v2 << endl;
//Output is:
// Input:
// 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
// Even: 0 2 4 6 8 10 12 14 16 18
一个人可以根据实际的存储顺序，使用适当的外步幅或内步幅来占据三列中的一列:
MatrixXf M1 = MatrixXf::Random(3, 8);
cout << "Column major input:" << endl
<< M1 << "\n";
Map<MatrixXf, 0, OuterStride<>> M2(M1.data(), M1.rows(), (M1.cols() + 2) / 3, OuterStride<>(M1.outerStride() * 3));
cout << "1 column over 3:" << endl
<< M2 << "\n";
typedef Matrix<float, Dynamic, Dynamic, RowMajor> RowMajorMatrixXf;
RowMajorMatrixXf M3(M1);
cout << "Row major input:" << endl
<< M3 << "\n";
Map<RowMajorMatrixXf, 0, Stride<Dynamic, 3>> M4(M3.data(), M3.rows(), (M3.cols() + 2) / 3,
Stride<Dynamic, 3>(M3.outerStride(), 3));
cout << "1 column over 3:" << endl
<< M4 << "\n";
//Column major input:
// 0.68 0.597 -0.33 0.108 -0.27 0.832 -0.717 -0.514
// -0.211 0.823 0.536 -0.0452 0.0268 0.271 0.214 -0.726
// 0.566 -0.605 -0.444 0.258 0.904 0.435 -0.967 0.608
// 1 column over 3:
// 0.68 0.108 -0.717
// -0.211 -0.0452 0.214
// 0.566 0.258 -0.967
// Row major input:
// 0.68 0.597 -0.33 0.108 -0.27 0.832 -0.717 -0.514
// -0.211 0.823 0.536 -0.0452 0.0268 0.271 0.214 -0.726
// 0.566 -0.605 -0.444 0.258 0.904 0.435 -0.967 0.608
// 1 column over 3:
// 0.68 0.108 -0.717
// -0.211 -0.0452 0.214
// 0.566 0.258 -0.967
2.9 混淆
这一节比较绕:
在Eigen中，混淆(aliasing)是指相同的矩阵（或数组或向量）出现在赋值运算符的左侧和右侧的赋值语句;
例如，A = AB , a = a^Tb A=A*A
产生混淆的原因是Eigen采用惰性求值。
混淆可能是有害的，也可能是无害的,有害的混淆，导致可能不正确的结果，无害的混淆可以产生正确的结果
-
有害混淆
使用.eval()方法，可以解决混淆问题。具体的说，.eval()方法生成临时对象，然后再执行赋值操作
如果Eigen中假定该操作是混淆的，那么它会自动使用.eval()方法，而不需要我们显示调用 -
无害混淆
无害混淆是我们不需要评估，它对系数计算无害。这包括标量乘法和矩阵或数组加法。
将两个矩阵相乘时，Eigen会假定发生混叠(注意，Eigen3.3以后的版本中需要区分目标矩阵大小是否改变了)
如果您知道没有混淆，则可以使用noalias()。
在所有其他情况下，Eigen假定不存在混叠问题，因此如果实际上发生混叠，则会给出错误的结果。
为防止这种情况，您必须使用eval()或xxxInPlace()函数之一。
2.9.1 混淆
在Eigen中，混淆(aliasing)是指相同的矩阵（或数组或向量）出现在赋值运算符的左侧和右侧的赋值语句。例如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;
// Output is:
// 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.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);该赋值具有混淆aliasing:系数mat(1,1)既出现在mat.bottomRightCorner(2,2)分配左侧的块中mat.topLeftCorner(2,2)，又出现在右侧的块中。分配后，右下角的(2,2)项应具有mat(1,1)分配前的值5。但是，输出显示mat(2,2)实际上为1。问题在于Eigen使用了惰性求值。结果类似于
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(2,2)分配了新值mat(1,1)而不是旧值。可以通过调用eval()解决此问题(注:eval()负责生成临时对象而避免混淆)， 尝试缩小矩阵时，混淆也会发生。例如，表达式vec = vec.head(n)和mat = mat.block(i,j,r,c)具有混淆。通常，在编译时无法检测到混淆：如果mat在第一个示例中稍大一点，则块将不会重叠，也不会出现混淆问题。但是Eigen确实会在运行时检测到一些混淆实例。Matrix和向量算术中提到了以下显示别名的示例：
// 注意：这段代码会报错～～～～～～
// 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;
// 输出显示混淆(alising)问题。
// 但是，默认情况下，Eigen使用运行时断言来检测到此情况并退出，并显示如下消息
// void Eigen::DenseBase<Derived>::checkTransposeAliasing(const OtherDerived&) const
// [with OtherDerived = Eigen::Transpose<Eigen::Matrix<int, 2, 2, 0, 2, 2> >, Derived = Eigen::Matrix<int, 2, 2, 0, 2, 2>]:
// Assertion `(!internal::check_transpose_aliasing_selector<Scalar,internal::blas_traits<Derived>::IsTransposed,OtherDerived>::run(internal::extract_data(derived()), other))
// && "aliasing detected during transposition, use transposeInPlace() or evaluate the rhs into a temporary using .eval()"' failed.
用户可以通过定义EIGEN_NO_DEBUG宏来关闭Eigen的运行时断言
2.9.2 解决混淆问题
解决方法：Eigen必须将右侧完全看作一个临时矩阵/数组，然后将其分配给左侧。函数**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;
相同的解决方案也适用于第二示例中，与转置:只需更换线a = a.transpose();用a = a.transpose().eval();但是在这种常见情况下有更好的解决方案。Eigen提供了专用函数transposeInPlace()，该函数通过其转置来替换矩阵。如下所示:
MatrixXf a(2, 3);
a << 1, 2, 3, 4, 5, 6;
cout << "Here is the initial matrix a:\n"
<< a << endl;
a.transposeInPlace();
cout << "and after being transposed:\n"
<< a << endl;
如果xxxInPlace()函数可用，则最好使用它，因为它可以更清楚地指示您正在做什么。
这也可以让Eigen更积极地进行优化。这些是提供的一些xxxInPlace()函数：
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()
在特殊情况下，矩阵或向量使用类似的表达式缩小vec = vec.head(n)，则可以使用conservativeResize()。
2.9.3 混淆和组件操作
如果同一矩阵或数组同时出现在赋值运算符的左侧和右侧，则可能很危险，因此通常有必要显示地评估右侧。但是，应用基于元素的操作(例如矩阵加法，标量乘法和数组乘法)是安全的。以下示例仅具有基于组件的操作。因此，即使相同的矩阵出现在赋值符号的两侧，也不需要eval()。
MatrixXf mat(2, 2);
mat << 1, 2, 4, 7;
cout << "Here is the matrix mat:\n"
<< mat << endl
<< endl;
mat = 2 * mat;
cout << "After 'mat = 2 * mat', mat = \n"
<< mat << endl
<< endl;
mat = mat - MatrixXf::Identity(2, 2);
cout << "After the subtraction, it becomes\n"
<< mat << endl
<< endl;
ArrayXXf arr = mat;
arr = arr.square();
cout << "After squaring, it becomes\n"
<< arr << endl
<< endl;
// Combining all operations in one statement:
mat << 1, 2, 4, 7;
mat = (2 * mat - MatrixXf::Identity(2, 2)).array().square();
cout << "Doing everything at once yields\n"
<< mat << endl
<< endl;
// Output is:
// Here is the matrix mat:
// 1 2
// 4 7
// After 'mat = 2 * mat', mat =
// 2 4
// 8 14
// After the subtraction, it becomes
// 1 4
// 8 13
// After squaring, it becomes
// 1 16
// 64 169
// Doing everything at once yields
// 1 16
// 64 169
通常，如果表达式右侧的(i，j)条目仅取决于左侧矩阵或数组的(i，j)条目,而不依赖于其他任何表达式，则赋值是安全的。在这种情况下，不必显示地评估右侧(.evl())。
2.9.4 混淆和矩阵乘法
在目标矩阵未调整大小的情况下，矩阵乘法是Eigen中唯一假定默认情况下为混淆的。若假定混淆，则会使用eval()生成临时对象,所以是安全的。因此，如果matA是平方矩阵，则该语句matA = matA * matA是安全的。Eigen中的所有其他操作都假定没有混淆问题，这是因为结果被分配给了不同的矩阵，或者因为它是逐个元素的操作。
MatrixXf matA(2, 2);
matA << 2, 0, 0, 2;
matA = matA * matA;
cout << matA << endl;
但是，这是有代价的。执行表达式时matA = matA * matA，Eigen会在计算后的临时矩阵中评估赋值给matA的乘积。虽然可以，但是在将乘积分配给其他矩阵（例如matB = matA * matA）时，Eigen会执行相同的操作。在这种情况下，直接评估matB，而不是先将matA*matA生成临时对象，然后评估临时对象，最后将临时对象赋值给矩阵matB更高效。我们可以使用noalias函数指示没有混淆，如下所示:matB.noalias() = matA * matA。这使Eigen可以matA * matA直接将矩阵乘积在matB中评估。
MatrixXf matA(2, 2), matB(2, 2);
matA << 2, 0, 0, 2;
// Simple but not quite as efficient
matB = matA * matA;
cout << matB << endl;
// More complicated but also more efficient
matB.noalias() = matA * matA;
cout << matB << endl;
// Output is:
// 4 0
// 0 4
// 4 0
// 0 4
当然，不应该在实际上发生混淆时使用noalias()，如果这样做，则可能会得到错误的结果:我的平台上没报错
MatrixXf matA(2, 2);
matA << 2, 0, 0, 2;
matA.noalias() = matA * matA;
cout << matA << endl;
// Output is:
// 4 0
// 0 4
此外，从Eigen 3.3 开始，如果调整了目标矩阵的大小(注意，上文中的操作假定目标矩阵大小不变)，并且未将乘积直接分配给目标，则不假定混淆。因此，以下示例也是错误的:
MatrixXf A(2, 2), B(3, 2);
B << 2, 0, 0, 3, 1, 1;
A << 2, 0, 0, -2;
A = (B * A).cwiseAbs(); // 由于不假定混淆，所以需要我们显示评价
cout << A << endl;
// Output is
// 4 0
// 0 6
// 2 2
对于任何混淆问题，您都可以通过在赋值之前显式评估表达式来解决它:
MatrixXf A(2, 2), B(3, 2);
B << 2, 0, 0, 3, 1, 1;
A << 2, 0, 0, -2;
A = (B * A).eval().cwiseAbs();
cout << A << endl;
// Output is
// 4 0
// 0 6
// 2 2
2.10 存储顺序
矩阵和二维数组有两种不同的存储顺序：列优先和行优先。本页说明了这些存储顺序以及如何指定应使用的存储顺序。矩阵的条目形成一个二维网格。但是，当矩阵存储在存储器中时，必须以某种方式线性排列条目。有两种主要方法可以做到这一点，按行和按列。我们说矩阵是按行优先存储。首先存储整个第一行，然后存储整个第二行，依此类推。另一方面，如果矩阵是逐列存储的，则以主列顺序存储，从整个第一列开始，然后是整个第二列，依此类推，可以通过Options为Matrix或Array指定模板参数来设置矩阵或二维数组的存储顺序。
由于Matrix类解释，Matrix类模板有六个模板参数，其中三个是强制性的（Scalar，RowsAtCompileTime和ColsAtCompileTime）三个是可选的（Options，MaxRowsAtCompileTime和MaxColsAtCompileTime）。如果Options参数设置为RowMajor，则矩阵或数组以行优先顺序存储：如果将其设置为ColMajor，则以列优先顺序存储。如果未指定存储顺序，则Eigen默认将条目存储在column-major中。如果方便的typedef（Matrix3f，ArrayXXd等）也是默认按列存储。
可以将使用一种存储顺序的矩阵和数组分配给使用另一种存储顺序的矩阵和数组，如一个按行存储的矩阵使用按列存储矩阵初始化。Eigen将自动对元素重新排序。更一般而言，按行存储矩阵和按列存储矩阵可以根据需要在表达式中混合使用。
选择哪个存储顺序？
那么，您应该在程序中使用哪个存储顺序？这个问题没有简单的答案。这取决于您的应用程序。请记住以下几点：
您的用户可能希望您使用特定的存储顺序。或者，您可以使用Eigen以外的其他库，并且这些其他库可能需要一定的存储顺序。在这些情况下，在整个程序中使用此存储顺序可能是最简单，最快的。当矩阵以行优先顺序存储时，逐行遍历矩阵的算法会更快，因为数据位置更好。同样，对于主要列矩阵，逐列遍历更快。可能需要尝试一下以找出对您的特定应用程序更快的方法。
Eigen中的默认值是column-major。自然，对Eigen库的大多数开发和测试都是使用列主矩阵完成的。这意味着，即使我们旨在透明地支持列主存储和行主存储顺序，Eigen库也最好与列主存储矩阵配合使用。
Matrix<int, 3, 4, ColMajor> Acolmajor;
Acolmajor << 8, 2, 2, 9,
9, 1, 4, 4,
3, 5, 4, 5;
cout << "The matrix A:" << endl;
cout << Acolmajor << endl
<< endl;
cout << "In memory (column-major):" << endl;
for (int i = 0; i < Acolmajor.size(); i++)
cout << *(Acolmajor.data() + i) << " ";
cout << endl
<< endl;
Matrix<int, 3, 4, RowMajor> Arowmajor = Acolmajor;
cout << "In memory (row-major):" << endl;
for (int i = 0; i < Arowmajor.size(); i++)
cout << *(Arowmajor.data() + i) << " ";
cout << endl;
// Output is
// The matrix A:
// 8 2 2 9
// 9 1 4 4
// 3 5 4 5
// In memory (column-major):
// 8 9 3 2 1 5 2 4 4 9 4 5
// In memory (row-major):
// 8 2 2 9 9 1 4 4 3 5 4 5
2.11 问题总
- 对齐错误
Eigen::internal::matrix_array<T, Size, MatrixOptions, Align>::internal::matrix_array()
[with T = double, int Size = 2, int MatrixOptions = 2, bool Align = true]:
Assertion `(reinterpret_cast<size_t>(array) & (sizemask)) == 0 && "this assertion
is explained here: http://eigen.tuxfamily.org/dox-devel/group__TopicUnalignedArrayAssert.html
READ THIS WEB PAGE !!! ****"' failed.
首先找到程序触发位置，
例如，
$gdb ./my_program
run
…
bt
然后与下面的4种原因对号入座，修改代码
- 四种原因
原因1:结构体中具有Eigen对象成员
请注意，此处Eigen :: Vector2d仅用作示例，更一般而言，所有固定大小的可矢量化Eigen类型都会出现此问题。固定大小的可矢量化Eigen类型是如果它具有固定的大小并且大小是16字节的倍数。
Eigen::Vector2d
Eigen::Vector4d
Eigen::Vector4f
Eigen::Matrix2d
Eigen::Matrix2f
Eigen::Matrix4d
Eigen::Matrix4f
Eigen::Affine3d
Eigen::Affine3f
Eigen::Quaterniond
Eigen::Quaternionf
首先, "固定大小"应该清楚：如果在编译时，Eigen对象的行数和列数是固定的，则其固定大小。因此，例如Matrix3f具有固定大小，但MatrixXf没有(固定大小的对立是动态大小)。固定大小的Eigen对象的系数数组是普通的“静态数组”，不会动态分配。例如，Matrix4f后面的数据只是一个“float array[16]”。固定大小的对象通常很小，这意味着我们要以零的运行时开销（在内存使用和速度方面）来处理它们。现在，矢量化（SSE和AltiVec）都可以处理128位数据包。此外，出于性能原因，这些数据包必须具有128位对齐。因此，事实证明，固定大小的Eigen对象可以向量化的唯一方法是，如果它们的大小是128位或16个字节的倍数。然后，Eigen将为这些对象请求16字节对齐，并且此后将依赖于这些对象进行对齐，因此不会执行运行时检查以进行对齐。
class Foo
{
//...
Eigen::Vector2d v;
//...
};
Foo *foo = new Foo;
Eigen需要Eigen :: Vector2d的数组（2个双精度）的128位对齐。对于GCC，这是通过属性 ( ( aligned(16) ) )完成的。Eigen重载了Eigen::Vector2d的operator new，因此它将始终返回128位对齐的指针。因此，通常情况下，您不必担心任何事情，Eigen会为您处理对齐。除了一种情况。当您具有上述的Foo类，并且如上所述动态分配新的Foo时，则由于Foo没有对齐“ operator new”，因此返回的指针foo不一定是128位对齐的。然后，成员v的alignment属性相对于类的开头foo。如果foo指针未对齐，则foo-> v也不会对齐！解决方案是让Foo类具有一致的Operator new
解决方法：
如果定义的结构具有固定大小的可矢量化Eigen类型的成员，则必须重载其operator new，以便它生成16字节对齐的指针。幸运的是，Eigen为您提供了一个宏EIGEN_MAKE_ALIGNED_OPERATOR_NEW来为您执行此操作。换句话说:您有一个类，该类具有固定大小的可矢量化Eigen对象作为成员，然后动态创建该类的对象。
很简单，您只需将EIGEN_MAKE_ALIGNED_OPERATOR_NEW宏放在您的类的public部分。
class Foo
{
Eigen::Vector2d v;
public : EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};
Foo *foo = new Foo;
该宏使new Foo始终返回对齐的指针。一个 Eigen::Vector2d有两个double类型，一个double为8字节=64位，则一个Eigen::Vector2d为128位。这恰好是SSE数据包的大小，这使得可以使用SSE对该向量执行各种操作。但是SSE指令（至少Eigen使用的是快速指令）需要128位对齐。否则会出现段错误。出于这个原因，Eigen自己通过执行以下两项工作自行要求对Eigen::Vector2d进行128位对齐。
原因2：STL容器或手动内存分配
如果您将STL容器(例如std :: vector，std :: map等)与Eigen对象或包含Eigen对象的类一起使用。
std::vector<Eigen::Matrix2f> my_vector;
struct my_class { ... Eigen::Matrix2f m; ... };
std::map<int, my_class> my_map;
请注意，此处Eigen::Matrix2f仅用作示例，更一般而言对于所有固定大小的可矢量化Eigen类型和具有此类Eigen对象作为member的结构，都会出现此问题。任何类/函数都会绕过new分配内存的运算符，也就是通过执行自定义内存分配，然后调用placement new运算符，也会出现相同的问题。例如，通常就是这种情况，std::make_shared或者std::allocate_shared解决方案是使用对齐的分配器，如STL容器解决方案中所述。
原因3：通过值传递Eigen对象
如果您代码中的某些函数正在通过值传递Eigen对象，例如这样，void func(Eigen::Vector4d v);那么您需要阅读以下单独的页面：将Eigen对象按值传递给函数。请注意，此处Eigen :: Vector4d仅用作示例，更一般而言，所有固定大小的可矢量化Eigen类型都会出现此问题。
原因4:编译器对堆栈对齐做出错误假设（例如Windows上的GCC），这是在Windows上使用GCC（例如MinGW或TDM-GCC）的人们的必读内容。如果在声明这样的局部变量的无辜函数中有此断言失败，请执行以下操作:
void foo()
{
Eigen::Quaternionf q;
//...
}
那么您需要阅读以下单独的页面：编译器对堆栈对齐做出了错误的假设。请注意，此处Eigen :: Quaternionf仅用作示例，更一般而言，所有固定大小的可矢量化Eigen类型都会出现此问题。
3.该断言的一般解释:
固定大小的矢量化Eigen对象必须绝对在16字节对齐的位置创建，否则寻址它们的SIMD指令将崩溃。Eigen通常会为您解决这些对齐问题，方法是在固定大小的可矢量化对象上设置对齐属性，并重载他们的“ operator new”。但是，在某些极端情况下，这些对齐设置会被覆盖：它们是此断言的可能原因。
4.我不在乎最佳矢量化，如何摆脱这些东西？
三种可能性：
使用Matrix，Array，Quaternion等对象的DontAlign选项会给您带来麻烦。这样，Eigen不会尝试对齐它们，因此不会采取任何特殊对齐方式。不利的一面是，您将为它们支付未对齐的加载/存储的成本，但是在现代CPU上，开销为null或边际的。定义EIGEN_DONT_ALIGN_STATICALLY。这将禁用所有16字节（或以上）的静态对齐代码，同时保持16字节（或以上）的堆对齐。这具有通过未对齐的存储区（由EIGEN_UNALIGNED_VECTORIZE控制）对固定大小的对象（如Matrix4d）进行矢量化的效果，同时保持动态大小的对象（如MatrixXd）的矢量化不变。但是请注意，这会破坏ABI与静态对齐方式默认行为的兼容性。 或同时定义EIGEN_DONT_VECTORIZE和EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT。这样可以保留16字节的对齐代码，从而保留ABI兼容性，但完全禁用向量化。如果您想知道为什么定义EIGEN_DONT_VECTORIZE本身并不能禁用16字节对齐和断言，则说明如下:
它不会禁用断言，因为如果不执行矢量化，则正常运行的代码将在启用矢量化时突然崩溃。它不会禁用16字节对齐，因为这将意味着矢量化和非矢量化的代码不相互兼容ABI。即使对于仅开发内部应用程序的人，这种ABI兼容性也非常重要，例如，可能希望在同一应用程序中同时具有矢量化路径和非矢量化路径。
本文主要是作者转载在github的学习教程，链接如下https://github.com/qixianyu-buaa/EigenChineseDocument，如果涉及到版权问题，请留言作者删除，谢谢！
242
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
