1 使用矩陣作為函數參數介紹
文章來源Writing Functions Taking %Eigen Types as Parameters
Eigen為了在函數中傳遞不同的類型使用了表達式模板技術。如果你傳遞一個表達式到函數時使用了Matrix作為參數，你的表達式會被隱含的作為Matrix模板傳遞給這個函數。這意味着你丟掉了表達式模板所帶來的好處，這有如下兩個缺點：
- 評估進模板的參數可能是無用或者無效的。
- 這只允許函數從表達式中讀取而無法寫入。
幸運的是這些表達式類型都有一個共同特點就是他們都繼承於幾個通用的、模板化的基類。通過讓你的函數使用這些基類作為模板參數，你可以讓它們很好的和表達式模板結合。
2 一些簡單的例子
這一段將會提供一些簡單的例子用來說明Eigen中提供的幾種不同類型對象。在開始介紹這些例子之前，我們首先需要概括如何使用這些基類（當然也可以查看官方的手冊類的層次）。
- MatrixBase: 它是所有稠密矩陣表達式的基類(和array表達式、稀疏矩陣表達式、特殊矩陣表達式相對應)。在函數中使用它作為參數意味着它只能傳遞稠密矩陣。
- ArrayBase: 它是稠密array表達式的基類（和matrix表達式相對應等）。在函數中使用它作為參數意味着它只能傳遞array.
- DenseBase: 它是所有密集表達式的基類，也就是說它包含了MatrixBase和ArrayBase.它可以作為函數參數傳遞matrices和arrays.
- EigenBase: 它是所有矩陣表達式統一的基類，它可以被應用在矩陣或者arrays，比如特殊矩陣（比如對角矩陣、轉置矩陣等）.它被應用在函數中作為參數意味着所有這種類型的參數都可以被傳遞。
2.1 EigenBase使用例程
打印Eigen中通用對象的維數。它可能是任何形式的矩陣表達式，任何稠密或稀疏或array的形式。
#include "Eigen/Core"
#include "iostream"
using namespace Eigen;
using namespace std;
template <typename Derived>
void print_size(const EigenBase<Derived>& b)
{
cout<<"size(rows,cols):"<<b.size()<<" ("<<b.rows()
<<","<<b.cols()<<")"<<std::endl;
}
int main()
{
Vector3f v;
print_size(v);
//v.asDiagonal()returns a 3x3 diagonal matrix pseudo-expression
print_size(v.asDiagonal());
}
執行結果如下：
size(rows,cols):3 (3,1)
size(rows,cols):9 (3,3)
2.2 DenseBase使用示例
打印一個稠密表達式的子塊，接收任何稠密matrix和array表達式。但是稀疏類型和特殊類型的matrix類(如DiagonalMatrix)不能接收。
#include "Eigen/Core"
#include "iostream"
using namespace Eigen;
using namespace std;
template<typename Derived>
void print_block(const DenseBase<Derived>& b,int x,int y,int r,int c)
{
cout<<"block:\n"<<b.block(x,y,r,c)<<endl;
}
int main()
{
Matrix3f v;
v<< 1,2,3,
4,5,6,
7,8,9;
print_block(v,1,1,2,2);
}
執行結果如下：
block:
5 6
8 9
2.3 ArrayBase使用示例
打印array或其表達式中的最大元素。
#include "Eigen/Core"
#include "iostream"
using namespace Eigen;
using namespace std;
template <typename Derived>
void PrintMaxCoeff(const ArrayBase<Derived>& a)
{
cout<<"max: "<<a.maxCoeff()<<endl;
}
int main()
{
Array33f v;
v<< 1,2,3,
4,5,6,
7,8,9;
PrintMaxCoeff(v);
}
輸出結果：
max: 9
2.4 MatrixBase使用示例
打印給定矩陣或矩陣表達式的條件數倒數(inverse condition number).
#include "Eigen/SVD"
#include "iostream"
using namespace Eigen;
using namespace std;
template <typename Derived>
void print_inv_cond(const MatrixBase<Derived>& a)
{
const typename JacobiSVD<typename Derived::PlainObject>::SingularValuesType&
sing_vals = a.jacobiSvd().singularValues();
std::cout << "inv cond: " << sing_vals(sing_vals.size()-1) / sing_vals(0) << std::endl;
}
int main()
{
Matrix3f m;
m<< 1,2,3,
4,5,6,
7,8,9;
print_inv_cond(m);
}
執行結果：
inv cond: 1.11994e-008
2.5 多模板參數使用示例
計算兩點之間的歐式距離。
#include "Eigen/Core"
#include "iostream"
using namespace Eigen;
using namespace std;
template <typename Derived1,typename Derived2>
typename Derived1::Scalar SquareDist(const MatrixBase<Derived1>& p1,
const MatrixBase<Derived2>& p2)
{
return (p1-p2).squaredNorm();
}
int main()
{
Vector3f p1;
Vector3f p2;
p1<<1,2,3;
p2<<4,5,6;
cout<<"p1和p2之間的歐式距離為: "<<SquareDist(p1,p2)<<endl;
}
執行結果如下：
p1和p2之間的歐式距離為： 27
上面的幾個簡單的例子給讀者留下一個第一印象——如何用一個普通的和不變的Matrix或Array作為函數的參數。也故意給讀者留下了一個印象，在函數參數通用基類的最佳選項是什么。在下面一節的例子中我們會着眼於更加細節的地方以及實現它的不同方式，同時討論各種實現過程中的存在的問題和優勢。對於下面的討論，Matrix和Array以及MatrixBase和ArrayBase可以被交換並且所有的參數仍然成立。
3 如何寫一個通用的但是不是模板的函數？
在之前的所有例子中，函數都被寫成了模板的形式。這種方法允許我們寫一下非常通用的代碼，但是如果能寫出一種函數不使用模板同時還能保證一定程度的通用性這會是我們更加滿意，因為這可以避免很多愚蠢的參數拷貝。這有一個典型的例子來寫一個函數能同時接受MatrixXf和MatrixXf的塊。這正是Ref類的用途。例程如下：
#include "Eigen/SVD"
#include "iostream"
using namespace Eigen;
using namespace std;
float InvCond(const Ref<const MatrixXf>& a)
{
const VectorXf sing_vals = a.jacobiSvd().singularValues();
return sing_vals(sing_vals.size() - 1)/sing_vals(0);
}
int main()
{
Matrix4f m = Matrix4f::Random();
cout<<"matrix m:"<<endl<<m<<endl<<endl;
cout<<"InvCond(m):"<<InvCond(m)<<endl;
cout<<"InvCond(m(1:3,1:3)): "<<InvCond(m.topLeftCorner(3,3))<<endl;
cout<<"InvCond(m+1): "<<InvCond(m+Matrix4f::Identity())<<endl;
}
執行結果如下：
matrix m:
-0.997497 0.170019 0.64568 0.421003
0.127171 -0.0402539 0.49321 0.0270699
-0.613392 -0.299417 -0.651784 -0.39201
0.617481 0.791925 0.717887 -0.970031
InvCond(m):0.150053
InvCond(m(1:3,1:3)): 0.163871
InvCond(m+1): 0.202112
上面的例子中前兩次調用InvCond並沒有產生類的拷貝，因為參數的存儲布局(memory layout) 和 Ref<MatrixXf>的存儲布局相匹配。但是，最后一次調用時，我們的表達式會自動的被計算同時將值存入到一個由Ref<>對象決定的MatrixXf臨時變量中。
一個Ref對象同時也是可以被寫入的。下面這個例子是計算兩個矩陣的協方差，其中每一列是一個隨機變量，每一行是一組測量數據。
#include "Eigen/SVD"
#include "iostream"
using namespace Eigen;
using namespace std;
void cov(const Ref<const MatrixXf> x,const Ref<const MatrixXf> y,Ref<MatrixXf> C)
{
const float num_observations = static_cast<float>(x.rows());
const RowVectorXf x_mean = x.colwise().sum()/num_observations;
const RowVectorXf y_mean = y.colwise().sum()/num_observations;
C = (x.rowwise() - x_mean).transpose()*(y.rowwise() - y_mean)/num_observations;
}
int main()
{
Matrix<float,4,3> m,n;
MatrixXf conMn;
m<< 1,2,3,
4,5,6,
7,8,9,
10,11,12;
n = m;
conMn.resize(3,3); //在函數內部不允許修改大小,必須在調用之前給出大小
cov(m,n,conMn);
cout<<conMn<<endl;
}
執行結果如下：
11.25 11.25 11.25
11.25 11.25 11.25
11.25 11.25 11.25
當然我們也可以使用矩陣的子塊作為參數conv(m.leftCols<3>(),n.leftCols<3>(),conMn.topLeftCorner<3,3>());.
Ref<>類還有其他兩個模板參數運行控制存儲布局(memory layout)以便接收參數的同時不進行任何拷貝。查看Ref文檔以便了解更多細節。
4 什么情況下函數傳遞一個普通的Matrix或Array參數能工作？
即不使用模板也不使用Ref類，上面例子的一個狹窄實現可能是下面這個樣子。
MatrixXf cov(const MatrixXf& x, const MatrixXf& y)
{
const float num_observations = static_cast<float>(x.rows());
const RowVectorXf x_mean = x.colwise().sum() / num_observations;
const RowVectorXf y_mean = y.colwise().sum() / num_observations;
return (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}
和一些人的第一直覺相反，這個實現能很好的工作除非你需要一個通用的實現讓double類型的matrix也工作或者除非你不關心模板對象。為什么呢？臨時變量放在哪呢？下面這段代碼又是怎么運行的呢？
MatrixXf x,y,z;
MatrixXf C = conv(x,y+z);
在這個特殊情況下，這個例子能正常工作是因為所有的參數都聲明為const reference.編譯器創建一個臨時變量並且計算並保存x+z.一旦程序運行完成，臨時變量就會被釋放掉同時將結果返回給C。
注意：在使用了臨時變量的代價下，函數對Matrix(或Array)使用const引用可以運行表達式。
5 什么情況下函數使用普通的Matrix或Array作為參數會失敗？
這里我們對上面的例子進行小小的修改，我們不想讓函數有返回值而是通過傳遞一個額外的非const參數來保存結果，一個狹窄的實現可能是下面這個樣子：
//Note: This code is flawed!
void cov(const MatrixXf& x, const MatrixXf& y, MatrixXf& C)
{
const float num_observations = static_cast<float>(x.rows());
const RowVectorXf x_mean = x.colwise().sum() / num_observations;
const RowVectorXf y_mean = y.colwise().sum() / num_observations;
C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}
接下來我們嘗試執行下面的語句：
MatrixXf C = MatrixXf::Zero(3,6);
cov(x,y, C.block(0,0,3,3));
這時會編譯失敗，因為不可能轉換一個非const的MatrixXf&表達式給MatrixXf::block().這是因為編譯器為了保護你的程序安全會避免你向一個臨時變量寫入結果。這里雖然我們很想向::block()里寫入返回結果。如何解決這個問題呢？
這里優先選擇的一個解決方案需要一點駭客手段。需要傳遞一個const引用這個matrix然后在內部使用去const的手段去掉const。在C98編譯器下一個正確的實現如下所示：
template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> const & C)
{
typedef typename Derived::Scalar Scalar;
typedef typename internal::plain_row_type<Derived>::type RowVectorType;
const Scalar num_observations = static_cast<Scalar>(x.rows());
const RowVectorType x_mean = x.colwise().sum() / num_observations;
const RowVectorType y_mean = y.colwise().sum() / num_observations;
const_cast< MatrixBase<OtherDerived>& >(C) =
(x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}
上面這個實現中不但允許向臨時表達式里寫入數據，同時也允許使用任意浮點類型矩陣作為參數的函數。
注意：去const的駭客手段只能在模板函數中工作，它不能在MatrixXf實現中工作因為不可能轉換一個Block表達式到MatrixXf的引用。
6 在一般的實現中如何修改matrices的大小？
到這里我們所有的工作都做完了嗎？前面我們例子中了解到在函數調用時不能改變matrix的大小，然后對於一般的應用，我們更習慣使用如下的方式調用前面的協方差函數。
MatrixXf x = MatrixXf::Random(100,3);
MatrixXf y = MatrixXf::Random(100,3);
MatrixXf C;
cov(x,y,C);
上面情況會因為C的維數改變而出現錯誤，在我們使用implementation作為MatrixBase的參數時。通常，Eigen支持自動的維數轉換但不能通過表達式這樣做。為什么應該允許改變一個矩陣塊的維數？它是一個到子矩陣的引用，我們真的不應該去修改它的維數（試想一下，對於一個2x2的矩陣，如果我們將它子矩陣(1,1)變成2x2那結果還能是矩陣嗎？）。那么如何不再MatrixBase的基礎上改變維數呢？在這個實現中辦法是改變derived object的維數。
#include "Eigen/SVD"
#include "iostream"
using namespace Eigen;
using namespace std;
template <typename Derived,typename OtherDerived>
void cov(const MatrixBase<Derived>& x,const MatrixBase<Derived>& y,
MatrixBase<OtherDerived> const& C_)
{
typedef typename Derived::Scalar Scalar;
typedef typename internal::plain_row_type<Derived>::type RowVectorType;
const Scalar num_observations = static_cast<Scalar>(x.rows());
const RowVectorType x_mean = x.colwise().sum()/num_observations;
const RowVectorType y_mean = y.colwise().sum()/num_observations;
MatrixBase<OtherDerived>& C = const_cast<MatrixBase<OtherDerived>&>(C_);
C.derived().resize(x.cols(),x.cols()); //resize the derived object
C = (x.rowwise() - x_mean).transpose()*(y.rowwise() - y_mean)/num_observations;
}
int main()
{
Matrix<float,4,3> m,n;
MatrixXf conMn;
m<< 1,2,3,
4,5,6,
7,8,9,
10,11,12;
n = m;
// conMn.resize(3,3); //在函數內部不允許修改大小,必須在調用之前給出大小
cov(m,n,conMn);
cout<<conMn<<endl;
}
執行結果如下：
11.25 11.25 11.25
11.25 11.25 11.25
11.25 11.25 11.25
上面這個實現中不但可以讓表達式作為參數，也可以讓有錯誤維數的矩陣作為參數。
上面的討論對Matrix,Array,MatrixBase,ArrayBase進行交換所有結論依然成立。
7 總結
- 實現使用不可寫(const referenced)對象作為參數的函數不是一個大問題，同時也不會在編譯和運行時引起任何問題。但是，一個狹窄的實現可能引入不必要中間對象。通過模板使用到MatrixBase或ArrayBase(const)引用可以避免上述情況。
- 函數使用了一個可寫參數必須使用const引用然后在函數體中去除const。
- 函數使用了MatrixBase(或ArrayBase)對象，如果有可能改變它們的維數，必須在derived類上調用resize()然后返回derived().
ps: 這篇隨便是邊看官方手冊邊實驗邊記錄的，有些地方當時也不太懂，另外限於本人中文表達能力和英文閱讀能力都有限，可能有些地方翻譯錯誤，建議英語水平高的同學還是直接看官網手冊比較好。