R语言学习RcppEigen进行矩阵运算

前面我们介绍了一些基本的rcpp的用法：让你的r代码更快——rcpp入门，但用基础的rcpp来进行矩阵运算还是非常麻烦，没有现成的函数来让我们使用。

这时我们就想到：是否可以调用别的库来解决矩阵运算的一些问题呢？这就需要我们的rcppeigen包，也就是c++中的eigen库。

这些矩阵的运算在进行模拟时会时常遇到，所以可以说是非常重要的一项技能，下面我们就给予一个现有的对矩阵处理的代码来说明其用法。

创建cpp文件

其创建方式可以参考上篇博客：让你的r代码更快——rcpp入门

代码示例

然后我们定义一个

来做矩阵乘法并求其迹（trace）的函数。

// [[rcpp::depends(rcppeigen)]]
#include <rcppeigen.h>
using namespace eigen;
using namespace std;

//[[rcpp::export]]
double myfun (matrixxd x, matrixxd y) {
  double z;
  
  z = (x.adjoint() * y).trace();
  cout << z << endl;
 
  return z;
}

前三行表示载入eigen库

// [[rcpp::depends(rcppeigen)]]
#include <rcppeigen.h>
using namespace eigen;

里面的转置函数adjoint()，求迹函数trace()，都需要用到这个库，如果不使用命名空间eigen后面库里面就要这样用eigen::adjoint()，eigen::trace()。

后面我们使用using namespace std;则是因为cout需要用到，这个可以在运行函数的时候展现我们的中间变量，也是一个比较有用的操作，当然如果不需要的话，就可以不用命名变量空间：std。

下面就是我们的函数：

//[[rcpp::export]]
double myfun (matrixxd x, matrixxd y) {
  double z;
  
  z = (x.adjoint() * y).trace();
  cout << z << endl;
 
  return z;
}

//[[rcpp::export]]为我们需要导出到r中的时候需要添加，double型的矩阵在eigen中命名为matrixxd，整型矩阵为matrixxi；类似，对应的向量命名方式为：vectorxd与vectorxi。

里面的内容就是我们按照公式敲的函数。

下面我们介绍一些eigen库中的其它一些矩阵操作。

其他矩阵操作

这部分原文：a simple quickref for eigen

命名

matrix<double, 3, 3> a;               // fixed rows and cols. same as matrix3d.
matrix<double, 3, dynamic> b;         // fixed rows, dynamic cols.
matrix<double, dynamic, dynamic> c;   // full dynamic. same as matrixxd.
matrix<double, 3, 3, rowmajor> e;     // row major; default is column-major.
matrix3f p, q, r;                     // 3x3 float matrix.
vector3f x, y, z;                     // 3x1 float matrix.
rowvector3f a, b, c;                  // 1x3 float matrix.
vectorxd v;                           // dynamic column vector of doubles
double s;                            

基础用法

// basic usage
// eigen          // matlab           // comments
x.size()          // length(x)        // vector size
c.rows()          // size(c,1)        // number of rows
c.cols()          // size(c,2)        // number of columns
x(i)              // x(i+1)           // matlab is 1-based
c(i,j)            // c(i+1,j+1)       //

a.resize(4, 4);   // runtime error if assertions are on.
b.resize(4, 9);   // runtime error if assertions are on.
a.resize(3, 3);   // ok; size didn't change.
b.resize(3, 9);   // ok; only dynamic cols changed.
                  
a << 1, 2, 3,     // initialize a. the elements can also be
     4, 5, 6,     // matrices, which are stacked along cols
     7, 8, 9;     // and then the rows are stacked.
b << a, a, a;     // b is three horizontally stacked a's.
a.fill(10);       // fill a with all 10's.

定义矩阵

// eigen                                    // matlab
matrixxd::identity(rows,cols)               // eye(rows,cols)
c.setidentity(rows,cols)                    // c = eye(rows,cols)
matrixxd::zero(rows,cols)                   // zeros(rows,cols)
c.setzero(rows,cols)                        // c = zeros(rows,cols)
matrixxd::ones(rows,cols)                   // ones(rows,cols)
c.setones(rows,cols)                        // c = ones(rows,cols)
matrixxd::random(rows,cols)                 // rand(rows,cols)*2-1            // matrixxd::random returns uniform random numbers in (-1, 1).
c.setrandom(rows,cols)                      // c = rand(rows,cols)*2-1
vectorxd::linspaced(size,low,high)          // linspace(low,high,size)'
v.setlinspaced(size,low,high)               // v = linspace(low,high,size)'
vectorxi::linspaced(((hi-low)/step)+1,      // low:step:hi
                    low,low+step*(size-1))  //

