关闭

Eigen初学相关介绍

标签: Eigen矩阵处理c++
1086人阅读 评论(0) 收藏 举报
分类:

最近在用Matlab处理图像,现在要做的是将其用C++语言进行翻译,由于要进行大量的矩阵计算,就研究了一下可以进行矩阵计算的开源库,详细的介绍可以参照http://my.oschina.net/cvnote/blog/165340,我从中选择了Eigen进行了一番学习,现在对里面一些基础知识做一下小结。以下内容可以看做它官方在线文档的一个学习笔记,粗略看看还是感觉很强大的,而且由于只包含头文件,方便跨平台使用,打算去使用一下。详细内容可以参照官方文档:http://eigen.tuxfamily.org/index.php?title=Main_Page 。

Dense matrix and array manipulation

The Matrix class

    创建语法:Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime, int Options, int MaxRowsAtCompileTime, int MaxColsAtCompileTime >, 其中Scalar是矩阵类型,有int、float、double、complex float、complex double等。

    RowsAtCompileTime和ColsAtCompileTime是编译时创建数组的大小。Eigen还定义了许多convenience typedefs 比如

Typedef Matrix<float, 4, 4> Matrix4f;

    Options可以定义一些矩阵的属性,不如存储顺序是按列存储还是按行存储;

    MaxRowsAtCompileTime 和MaxColsAtCompileTime 是规定矩阵的最大范围,以防动态分配(下面介绍)时分配过大。

    向量是一种特殊的矩阵,也有很多convenience typedef,如

Typedef Matrix<float, 3, 1> Vector3f;

    这种定义方式得到的是列向量,数据是按列存储,还可以定义行向量:

Typedef Matrix<int, 1, 2> RowVector2i;

    矩阵大小也可以不在编译时设置,而是动态设定,如

Matrix(double, Dynamic, Dynamic)

    向量也是如此:

Matrix(int, Dynamic, 1)

    定义时可以对其进行初始化:

MatrixXf a(10, 15); //定义了一个10*15的矩阵

Vector2d a(5.0, 6.0);

    矩阵变量允许下标访问,不过是采用括号,例如定义一个变量

Matri3dx m = Matrix3d::Random();    //生成一个3*3的随机矩阵。

    访问第二行第三个元素,可以使用m(1, 2)来进行访问。

rows(), cols(), size()函数可以用以返回矩阵的行数、列数和元素个数。

    Resize()函数可以用来改变矩阵或者向量的大小,但对于行数和列数固定的矩阵或向量进行修改会导致错误产生。

    对于编译时未规定数组,如果在定义时两个矩阵长度不同,也可以进行赋值操作,只不过会将被复制的矩阵size改变,如

MatrixXf a(2, 2); MatrixXf b(3, 3);

b = a;

    则此时b是一个2*2的矩阵。

    需要说明的是,如果声明一个矩阵用固定大小的方式(e.g. Matrix4f),效率比用动态动态声明(MatrixXf a(4, 4))高很多。

Matrix and Vector Arithmetic

    定义了矩阵相加、相减、相乘等操作

 

    需要注意的是,在做矩阵计算是,用的是循环语句,所以是串行计算的方式。

    Eigen还提供了一些常用的变换函数:transpose(),转置操作:conjugata(),共轭操作;adjoint(),共轭转置操作;需要注意的是这些是非原位操作,例如a = [1, 2;, 3, 4] 那么a = a.transpose() 得到的a为[1, 2, 2, 4],这样就产生了所粗。原位转置操作可以使用transposeInPlace() 函数。

    矩阵乘法是一个特殊的运算,为了防止上述的错误,矩阵与矩阵相乘都会使用临时变量,例如m = m* m; 在实际编译的时候是tem = m* m; m = tem; 如果在乘法中你可以确定不会有错误的产生,那么可以使用noaliasd()函数来避免产生临时变量,例如:c.noalias() += a * b;

    Eigen还提供了dot()点乘函数和cross()叉乘函数。其中dot可以针对任意两个长度相等的向量,而叉乘只能是两个三维向量,例如Vector3d v(1, 2, 3); Vector3d w(0, 1, 2); 那么v.dot(w) 得到的结果是8(等价于v.adjoint() * w),v.corss(w)得到的结果是(1;-2;1)。

    Eigen提供了一些基础的算法,如sum(),求所有元素和;prod(),求所有元素之积;mean(),求元素平均值;minCoeff(),最小元素;maxCoeff(),最大元素;trace(),求迹等。下面给出一个实例:

    上述例子不仅求出了最小值,也求得了最小值所在的位置,类似于Matlab中find的用法。

