eigen 笔记2

4. Block operations

1) Using block operations

Block of size(p, q), starting at (i, j)

dynamic-size block expression: matrix.block(i, j, p, q);

fixed-size block expression: matrix.block<p, q>(i, j);

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;
  }
}

Output:

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

.block() 的返回值可以是左值, 也可以是右值.

2) Columns and rows

matrix.row(i), matrix.col(j)

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:
";
  cout << m << endl;
}

3) Corner-related operations

太复杂, 没兴趣.

4) Block operations for vectors

同上, 用时再查.

5. Advanced initialization

1) The comma initializer

RowVectorXd vec1(3);
vec1 << 1, 2, 3;

RowVectorXd vec2(4);
vec2 << 1, 4, 9, 16;

RowVectorXd joined(7);
joined << vec1, vec2;

MatrixXf matA(2, 2);
matA << 1, 2, 3, 4;

MatrixXf matB(4, 4);
matB << matA, matA/10, matA/10, matA;

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;   

2) Special matrices and arrays

Zero(), Constant(value), Constant(rows, cols, value), Random(), Identity(), LinSpaced(size, low, high), 

Identity() only for matrix;

LinSpaced only for vectors and one-dimensional arrays; high也包含在数组中. 默认都是列向量.

std::cout << "A fixed-size array:
";
Array33f a1 = Array33f::Zero();

std::cout << "A one-dimensional dynamic-size array:
";
ArrayXf a2 = ArrayXf::Zero(3);

std::cout << "A two-dimensional dynamic-size array:
";ArrayXXf a3 = ArrayXXf::Zero(3, 4);

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();

Output:
  Degrees   Radians      Sine    Cosine
        0         0         0         1
       10     0.175     0.174     0.985
       20     0.349     0.342      0.94
       30     0.524       0.5     0.866
       40     0.698     0.643     0.766
       50     0.873     0.766     0.643
       60      1.05     0.866       0.5
       70      1.22      0.94     0.342
       80       1.4     0.985     0.174
       90      1.57         1 -4.37e-08

类似的, 有 setZero(), MatrixBase::setIdentity(), DenseBase::setLinSpaced()

3) Usage as temporary objects

int main()
{
  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;
}

comma-initializer

MatrixXf mat = MatrixXf::Random(2, 3);
std::cout << mat << std::endl << std::endl;
mat = (MatrixXf(2,2) << 0, 1, 1, 0).finished() * mat;

注意到, .finished() 是必须的.

6. Reductions, visitors and broadcasting

1) Reducations

sum(), prod(), mean(), minCoeff(), maxCoeff(), trace()

2) Norm computations

应该用不到.

3) Boolean reductions

all(), any(), count()

int main()
{
  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;
}

Output:
(a > 0).all()   = 1
(a > 0).any()   = 1
(a > 0).count() = 4

(a > 2).all()   = 0
(a > 2).any()   = 1
(a > 2).count() = 2

4) Visitors

maxCoeff(&x, &y), minCoeff(&x, &y)

x, y 即为行索引和列索引, MatrixXf::Index, 注意, 不是 int

  //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);

5) Partial reductions

colwise(), rowwise()

int main() {
    Eigen::MatrixXf mat(2, 4);
    mat << 1, 2, 6, 9,
           3, 1, 7, 2;
    cout << "Column's maximum: " << endl
        << mat.colwise().maxCoeff() << endl;
    cout << "Row's maximum: " << endl
        << mat.rowwise().maxCoeff() << endl;
}

Output:
Column's maximum:
3 2 7 9

Row's maximum:
9
7

combining partial reductions with other operations

int main() {
    MatrixXf mat(2, 4);
    mat << 1, 2, 6, 9,
           3, 1, 7, 2;
    MatrixXf::Index maxIndex;
    float maxNorm = mat.colwise().sum().maxCoeff(&maxIndex);

    cout << "Maximum sum at position: " << maxIndex << endl;
    cout << "The corresponding vector is: " << endl;
    cout << mat.col(maxIndex) << endl;
    cout << "And its sum is: " << maxNorm << endl;
}

Output:
Maximum sum at position 2
The corresponding vector is: 
6
7
And its sum is is: 13

6) Broadcasting

int main() {
    MatrixXf mat(2, 4);
    VectorXf v(2);

    mat << 1, 2, 6, 9,
           3, 1, 7, 2;
    v << 0, 1;
    w << 0, 1, 2, 3;

    mat.colwise() += v;
   
    cout << "Broadcasting result: " << endl;
    cout << mat << endl;

    mat.rowwise() += w.transpose();
    cout << "Broadcasting result: " << endl;
    cout << mat << endl;
}

Output:

Broadcasting result: 
1 2 6 9
4 2 8 3

Broadcasting result: 
1 3 8 12
4 3 10 6

7) Combining broadcasting with other operations

略.

7. Interfacing with raw buffers: the Map class

1) Map types and declaring Map variables

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

RowsAtCompileTime 和 ColsAtCompileTime 是生成矩阵的行数和列数

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