对矩阵的一些基础操作1

// matrix slicing and blocks. all expressions listed here are read/write.
// templated size versions are faster. note that matlab is 1-based (a size n
// vector is x(1)...x(n)).
// eigen                           // matlab
x.head(n)                          // x(1:n)
x.head<n>()                        // x(1:n)
x.tail(n)                          // x(end - n + 1: end)
x.tail<n>()                        // x(end - n + 1: end)
x.segment(i, n)                    // x(i+1 : i+n)
x.segment<n>(i)                    // x(i+1 : i+n)
p.block(i, j, rows, cols)          // p(i+1 : i+rows, j+1 : j+cols)
p.block<rows, cols>(i, j)          // p(i+1 : i+rows, j+1 : j+cols)
p.row(i)                           // p(i+1, :)
p.col(j)                           // p(:, j+1)
p.leftcols<cols>()                 // p(:, 1:cols)
p.leftcols(cols)                   // p(:, 1:cols)
p.middlecols<cols>(j)              // p(:, j+1:j+cols)
p.middlecols(j, cols)              // p(:, j+1:j+cols)
p.rightcols<cols>()                // p(:, end-cols+1:end)
p.rightcols(cols)                  // p(:, end-cols+1:end)
p.toprows<rows>()                  // p(1:rows, :)
p.toprows(rows)                    // p(1:rows, :)
p.middlerows<rows>(i)              // p(i+1:i+rows, :)
p.middlerows(i, rows)              // p(i+1:i+rows, :)
p.bottomrows<rows>()               // p(end-rows+1:end, :)
p.bottomrows(rows)                 // p(end-rows+1:end, :)
p.topleftcorner(rows, cols)        // p(1:rows, 1:cols)
p.toprightcorner(rows, cols)       // p(1:rows, end-cols+1:end)
p.bottomleftcorner(rows, cols)     // p(end-rows+1:end, 1:cols)
p.bottomrightcorner(rows, cols)    // p(end-rows+1:end, end-cols+1:end)
p.topleftcorner<rows,cols>()       // p(1:rows, 1:cols)
p.toprightcorner<rows,cols>()      // p(1:rows, end-cols+1:end)
p.bottomleftcorner<rows,cols>()    // p(end-rows+1:end, 1:cols)
p.bottomrightcorner<rows,cols>()   // p(end-rows+1:end, end-cols+1:end)

基础操作2

// of particular note is eigen's swap function which is highly optimized.
// eigen                           // matlab
r.row(i) = p.col(j);               // r(i, :) = p(:, j)
r.col(j1).swap(mat1.col(j2));      // r(:, [j1 j2]) = r(:, [j2, j1])

// views, transpose, etc;
// eigen                           // matlab
r.adjoint()                        // r'
r.transpose()                      // r.' or conj(r')       // read-write
r.diagonal()                       // diag(r)               // read-write
x.asdiagonal()                     // diag(x)
r.transpose().colwise().reverse()  // rot90(r)              // read-write
r.rowwise().reverse()              // fliplr(r)
r.colwise().reverse()              // flipud(r)
r.replicate(i,j)                   // repmat(p,i,j)

矩阵基础运算1

// all the same as matlab, but matlab doesn't have *= style operators.
// matrix-vector.  matrix-matrix.   matrix-scalar.
y  = m*x;          r  = p*q;        r  = p*s;
a  = b*m;          r  = p - q;      r  = s*p;
a *= m;            r  = p + q;      r  = p/s;
                   r *= q;          r  = s*p;
                   r += q;          r *= s;
                   r -= q;          r /= s;

矩阵基础运算2

// vectorized operations on each element independently
// eigen                       // matlab
r = p.cwiseproduct(q);         // r = p .* q
r = p.array() * s.array();     // r = p .* s
r = p.cwisequotient(q);        // r = p ./ q
r = p.array() / q.array();     // r = p ./ q
r = p.array() + s.array();     // r = p + s
r = p.array() - s.array();     // r = p - s
r.array() += s;                // r = r + s
r.array() -= s;                // r = r - s
r.array() < q.array();         // r < q
r.array() <= q.array();        // r <= q
r.cwiseinverse();              // 1 ./ r
r.array().inverse();           // 1 ./ r
r.array().sin()                // sin(r)
r.array().cos()                // cos(r)
r.array().pow(s)               // r .^ s
r.array().square()             // r .^ 2
r.array().cube()               // r .^ 3
r.cwisesqrt()                  // sqrt(r)
r.array().sqrt()               // sqrt(r)
r.array().exp()                // exp(r)
r.array().log()                // log(r)
r.cwisemax(p)                  // max(r, p)
r.array().max(p.array())       // max(r, p)
r.cwisemin(p)                  // min(r, p)
r.array().min(p.array())       // min(r, p)
r.cwiseabs()                   // abs(r)
r.array().abs()                // abs(r)
r.cwiseabs2()                  // abs(r.^2)
r.array().abs2()               // abs(r.^2)
(r.array() < s).select(p,q );  // (r < s ? p : q)
r = (q.array()==0).select(p,r) // r(q==0) = p(q==0)
r = p.unaryexpr(ptr_fun(func)) // r = arrayfun(func, p)   // with: scalar func(const scalar &x);

