eigen：矩阵处理工具

2021年3月1日发(作者：季札挂剑翻译)

C++

开源矩阵计算工具

Eigen

< p>
简单用法（一）

矩阵转置函数

ose();

矩阵的特征值

alues();

矩阵 求逆矩阵

e();

特征向量

ctors ();

1

、

矩阵的定义

Eigen

中关于矩阵类的模板函数中，共有

6

个模板参数，但 是目前常用的只有前三个，如下所

示：

[cpp]

view plaincopy

1.

template

<

typename

_Scalar,

int

_Rows,

int

_Cols,

int

_Options,

in

t

_MaxRows,

int

_MaxCols>

2.

struct

traits

_Rows,

_Cols,

_Options,

_MaxRows,

_Max

Cols>

>

3.

.......

其前三个参数分别表示矩阵元素的类型，行数和列数。

矩阵定义时可以使用

Dynamic

来表示矩 阵的行列数为未知，例如：

typedef

Matrix

Dynamic, Dynamic

>

MatrixXd

;

在

Eigen

中也提供了很多常见的简化定义形式，例 如：

typedef Matrix< double , 3 , 1> Vector3d

注意：

（

1

）

Eigen

中无论是矩阵还是数组、

向量，

无论是静态矩阵还是动态矩阵都提供默认构造函数，

也就是你定义这些数据结构时都可以不用提供任何参数，其大小均由运行时来确定。

（

2

）矩阵的构造函 数中只提供行列数、元素类型的构造参数，而不提供元素值的构造，对于

比

较小

的、

固定长度

的

向量

提供初始化元素的定义，例如：

[cpp]

view plaincopy

1.

Vector2d

a(5.0,

6.0);

2.

Vector3d

b(5.0,

6.0,

7.0);

3.

Vector4d

c(5.0,

6.0,

7.0,

8.0);

2

、动态矩阵和静态矩阵

< p>
动态矩阵是指其大小在运行时确定，静态矩阵是指其大小在编译时确定，在

Eigen

中并未这样

称呼矩阵。具体可见如下两段代码：

代码段

1

：

[cpp]

view plaincopy

1.

#include

2.

#include

3.

using

namespace

Eigen;

4.

using

namespace

std;

5.

int

main()

6.

{

7.

MatrixXd

m

=

MatrixXd::Random(3,3);

8.

m

=

(m

+

MatrixXd::Constant(3,3,1.2))

*

50;

9.

cout

<<

=

<<

endl

<<

m

<<

endl;

10.

VectorXd

v(3);

11.

v

<<

1,

2,

3;

12.

cout

<<

*

v

=

<<

endl

<<

m

*

v

<<

endl;

13.

}

代码段

2

：

[cpp]

view plaincopy

1.

#include

2.

#include

3.

using

namespace

Eigen;

4.

using

namespace

std;

5.

int

main()

6.

{

7.

Matrix3d

m

=

Matrix3d::Random();

8.

m

=

(m

+

Matrix3d::Constant(1.2))

*

50;

9.

cout

<<

=

<<

endl

<<

m

<<

endl;

10.

Vector3d

v(1,2,3);

11.

cout

<<

*

v

=

<<

endl

<<

m

*

v

<<

endl;

12.

}

说明：

1

）

代码段

1

中

M atrixXd

表示任意大小的元素类型为

double

的矩阵变量，

其大小只有在运行时

被赋值之后 才能知道；

MatrixXd::Random

(3,3)

表示产生一个元素类型为

double< /p>

的

3*3

的临

时 矩阵对象。

2

）

代码段

2

中

Matrix3d

表示元素类型为

double

大小为

< br>3*3

的矩阵变量，其大小在编译时

就知道；

< p>

3

）

上例中向量的定义 也是类似，

不过这里的向量时列优先，

在

Eigen

中行优先的矩阵会在其名

字中包含有

row

，否则就是列优先

。

4

）向量只是一个特殊的矩阵，其一个维度为

< br>1

而已，如：

typedef Matrix< double , 3 , 1>

Vector3d

3

、矩阵元素的访问

在矩阵的访问中，

行索引总是作为第一个参数，

需注意

Eigen

中遵循大家的习惯让矩阵、

数组、

向量的下标都是从

0