The Array class and coefficient-wise operations

    Array也是一个模板类,可以跟Matrix进行对比,声明一个Array对象的方式与Matrix相似:

Array<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

    参数与Matrix一致,这里不多做介绍,详细可看官方文档。Array也定义了许多conventional typedef,如

    Array的赋值、加法、减法也与Matrix类似,需要说明的是Array相乘,这个操作类似于Matlab中的".*",对应元素做计算,所以两个Array相乘,只能是大小相同的Array。

    Array提供了一些操作函数,abs(),求每个元素的绝对值;sqrt(),求每个元素的算术平方根;a.min(b),求a和b中每个位置较小的元素。

    Matrix和Array可以进行互换,利用matrix中的Array()函数和Array中的Matrix()函数,下面给出实例

    需要说明的是,将array表达式赋值给Matrix变量是允许的;cwiseProduct()函数允许Matrix直接进行点对点乘法,而不用转换至Array。

Block Operations

    利用block()函数,可以从Matrix中取出一个小矩阵来进行处理,使用的语法为matrix.block(i, j, p, q); 和 matrix.block<p, q>(i, j)。其中i,j是block左上角元素位于矩阵中的位置,p、q是block的大小,如果是小块的话,那么采用第二种方法,在编译时就固定p、q,速度会比第一种方法快一些。

    需要特别注意的是,block利用非常灵活,原因在于它不仅仅可以是右值表达式,它也可以作为左值对matrix进行修改。同样Array也有block函数,如下面实例

    Block函数可以从矩阵中提取一块,特殊的如提取一整行或者一整列,如想达到如此操作,可以使用row()函数和col()函数,已达到更高效率。这两个函数很类似与Matlab中":"操作符。

    如果想在矩阵中操作角落的block,Eigen还提供了一些特殊的函数:

    最后一列如果在操作大小固定的小block时,在编译时提供size,可以达到更快的效果。同样对于vector,也有类似的操作函数:

Advanced Initialization

    简单初始一个矩阵的方法是利用逗号表达式,初始顺序是从左到右、从上到下。可以一个元素一个元素的初始矩阵,也可以一次用一个矩阵去初始大矩阵的一部分,如

    同样利用上面讲的block去初始矩阵的一部分:

    Eigen提供了一些函数来初始特殊的矩阵,如Zero()用以初始全领的矩阵,其size可以通过参数或者命名空间设置,如Array33f::Zero(), ArrayXf::Zero(3), ArrayXXf::Zero(3, 4);Constant(rows, cols, value),用以初始元素相等的矩阵;Random(),用以初始一个随机矩阵;Identity(),用以初始一个单位阵,但这个函数只针对Matrix类型,不针对Array类型。LinSpaced(size, low, high)是一个特殊的函数,功能类似Matlab中的冒号用法,如下面实例:

    Eigen也提供了一个函数来进行初始设置,如setZero(), MatrixBase::setIdentity(), DenseBase::setLinSpaced()等,如我们要生成一个矩阵[0II0][0II0],可以有下面几种方法:

    上面介绍的例子中,一些方法如Zero(), Constant()等用以初始一个矩阵,我们可以认为这都是生成一个临时矩阵用以初始化,但事实是他们返回了一个所谓expression objects,而这不会引起多余的开销。例如我们想令一个2×32×3的矩阵行互换,即左乘一个

[0110][0110]

    上述的finished()是必要的,因为要获取一个实际的matrix来进行乘法计算,而不是一个用于初始某个矩阵的临时矩阵,它表明要在这个矩阵形成之后在用以做乘法计算。

Reductions,Visitors and Broadcasting

    Reductions可以看做一个作用在matrix或者vector上的函数,返回结果是一个值,比如之前介绍的sum(), prod(), trace(), minCoeff()等,这里再介绍一类Reductions—Norm computations。Eigen提供了范数类函数,norm(),如l2l2范数,squaredNorm(), lplp范数,lpNnorm<p>(),这里的p可以是1或者无穷。

    Reductions还有一类波尔类型,即一个矩阵可以对其中每一个元素进行大小比较。all(),如果每个元素都为真,则为真;any(),如果有一个元素为真,返回真;count(),为真元素的个数,如下面的例子所述:

    在使用maxCoeff()和minCoeff()函数是,可以寻求最大元素和最小元素,但如果我们想返回其位置,则需要给定相关参数,而参数的类型我们要使用Index类型,例如:

    Reductions还可以部分使用,即返回的不是一个值,而是一组值,可以看做是一种降维操作,所使用的函数为colwise(), rowwise()。例如

    Mat.colwise()理解为分别去看矩阵的每一列,然后再作用maxCoeff()函数,即求每一列的最大值。

