Eigen矩阵库——初步用法（三）

星期四, 05. 一月 2017 12:49下午

本节主要涉及Eigen的块操作以及QR分解。Eigen的QR分解非常绕人，搞了很久才搞明白是怎么回事，最后是一个使用Eigen的矩阵操作完成二维高斯拟合求取光点的代码例子，关于二维高斯拟合求取光点的详细内容可参考：链接

矩阵的块操作

矩阵的块表示

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

matrix.block(i,j,p,q);        	   (1)
matrix.block<p,q>(i,j)   		   (2)

• 定义（1）表示返回从矩阵的(i, j)开始，每行取p个元素，每列取q个元素所组成的临时新矩阵对象，原矩阵的元素不变。

• 定义（2）中block<p, q>可理解为一个p行q列的子矩阵，该定义表示从原矩阵中第(i, j)开始，获取一个p行q列的子矩阵，返回该子矩阵组成的临时 矩阵对象，原矩阵的元素不变。

演示代码如下：

#include<iostream>
#include<eigen3/Eigen/Dense>
using namespace std;
int main()
{
  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;
  }
}

运行结果：

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

通过上述方式获取的子矩阵即可以作为左值也可以作为右值，也就是即可以用这个子矩阵给其他矩阵赋值，也可以给这个子矩阵对象赋值。

获取其指定行/列

矩阵也提供了获取其指定行/列的函数，其实获取某行/列也是一种特殊的获取子块。可以通过 .col()和 .row()来完成获取指定列/行的操作，参数为列/行的索引。 注意：

• 需与获取矩阵的行数/列数的函数 rows(),cols()的进行区别，不要弄混淆。

• 函数参数为响应行/列的索引，需注意矩阵的行列均以0开始。

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

#include<eigen3/Eigen/Dense> 
#include <iostream>  
using namespace std; 
int main() 
{ 
	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; 
} 

输出结果为：

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

向量的块操作

向量的块操作，其实向量只是一个特殊的矩阵，但是Eigen也为它单独提供了一些简化的块操作，如下三种形式：

• 获取向量的前n个元素：vector.head(n);
• 获取向量尾部的n个元素：vector.tail(n);
• 获取从向量的第i个元素开始的n个元素：vector.segment(i,n);

其用法可参考如下代码段

#include<eigen3/Eigen/Dense> 
#include <iostream> 
using namespace std; 
int main() 
{ 
    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; 
} 

输出结果为：

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

QR分解

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

由于作者只用到了QR分解，而且Eigen的QR分解开始使用时确实不容易入手，因此这里只提供了householderQR的分解方式的演示代码：

void QR2() 
{ 
    Matrix3d A; 
    A<<1,1,1, 
        2,-1,-1,  
        2,-4,5;  
  
    HouseholderQR<Matrix3d> qr; 
    qr.compute(A); 
    MatrixXd R = qr.matrixQR().triangularView<Upper>(); 
    MatrixXd Q =  qr.householderQ(); 
    std::cout << "QR2(): HouseholderQR---------------------------------------------"<< std::endl; 
    std::cout << "A "<< std::endl <<A << std::endl << std::endl; 
    std::cout <<"qr.matrixQR()"<< std::endl << qr.matrixQR() << std::endl << std::endl; 
    std::cout << "R"<< std::endl <<R << std::endl << std::endl; 
    std::cout << "Q "<< std::endl <<Q << std::endl << std::endl; 
    std::cout <<"Q*R" << std::endl <<Q*R << std::endl << std::endl; 
} 

运行结果如下：

QR2(): HouseholderQR---------------------------------------------
A 
 1  1  1
 2 -1 -1
 2 -4  5

qr.matrixQR()
-3  3 -3
0.5 -3  3
0.5 -1 -3

R
-3  3 -3
 0 -3  3
 0  0 -3

Q 
-0.333333 -0.666667 -0.666667
-0.666667 -0.333333  0.666667
-0.666667  0.666667 -0.333333

Q*R
 1  1  1
 2 -1 -1
 2 -4  5

原文来源

C++开源矩阵计算工具——Eigen的简单用法（三） 逍遥子_