开始。矩阵 元素的访问可以通过（）操作符完成，例如

m(2,3)

即是获 取矩

阵

m

的第

2

行第

3

列元素（注意行列数从

0

开始）。可参看如下代码：

[cpp]

view plaincopy

1.

#include

2.

#include

3.

using

namespace

Eigen;

4.

int

main()

5.

{

6.

MatrixXd

m(2,2);

7.

m(0,0)

=

3;

8.

m(1,0)

=

2.5;

9.

m(0,1)

=

-1;

10.

m(1,1)

=

m(1,0)

+

m(0,1);

11.

std::cout

<<

is

the

matrix

m:n

<<

m

<<

std::endl;

12.

VectorXd

v(2);

13.

v(0)

=

4;

14.

v(1)

=

v(0)

-

1;

15.

std::cout

<<

is

the

vector

v:n

<<

v

<<

std::endl;

16.

}

其输出结果为：

Here is the matrix m:

3 -1

2.5 1.5

Here is the vector v:

4

3

针对向量还提供

[]

操作符，注意矩阵则不可如此使用，原因为：在

C++

中

m[i,

j]

中逗号表达式

“i, j”的值始终都是“j”的值，即

m[i, j]

对于

C++

来讲就是

m[j];

4

、设置矩阵的元素

在

Eigen

中重载了

操作符，通过该操作符即可以一个一个元素的进行赋值，也可以一块一

块的赋值。另外也可以使用下标进行复制，例如下面两段代码：

代码段

1

[cpp]

view plaincopy

1.

Matrix3f

m;

2.

m

<<

1,

2,

3,

3.

4,

5,

6,

4.

7,

8,

9;

5.

std::cout

<<

m;

输出结果为：

1 2 3

4 5 6

7 8 9

代码段二（使用下标进行复制）

[cpp]

view plaincopy

1.

VectorXf

m_Vector_A;

2.

MatrixXf

m_matrix_B;

3.

int

m_iN

=-1;

4.

5.

bool

InitData(

int

pSrc[100][100],

int

iWidth,

int

iHeight)

6.