求最小最大值、迹等

// reductions.
int r, c;
// eigen                  // matlab
r.mincoeff()              // min(r(:))
r.maxcoeff()              // max(r(:))
s = r.mincoeff(&r, &c)    // [s, i] = min(r(:)); [r, c] = ind2sub(size(r), i);
s = r.maxcoeff(&r, &c)    // [s, i] = max(r(:)); [r, c] = ind2sub(size(r), i);
r.sum()                   // sum(r(:))
r.colwise().sum()         // sum(r)
r.rowwise().sum()         // sum(r, 2) or sum(r')'
r.prod()                  // prod(r(:))
r.colwise().prod()        // prod(r)
r.rowwise().prod()        // prod(r, 2) or prod(r')'
r.trace()                 // trace(r)
r.all()                   // all(r(:))
r.colwise().all()         // all(r)
r.rowwise().all()         // all(r, 2)
r.any()                   // any(r(:))
r.colwise().any()         // any(r)
r.rowwise().any()         // any(r, 2)

点乘等

// dot products, norms, etc.
// eigen                  // matlab
x.norm()                  // norm(x).    note that norm(r) doesn't work in eigen.
x.squarednorm()           // dot(x, x)   note the equivalence is not true for complex
x.dot(y)                  // dot(x, y)
x.cross(y)                // cross(x, y) requires #include <eigen/geometry>

特征值与特征向量

// eigenvalue problems
// eigen                          // matlab
a.eigenvalues();                  // eig(a);
eigensolver<matrix3d> eig(a);     // [vec val] = eig(a)
eig.eigenvalues();                // diag(val)
eig.eigenvectors();               // vec
// for self-adjoint matrices use selfadjointeigensolver<>

形式转换

 type conversion
// eigen                  // matlab
a.cast<double>();         // double(a)
a.cast<float>();          // single(a)
a.cast<int>();            // int32(a)
a.real();                 // real(a)
a.imag();                 // imag(a)
// if the original type equals destination type, no work is done

矩阵初始化0

// note that for most operations eigen requires all operands to have the same type:
matrixxf f = matrixxf::zero(3,3);
a += f;                // illegal in eigen. in matlab a = a+f is allowed
a += f.cast<double>(); // f converted to double and then added (generally, conversion happens on-the-fly)

map等操作

// eigen can map existing memory into eigen matrices.
float array[3];
vector3f::map(array).fill(10);            // create a temporary map over array and sets entries to 10
int data[4] = {1, 2, 3, 4};
matrix2i mat2x2(data);                    // copies data into mat2x2
matrix2i::map(data) = 2*mat2x2;           // overwrite elements of data with 2*mat2x2
matrixxi::map(data, 2, 2) += mat2x2;      // adds mat2x2 to elements of data (alternative syntax if size is not know at compile time)

求解ax = b

// solve ax = b. result stored in x. matlab: x = a \ b.
x = a.ldlt().solve(b));  // a sym. p.s.d.    #include <eigen/cholesky>
x = a.llt() .solve(b));  // a sym. p.d.      #include <eigen/cholesky>
x = a.lu()  .solve(b));  // stable and fast. #include <eigen/lu>
x = a.qr()  .solve(b));  // no pivoting.     #include <eigen/qr>
x = a.svd() .solve(b));  // stable, slowest. #include <eigen/svd>
// .ldlt() -> .matrixl() and .matrixd()
// .llt()  -> .matrixl()
// .lu()   -> .matrixl() and .matrixu()
// .qr()   -> .matrixq() and .matrixr()
// .svd()  -> .matrixu(), .singularvalues(), and .matrixv()

以上就是r语言学习rcppeigen进行矩阵运算的详细内容，更多关于rcppeigen矩阵运算的资料请关注其它相关文章！