需要注意的是,colwise返回的是一个行向量(列方向降维),rowwise返回的是一个列向量(行方向降维)。

下面说一下Broadcasting,他与Reductions的部分使用类似,但是他作用的类型是vector,相当于把vector按照某种方式进行扩充,再与矩阵进行计算。如

    Mat.colwise() += v;是核心语句,相当于[13216792]+[01010101][12693172]+[00001111]

    经过上面介绍,我们来一个实例:

    m.colwise() – v得到的是 [102184471][−12147084−1],再将得到的矩阵取列,再求l2l2范数的平方,即(m.colwise()-v).colwise().squaredNorm()得到的是[15053250][15053250],最后再求这一列的最小值,并记录位置。输出时m.col(index)便将最小的那一列输出。

Interfacing with Raw Buffers: the Map class

    Map类型,给用户一种选择:在已经存在的矩阵或者向量中,不必拷贝对象,而是直接在该对象的内存上进行操作。创建语句如下:

Map<Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>>;

    由此可见Map需要一个模板参数,还要确保给予一个指向该类型的指针,如

Map<MatrixXf> mf(pf, rows, cols);

    其中pf是有一个float指针。

    Map还可以接受两个参量,即Map<typename MatrixType, int MapOptions, tyepname StrideType>,MapOptions可以有Aligned和Unaligned两个选择;StrideType是对如何具体选择数据的一个选择,看下面两个实例

    上面实例解释了在创建Map的时候,如何决定内存的存储顺序;下面的实例将重点解释StrideType

 

    class Eigen::Stride< _OuterStrideAtCompileTime, _InnerStrideAtCompileTime >,这是stride的声明语法,有两个int类型的参数:OuterStride是按行存储矩阵每行之间的增量或者按列存储矩阵每列之间的增量,InnerStreide是按行存储矩阵每列之间的增量或者按列存储矩阵每行之间的增量。具体增量可以编译时给定,也可以运行时给定,即设置为Dynamic。上面的例子中,MatrixXi是按列存储,运行时stride的第一个参数是8,即每列之间相差8个位置;第二个参数是2,即每行直接相差2个位置;所以该矩阵存储为

array[0]array[2]array[4]array[8]array[10]array[12]array[16]array[18]array[20][array[0]array[8]array[16]array[2]array[10]array[18]array[4]array[12]array[20]]

    下面给一个使用Map变量的实例

    需要注意的是,这里m2,m2map和m2mapconst指向的是同一块内存。

    同样,Map也可以改变指向的内存,利用New语句,方法如下:

    这里v改变的只是指向内存的位置,并没有为此新开辟内存,而且在声明v的时候,应该已经有存在的内存区域被指向,即不会使用如下语句:

Map<Matrix3f> A(NULL);

Aliasing

    混淆往往出现在赋值号左右都用同一变量的时候,如mat = mat.transpose();解决方法之一是使用eval()函数,看下面两个实例

    如果有xxxInplace()函数,用这个函数可达到最高效的效果,如transposeInPlace()

    通常情况下,如果元素的表达只与该元素自己有关,那么一般不会引起混淆,如语句:mat = (2 * mat – MatrixXf::Identity(2, 2)).array().square();

    矩阵乘法是个特殊的操作,编译器会默认为其做避免混淆的处理,尽管这会引起额外的开销。如果你确定不会引起混淆,可以使用noalias()函数来避免额外开销。

    总结一下:

  1. 对于单个元素操作,如单个元素成绩或者加减操作一般不会引起混淆。
  2. 两个矩阵相乘时会做避免混淆的处理,如果你能确定不产生混淆,可以使用noalias()函数来避免额外开销。
  3. 其他情况下,都认为没有混淆存在,为了避免混淆产生,使用eval()函数或者xxxInPlace()函数,后者效率更高。

Storage orders

    Matrix 创建的矩阵默认是按列存储,如果想修改可以在创建矩阵的时候加入参数,如Matrix<int, 3, 4, ColMajor> Acolmajor; 或者 Matrix<int, 3, 4, RowMajor> Arowmajor;

    在选择存储顺序的时候,要有如下考虑:

  1. 你的使用者希望你用哪种方式存储,这可能跟其他处理库所要求的存储顺序有关。
  2. 与你所处理的实际算法有关。
  3. Eigen库在处理按列存储的矩阵时会更加高效。
0
0

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:5293次
    • 积分:92
    • 等级:
    • 排名:千里之外
    • 原创:3篇
    • 转载:6篇
    • 译文:0篇
    • 评论:0条
    文章分类