{

7.

if

(NULL

==

pSrc

||

iWidth

<=0

||

iHeight

<=

0)

8.

return

false

;

9.

m_iN

=

iWidth*iHeight;

10.

VectorXf

tmp_A(m_iN);

11.

MatrixXf

tmp_B(m_iN,

5);

12.

int

i

=0,

j=0,

iPos

=0;

13.

while

(i

14.

{

15.

j=0;

16.

while

(j

17.

{

18.

tmp_A(iPos)

=

pSrc[i][j]

*

log( (

float

)pSrc[i][j]

);

19.

20.

tmp_B(iPos,0)

=

pSrc[i][j]

;

21.

tmp_B(iPos,1)

=

pSrc[i][j]

*

i;

22.

tmp_B(iPos,2)

=

pSrc[i][j]

*

j;

23.

tmp_B(iPos,3)

=

pSrc[i][j]

*

i

*

i;

24.

tmp_B(iPos,4)

=

pSrc[i][j]

*

j

*

j;

25.

++iPos;

26.

++j;

27.

}

28.

++i;

29.

}

30.

m_Vector_A

=

tmp_A;

31.

m_matrix_B

=

tmp_B;

32.

}

5

、重置矩阵大小

< br>当前矩阵的行数、列数、大小可以通过

rows()

，< /p>

cols()

和

size()

< p>
来获取，对于动态矩阵可以通

过

resize()

函数来动态修改矩阵的大小

.

需注意：

(1)

< br>固定大小的矩阵是不能使用

resize

（）来修改矩阵 的大小；

(2)

resize

（）函数会析构掉原来的数据，因此调用

resize

（）函数之后将不能保证元素的

值不改变。

(3)

使用“=”操作符操作动态矩阵时，

< br>如果左右边的矩阵大小不等，

则左边的动态矩阵的大小

会 被修改为右边的大小。例如下面的代码段：

[cpp]

view plaincopy

1.

MatrixXf

a(2,2);

2.

std::cout

<<

is

of

size

<<

()

<<

<<

()

<<

std::endl;

3.

MatrixXf

b(3,3);

4.

a

=

b;

5.

std::cout

<<

is

now

of

size

<<

()

<<

<<

()

<<

std::endl;

输出结果为：

a is of size 2x2

a is now of size 3x3

6

、如何选择动态矩阵和静态矩阵？

Eigen

对于这问题的答案是：对于小矩阵（一般大小小于< /p>

16

）的使用固定大小的静态矩阵，它

可 以带来比较高的效率，对于大矩阵（一般大小大于

32

）建议使 用动态矩阵。

还需特别注意的是：如果特别大的矩阵使用了固 定大小的静态矩阵则可能造成栈溢出的问题

本文主要是

Eigen

中矩阵和向量的算术运算，

在

Eigen

中的这些算术运算重载了

C++

的

+

，

-

，

*

，

所以使用起 来非常方便。

1

、矩阵的运算

Eigen

提供

+

、

-

、一元操作符“

-

”、< /p>

+=

、

-=

，例 如：

二元操作符

+/-

表示两矩阵相加

（矩阵中对应元素相加

/

< p>
减，

返回一个临时矩阵）

：

B+C

或

B-C

；

一 元操作符

-

表示对矩阵取负（矩阵中对应元素取负，返回一个临 时矩阵）：

-C;

组合操作法

+=

或者

-=< /p>

表示（对应每隔元素都做相应操作）：

A += B

或者

A-=B

代码段

1

为矩阵的加减操作，代码如下：

[cpp]

view plaincopy

1.

#include

2.

#include

3.

using

namespace

Eigen;

4.

int

main()

5.

{

6.

Matrix2d

a;

7.

a

<<

1,

2,

8.

3,

4;

9.

MatrixXd

b(2,2);

10.

b

<<

2,

3,

11.

1,

4;

12.

std::cout

<<

+

b

=n

<<

a

+

b

<<

std::endl;

13.

std::cout

<<

-

b

=n

<<

a

-

b

<<

std::endl;

14.

std::cout

<<

a

+=

b;

<<

std::endl;

15.

a

+=

b;

16.

std::cout

<<

a

=n

<<

a

<<

std::endl;

17.

Vector3d

v(1,2,3);

18.

Vector3d

w(1,0,0);

19.

std::cout

<<

+

w

-

v

=n

<<

-v

+

w

-

v

<<

std::endl;

20.

}

输出结果为：

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

另外，

矩阵还提供与标量

（单一个数字）

的乘除操作，

表示每个元素都与该标量进行乘除操作。

例如：

二元操作符

*

在：

A*a

中表示矩阵

A

中的每隔元素都与数字

a

相乘，结 果放在一个临时矩阵中，

矩阵的值不会改变。

对于

a*A

、

A/a

、

A*=a

、

A /=a

也是一样，例如下面的代码：

[cpp]

view plaincopy

1.

#include

2.

#include

3.

using

namespace

Eigen;

4.

int

main()

5.

{

6.

Matrix2d

a;

7.

a

<<

1,

2,

8.

3,

4;

9.

Vector3d

v(1,2,3);

10.

std::cout

<<

*

2.5

=n

<<

a

*

2.5

<<

std::endl;

11.

std::cout

<<

*

v

=n

<<

0.1

*

v

<<

std::endl;

12.

std::cout

<<

v

*=

2;

<<

std::endl;

13.

v

*=

2;

14.

std::cout

<<

v

=n

<<

v

<<

std::endl;

15.

}

输出结果为：

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

中，算术操作例如

“操 作符

+

”并不会自己执行计算操作，他们只是返回一个“算术< /p>

表达式对象”，而实际的计算则会延迟到后面的赋值时才进行。这些不影响你的使用，它只 是

为了方便

Eigen

的优化。

2

、求矩阵的转秩、共轭矩阵、伴随矩阵。

可以通过

成员函数

transpose()

,

conjugate()

,

和

adjoint()

来完成，注意这些函数返 回

操作后的结果，而不会对原矩阵的元素进行直接操作，如果要让原矩阵的进行转换，则 需要使

用响应的

InPlace

函数， 例如：

transposeInPlace()

、

adjointInPlace()

之类。

例如下面的代码所示：

[cpp]

view plaincopy

1.

MatrixXcf

a

=

MatrixXcf::Random(2,2);

2.

cout

<<

is

the

matrix

an

<<

a

<<

endl;

3.

cout

<<

is

the

matrix

a^Tn

<<

ose()

<<

endl;

4.

cout

<<

is

the

conjugate

of

an

<<

ate()

<<

endl;

5.

cout

<<

is

the

matrix

a^*n

<<

t()

<<

endl;

输出结果为：

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)

3

、矩阵相乘、矩阵向量相乘

矩阵的相乘，矩阵与向量的相乘也是使用操作符

*

，共有

*

和

*=

< br>两种操作符，其用法可以参考如

下代码：

[cpp]

view plaincopy

1.

#include

2.

#include

3.

using

namespace

Eigen;

4.

int

main()

5.

{

6.

Matrix2d

mat;

7.

mat

<<

1,

2,

8.

3,

4;

9.

Vector2d

u(-1,1),

v(2,0);

10.

std::cout

<<

is

mat*mat:n

<<

mat*mat

<<

std::endl;

11.

std::cout

<<

is

mat*u:n

<<

mat*u

<<

std::endl;

12.

std::cout

<<

is

u^T*mat:n

<<

ose()*mat

<<

std::endl;

13.

std::cout

<<

is

u^T*v:n

<<

ose()*v

<<

std::endl;

14.

std::cout

<<

is

u*v^T:n

<<

u*ose()

<<

std::endl;

15.

std::cout

<<

multiply

mat

by

itself

<<

std::endl;

16.

mat

=

mat*mat;

17.

std::cout

<<

mat

is

mat:n

<<

mat

<<

std::endl;

18.

}

输出结果为：

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

本节主要涉 及

Eigen

的块操作以及

QR

分解

1

、矩阵的块操作

< br>1

）矩阵的块操作有两种使用方法，其定义形式为：

[cpp]

view plaincopy

1.

(i,j,p,q);

（

1

）

2.

3.

(i,j);

（

2

）

定义（

1

）表示返回从矩阵的

(i,

j)

开始，每行取

p

个元素，每列取

< p>
q

个元素所组成的临时新矩

阵对象，原矩阵的元素 不变。

定义（

2

）中

block

（

p, q

）可理解为一个

p

行

q

列的子矩阵，该定义表示从原矩阵中第

(i, j)

开始，获取一个

p

行

q

列的子矩阵，返回该子矩阵组成的临时

矩阵对象，原矩阵的元素不变。

详细使用情况，可参考下面的代码段：

[cpp]

view plaincopy

1.

#include

2.

#include

3.

using

namespace

std;

4.

int

main()

5.

{

6.

Eigen::MatrixXf

m(4,4);

7.

m

<<

1,

2,

3,

4,

8.

5,

6,

7,

8,

9.

9,10,11,12,

10.

13,14,15,16;

11.

cout

<<

in

the

middle

<<

endl;

12.

cout

<<

<2,2>(1,1)

<<

endl

<<

endl;

13.

for

(

int

i

=

1;

i

<=

3;

++i)

14.

{

15.

cout

<<

of

size

<<

i

<<

<<

i

<<

endl;

16.

cout

<<

(0,0,i,i)

<<

endl

<<

endl;

17.

}

18.

}

输出的结果为：

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

）

矩阵也提供了获取其指定行

/

列的函数，

其实获取某行

/

列也是一种特殊的获取子块。

可以通

过

.col()

和

.row()

来完成 获取指定列

/

行的操作，参数为列

/< /p>

行的索引。

注意：

（

1

）需与获取矩阵的行数

/

列数的函数（

rows(), cols()

）的进行区别，不要弄混淆。

（

2

）函数参数为响应行

/

< br>列的索引，需注意矩阵的行列均以

0

开始。

下面的代码段用于演示获取矩阵的指定行列：

[cpp]

view plaincopy

1.

#include

2.

#include

3.

using

namespace

std;

4.

int

main()

5.

{

6.

Eigen::MatrixXf

m(3,3);

7.

m

<<

1,2,3,

8.

4,5,6,

9.

7,8,9;

10.

cout

<<

is

the

matrix

m:

<<

endl

<<

m

<<

endl;

11.

cout

<<

Row:

<<

(1)

<<

endl;

12.

(2)

+=

3

*

(0);

13.

cout

<<

adding

3

times

the

first

column

into

the

third

co

lumn,

the

matrix

m

is:n

;

14.

cout

<<

m

<<

endl;

15.

}

输出结果为：

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

3

）

向量的块操作，

其实向量只是一个特殊的矩阵，

但是

Eige n

也为它单独提供了一些简化的块

操作，如下三种形式：

获取向量的前

n

个元素：

(n);

获取向量尾部的

n

个元素：

(n);

获取从向量的第

i

个元素开始的

n

个元素：

t( i,n);

其用法可参考如下代码段：

[cpp]

view plaincopy

1.

#include

2.

#include

3.

using

namespace

std;

4.

int

main()

5.

{

6.

Eigen::ArrayXf

v(6);

7.

v

<<

1,

2,

3,

4,

5,

6;

8.

cout

<<

=

<<

endl

<<

(3)

<<

endl

<<

endl;

9.

cout

<<

=

<<

endl

<<

<3>()

<<

endl

<<

endl;

10.

t(1,4)

*=

2;

11.

cout

<<

't(1,4)

*=

2',

v

=

<<

endl

<<

v

<<

endl

;

12.

}

输出结果为：

(3) =

1

2

3

<3>() =

4

5

6

after 't(1,4) *= 2', v =

1

4

6

8

10

6

2

、

QR

分解

Eigen

的

< br>QR

分解非常绕人，它总共提供了下面这些矩阵的分解方式：

Decomposition

PartialPivL U

FullPivLU

Method

partialPivLu()

fullPivLu()

Requirements on

the matrix

Invertible

None

None

None

Speed

++

-

++

+

-

Accuracy

+

+++

+

++

+++

+

++

HouseholderQR

householderQr()

ColPivHouseholderQR

colPivHou seholderQr()

FullPivHouseholderQR

LLT

fullPivHouseholderQr()

None

llt()

Positive definite

+++

Positive or

negative

semidefinite

+++

LDLT

ldlt()

由于我只用到 了

QR

分解，而且

Eigen

的

QR

分解开始使用时确实不容易入手，因此这 里只提供

了

householderQR

的分解方式的演示代码：

[cpp]

view plaincopy

1.

void

QR2()

2.

{

3.

Matrix3d

A;

4.

A<<1,1,1,

5.

2,-1,-1,

6.

2,-4,5;

7.

8.

HouseholderQR

qr;

9.

e(A);

10.

MatrixXd

R

=

QR().triangularView();

11.

MatrixXd

Q

=

olderQ();

12.

std::cout

<<

HouseholderQR- ------------------------------

--------- -----

<<

std::endl;

13.

std::cout

<<

<<

std::endl

<

<<

std::endl

<<

std::endl

;

14.

std::cout

<<

<<

std::endl

<<

QR()

<<

st

d::endl

<<

std::endl;

15.

std::cout

<<

<<

std::endl

<

<<

std::endl

<<

std::endl;

16.

std::cout

<<

<<

std::endl

<

<<

std::endl

<<

std::endl

;

17.

std::cout

<<

<<

std::endl

<

<<

std::endl

<<

std::en

dl;

18.

}

输出结果为：

3

、一个矩阵使用的例子：用矩阵操 作完成二维高斯拟合，并求取光斑中心

下面的代码段是一个使 用

Eigen

的矩阵操作完成二维高斯拟合求取光点的代码例子 ，关于二维

高斯拟合求取光点的详细内容可参考：

/hjx_1 000/article/details/8490653

[cpp]

view plaincopy

1.

bool

GetCentrePoi nt(

float

&

x0,

float

&

y0)

2.

{

3.

if

(m_iN<=0)

4.

return

false

;

5.

//QR

分解

6.

HouseholderQR

qr;

7.

e(m_matrix_B);

8.

MatrixXf

R

=

QR().triangularView();

9.

MatrixXf

Q

=

olderQ();

10.

11.

//

块操作，获取向量或矩阵的局部

12.

VectorXf

S;

13.

S

=

(ose()*

m_Vector_A).head(5);

14.

MatrixXf

R1;

15.

R1

=

(0,0,5,5);

16.

17.

VectorXf

C;

18.

C

=

e()

*

S;

19.

20.

x0

=

-0.5

*

C[1]

/

C[3];

21.

y0

=

-0.5

*

C[2]

/

C[4];