pf 指向一个定义 coefficient array 的 region 的起始, 这里注意 pf 是一个指针, 指向的是一个 coefficient, rows 和 columns 是该 mf 的行和列, 默认全部.

这里有一个问题, Map 模板参数指定的 RowsAtCompileTime 和 ColsAtCompileTime 与括号中的 rows 和 columns 有什么关系?

a. Matrix

Dynamic

Map<Matrix<int, Dynamic, Dynamic> > 即 Map<MatrixXi>, 此时生成 matrix 的行和列由 rows 和 columns 指定

这里还要注意一下, pf 的属性, 如果 pf 是指向一个内存区域, 内存区域的数字怎么获取, 要看这个区域存储的是什么.

存储 array, 假设 array 长度为 15, 要转换成3行4列矩阵, 那么取前12个数

存储 matrix, 注意 matrix 默认是按列存储的, 所以会取按列排布的前12个数, 也就是前4列

fixed-size

Map<Matrix<int, 3, 4> > , 此时生成 Matrix 的行和列已定, rows 和 columns 要与 3, 4 保持一致, 或者省略否则报错.

b. Vector

Dynamic

Map<Matrix<int, 1, Dynamic> > or Map<Matrix<int, Dynamic, 1> > 

Vector 有一个特别地方, rows 和 columns 可以只写一个, 比如 Map<VectorXi> m2map(p, m2.size()), 也可以全写, 但是维度要与 Map 中的一致.

fixed-size

注意维度一致.

Map<const Vector4i> mi(pi);

声明一个 fixed-size read-only vector.

Map 还有几个可选的模板参数

Map<typename MatrixType, int MapOptions, typename StrideType>

MapOptions: Aligned or Unaligned, 默认 Unaligned

StrideType: 决定 array 的排布, 使用 Stride class

int array[8];
for(int i = 0; i < 8; ++i) array[i] = i;
cout << "Column-major:
" << Map<Matrix<int,2,4> >(array) << endl;
cout << "Row-major:
" << Map<Matrix<int,2,4,RowMajor> >(array) << endl;
cout << "Row-major using stride:
" <<
  Map<Matrix<int,2,4>, Unaligned, Stride<1,4> >(array) << endl;

Output:
    
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) Using Map variables

注意: a Map type is not identical to its Dense equivalent.

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:

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

这里注意一下, m2map, m2mapconst 与 m2 共享内存空间, 所以对 m2map 的修改, 会同步到 m2 和 m2mapconst.

但是, 如果如果 *p 指向的是 array 中的地址, 那么 map 转化出的 matrix 不会与 array 共享内存空间.

3) Changing the mapped array

用 new 可以对声明后的 Map object 进行修改

int data[] = {1,2,3,4,5,6,7,8,9};
Map<RowVectorXi> v(data,4);
cout << "The mapped vector v is: " << v << "
";
new (&v) Map<RowVectorXi>(data+4,5);
cout << "Now v is: " << v << "
";

Output:
The mapped vector v is: 1 2 3 4
Now v is: 5 6 7 8 9

所以也可以声明未知 array 的 map object

Map<Matrix3f> A(NULL);  // don't try to use this matrix yet!
VectorXf b(n_matrices);
for (int i = 0; i < n_matrices; i++)
{
  new (&A) Map<Matrix3f>(get_matrix_pointer(i));
  b(i) = A.trace();
}

8. Reshape and Slicing

1) Reshape

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:
v1:
1 4 7 2 5 8 3 6 9
v2:
1 2 3 4 5 6 7 8 9

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:
M2:
 1  4
 7 10
 2  5
 8 11
 3  6
 9 12

注意这里, 存储方式的不同, 会导致输出不同.

2) Slicing

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:
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 << "
";
Map<MatrixXf,0,OuterStride<> > M2(M1.data(), M1.rows(), (M1.cols()+2)/3, OuterStride<>(M1.outerStride()*3));
cout << "1 column over 3:" << endl << M2 << "
";
typedef Matrix<float,Dynamic,Dynamic,RowMajor> RowMajorMatrixXf;
RowMajorMatrixXf M3(M1);
cout << "Row major input:" << endl << M3 << "
";
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 << "
";

Output:
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

注意 M1.data(), 返回一个指针.

9. Aliasing

m = m.transpose() 会造成 aliasing

a) m = m.transpose().eval();

b) m = m.transposeInPlace();

component-wise operations is safe.

mat = 2 * mat;

mat = mat - MatrixXf::Identity(2, 2);

matB = matA * matA; // Simple bu not quite as efficient

matB.noalias() = matA * matA; // Use this

10. Storage orders

默认 column-major, 也可以选择 row-major

11. Quick reference guide

Eigen/LU    Inverse, determinant, LU decompositions with solver (FullPivLU, PartialPivLU)

Eigen/SVD    SVD decompositions with least-squares solver (JacobiSVD, BDCSVD)

Eigen/QR    decomposition with solver (HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR)

...

Finished.