c++# #eigen# #线性代数# #矩阵# #lib#

• http://rswiki.org/doku.php?id=%E7%A8%8B%E5%BA%8F%E8%AE%BE%E8%AE%A1:%E7%9F%A9%E9%98%B5%E5%BA%93:eigen
• Eigen是一个轻量级的矩阵库,除了稀疏矩阵不成熟以外,其他的矩阵和向量操作都比较完善,而且速度不错.
• 2021/11/19 11:10:19

Eigen是一个轻量级的矩阵库,除了稀疏矩阵不成熟以外,其他的矩阵和向量操作都比较完善,而且速度不错.

#include <Eigen/Dense>
using namespace Eigen;
int main()
{
	MatrixXf m1(3,4);   //动态矩阵，建立3行4列。
	MatrixXf m2(4,3);	//4行3列，依此类推。
	MatrixXf m3(3,3);
 
	Vector3f v1;		//若是静态数组，则不用指定行或者列
	/* 初始化 */
	Matrix3d m = Matrix3d::Random();
	m1 = MatrixXf::Zero(3,4);		//用0矩阵初始化,要指定行列数
	m2 = MatrixXf::Zero(4,3);
	m3 = MatrixXf::Identity(3,3);	//用单位矩阵初始化
	v1 = Vector3f::Zero();			//同理，若是静态的，不用指定行列数
 
	m1 << 1,0,0,1,		//也可以以这种方式初始化
		1,5,0,1,
		0,0,9,1;
	m2 << 1,0,0,
		0,4,0,
		0,0,7,
		1,1,1;
}

使用Map函数，可以实现Eigen的矩阵和c++中的数组直接转换,语法如下：

//@param MatrixType 矩阵类型
//@param MapOptions 指的是指针是否对齐，Aligned, or Unaligned. The default is Unaligned.
//@param StrideType 
/*
	Map<typename MatrixType,
		int MapOptions,
		typename StrideType>
*/
//pf是float数组 mf是flost矩阵
Map<MatrixXf> mf(pf,rows,columns);
//mi是const Vector4i，pi是int*
Map<const Vector4i> mi(pi);

eigen重载了基础的+ - * / += -= *= /=
*可以表示标量和矩阵或者矩阵和矩阵

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix2d a;
 a << 1, 2,
 3, 4;
 MatrixXd b(2,2);
 b << 2, 3,
 1, 4;
 std::cout << "a + b =\n" << a + b << std::endl;
 std::cout << "a - b =\n" << a - b << std::endl;
 std::cout << "Doing a += b;" << std::endl;
 a += b;
 std::cout << "Now a =\n" << a << std::endl;
 Vector3d v(1,2,3);
 Vector3d w(1,0,0);
 std::cout << "-v + w - v =\n" << -v + w - v << std::endl;
 a << 1, 2,
 3, 4;
 Vector3d v(1,2,3);
 std::cout << "a * 2.5 =\n" << a * 2.5 << std::endl;
 std::cout << "0.1 * v =\n" << 0.1 * v << std::endl;
 std::cout << "Doing v *= 2;" << std::endl;
 v *= 2;
 std::cout << "Now v =\n" << v << std::endl;
}
/*
	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 * 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
*/

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
	Vector3d v(1,2,3);
	Vector3d w(0,1,2);
	 cout << "Dot product: " << v.dot(w) << endl;
	double dp = v.adjoint()*w; // automatic conversion of the inner product to a scalar
	 cout << "Dot product via a matrix product: " << dp << endl;
	 cout << "Cross product:\n" << v.cross(w) << endl;
}
/*
Dot product: 8
Dot product via a matrix product: 8
Cross product:
1
-2
1
*/

解 ax = b

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
Matrix3f A;
 Vector3f b;
 A << 1,2,3, 4,5,6, 7,8,10;
 b << 3, 3, 4;
 cout << "Here is the matrix A:\n" << A << endl;
 cout << "Here is the vector b:\n" << b << endl;
 Vector3f x = A.colPivHouseholderQr().solve(b);
 cout << "The solution is:\n" << x << endl;
}
/* output
	Here is the matrix A:
	 1  2  3
	 4  5  6
	 7  8 10
	Here is the vector b:
	3
	3
	4
	The solution is:
	-2
	1
	1
*/

除了colPivHouseholderQr，还有以下的函数：

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
Matrix2f A;
 A << 1, 2, 2, 3;
 cout << "Here is the matrix A:\n" << A << endl;
SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
if (eigensolver.info() != Success) abort();
 cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl;
 cout << "Here's a matrix whose columns are eigenvectors of A \n"
 << "corresponding to these eigenvalues:\n"
 << eigensolver.eigenvectors() << endl;
}
/*
	Here is the matrix A:
	1 2
	2 3
	The eigenvalues of A are:
	-0.236
	4.24
	Here's a matrix whose columns are eigenvectors of A 
	corresponding to these eigenvalues:
	-0.851 -0.526
	 0.526 -0.851
*/

小矩阵（4 * 4及以下）eigen会自动优化，默认采用LU分解，效率不高

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
Matrix3f A;
 A << 1, 2, 1,
 2, 1, 0,
 -1, 1, 2;
 cout << "Here is the matrix A:\n" << A << endl;
 cout << "The determinant of A is " << A.determinant() << endl;
 cout << "The inverse of A is:\n" << A.inverse() << endl;
}
/*
	Here is the matrix A:
	1 2 1
	2 1 0
	-1 1 2
	The determinant of A is -3
	The inverse of A is:
	-0.667      1  0.333
	  1.33     -1 -0.667
		-1      1      1
*/

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
MatrixXf A = MatrixXf::Random(3, 2);
 cout << "Here is the matrix A:\n" << A << endl;
VectorXf b = VectorXf::Random(3);
 cout << "Here is the right hand side b:\n" << b << endl;
 cout << "The least-squares solution is:\n"
 << A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b) << endl;
}
/*
	Here is the matrix A:
	  0.68  0.597
	-0.211  0.823
	 0.566 -0.605
	Here is the right hand side b:
	-0.33
	0.536
	-0.444
	The least-squares solution is:
	-0.67
	0.314
*/

稀疏矩阵的头文件包括：

  #include <Eigen/SparseCore>
  #include <Eigen/SparseCholesky>
  #include <Eigen/IterativeLinearSolvers>
  #include <Eigen/Sparse>

初始化有两种方式：
1.使用三元组插入

typedef Eigen::Triplet<double> T;
std::vector<T> tripletList;
triplets.reserve(estimation_of_entries); //estimation_of_entries是预估的条目
for(...)
{
	tripletList.push_back(T(i,j,v_ij));//第 i,j个有值的位置的值
}
SparseMatrixType mat(rows,cols);
mat.setFromTriplets(tripletList.begin(), tripletList.end());
// mat is ready to go!

2.直接将已知的非0值插入

SparseMatrix<double> mat(rows,cols);
mat.reserve(VectorXi::Constant(cols,6));
for(...)
{
	// i,j 个非零值 v_ij != 0
	mat.insert(i,j) = v_ij;
}
mat.makeCompressed(); // optional

稀疏矩阵支持大部分一元和二元运算:

sm1.real() sm1.imag() -sm1 0.5*sm1
sm1+sm2 sm1-sm2 sm1.cwiseProduct(sm2)

二元运算中，稀疏矩阵和普通矩阵可以混合使用

//dm表示普通矩阵
dm2 = sm1 + dm1;

也支持计算转置矩阵和伴随矩阵