1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
|
namespace Eigen {
/** \page TopicFunctionTakingEigenTypes Writing Functions Taking Eigen Types as Parameters
Eigen's use of expression templates results in almost all operations being objects of a different type. The effect of this design is that without writing template functions it is not possible to provide \e generically applicable methods taking simple references to Matrices or Arrays.
This page explains how to write generic functions taking Eigen types and potential pitfalls one should be aware of in order to write safe and efficient code.
<b>Table of contents</b>
- \ref TopicFirstExamples
- \ref TopicSimpleFunctionsWorking
- \ref TopicSimpleFunctionsFailing
- \ref TopicResizingInGenericImplementations
- \ref TopicSummary
\section TopicFirstExamples Some First Examples
This section will provide simple examples for different types of objects Eigen is offering. Before starting with the actual examples, we need to recapitulate which base objects we can work with (see also \ref TopicClassHierarchy).
- MatrixBase: A dense matrix-expression which is any object that can be used for matrix arithmetic or on matrix decompositions.
- ArrayBase: A dense array-expressions which is any object that can be used in basic and component-wise arithmetic expressions.
- DenseBase: The base class for \c MatrixBase and \c ArrayBase. It can be used in functions that are meant to work on both matrices and arrays.
- EigenBase: The base unifying all types of objects that can be evaluated into dense matrices or arrays.
<b> %EigenBase Example </b><br/><br/>
Prints the dimensions of the most generic object present in Eigen. It coulde be any matrix expressions, any dense or sparse matrix and any array.
\code
template <typename Derived>
void print_size(const EigenBase<Derived>& b)
{
std::cout << "size (rows, cols): " << b.size() << " (" << b.rows() << ", " << b.cols() << ")" << std::endl;
}
\endcode
<b> %DenseBase Example </b><br/><br/>
Prints a sub-block of the dense expressions. Could all of the above but no sparse objects.
\code
template <typename Derived>
void print_block(const DenseBase<Derived>& b, int x, int y, int r, int c)
{
std::cout << "block: " << b.block(x,y,r,c) << std::endl;
}
\endcode
<b> %ArrayBase Example </b><br/><br/>
Prints the maximum coefficient from both arrays or array-expressions.
\code
template <typename Derived>
void print_max(const ArrayBase<Derived>& a, const ArrayBase<Derived>& b)
{
std::cout << "max: " << (a.max(b)).maxCoeff() << std::endl;
}
\endcode
<b> %MatrixBase Example </b><br/><br/>
Prints the inverse condition number of the given matrix or matrix-expression.
\code
template <typename Derived>
void print_inv_cond(const MatrixBase<Derived>& a)
{
const typename SVD<typename Derived::PlainObject>::SingularValuesType& sing_vals = a.svd().singularValues();
std::cout << "inv cond: " << sing_vals(sing_vals.size()-1) / sing_vals(0) << std::endl;
}
\endcode
These examples are just intended to give the reader a first impression of how simple functions taking working on constant objects can be written. They are also intended to give the reader an idea about the most common base classes being the optimal candidates for functions. In the next section we will look in more detail at an example and the different ways it can be implemented, while discussing each implementation's problems and advantages. For the discussion below, Matrix and Array as well as MatrixBase and ArrayBase can be exchanged and all arguments still hold.
\section TopicSimpleFunctionsWorking In which cases do simple functions work?
Let's assume one wants to write a function computing the covariance matrix of two input matrices where each row is an observation. The implementation of this function might look like this
\code
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;
}
\endcode
and contrary to what one might think at first, this implementation is fine unless you require a genric implementation that works with double matrices too and unless you do not care about temporary objects. Why is that the case? Where are temporaries involved? How can code as given below compile?
\code
MatrixXf x,y,z;
MatrixXf C = cov(x,y+z);
\endcode
In this special case, the example is fine and will be working because both parameters are declared as \e const references. The compiler creates a temporary and evaluates the expression x+z into this temporary. Once the function is processed, the temporary is released and the result is assigned to C.
\b Note: Functions taking \e const references to Matrix (or Array) can process expressions at the cost of temporaries.
\section TopicSimpleFunctionsFailing In which cases do simple functions fail?
Here, we consider a slightly modified version of the function given above. This time, we do not want to return the result but pass an additional non-const paramter which allows us to store the result. A first naive implementation might look as follows.
\code
// 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;
}
\endcode
When trying to execute the following code
\code
MatrixXf C = MatrixXf::Zero(3,6);
cov(x,y, C.block(0,0,3,3));
\endcode
the compiler will fail, because it is not possible to convert the expression returned by \c MatrixXf::block() in a non-const \c MatrixXf&. This is the case because the compiler wants to protect you from writing your result to a temporary object. In this special case this protection is not intended -- we want to write to a temporary object. So how can we overcome this problem? There are two possible solutions depending on the type of compiler you are using.
<b>Solution A)</b>
Assuming you are using a compiler following the C98 standard, the only thing you can do is to use a little \em hack. You need to pass a const reference and internally the constness needs to be cast away. The correct implementation for C98 compliant compilers would be
\code
template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY C)
{
typedef typename ei_traits<Derived>::Scalar Scalar;
typedef typename ei_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;
}
\endcode
The implementation above does now not only work with temporary expressions but it also allows to use the function with matrices of arbitrary floating point scalar types.
\b Note: The const cast hack will only work with templated functions. It will not work with the MatrixXf implementation because it is not possible to cast a Block expression to a Matrix reference!
<b>Solution B)</b>
In the next solution we are going to utilize a new feature introduced with C++0x compliant compilers -- so called rvalue references. Rvalue references allow to explicitly tell the compiler that the object we are going to pass to a function is a temporary object that is writeable. The C++0x compliant implementation of the covariance function will be
\code
template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived>&& C)
{
typedef typename ei_traits<Derived>::Scalar Scalar;
typedef typename ei_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;
C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}
\endcode
\b Note: Currently, this fails to work on Visual Studio 2010 even though the compiler is assumed to be supporting rvalue references. So stick to the EIGEN_REF_TO_TEMPORARY \em hack.
\section TopicResizingInGenericImplementations How to resize matrices in generic implementations?
One might think we are done now, right? This is not completely true because in order for our covariance function to be generically applicable, we want the follwing code to work
\code
MatrixXf x = MatrixXf::Random(100,3);
MatrixXf y = MatrixXf::Random(100,3);
MatrixXf C;
cov(x, y, C);
\endcode
This is not the case anymore, when we are using an implementation taking MatrixBase as a parameter. In general, Eigen supports automatic resizing but it is not possible to do so on expressions. Why should resizing of a matrix Block be allowed? It is a reference to a sub-matrix and we definitely don't want to resize that. So how can we incorporate resizing if we cannot resize on MatrixBase? The solution is to resize the derived object as in this implementation.
\code
template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY C_)
{
typedef typename ei_traits<Derived>::Scalar Scalar;
typedef typename ei_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;
}
\endcode
This implementation is now working for paramters being expressions and for parameters being matrices and having the wrong size. Resizing the expressions does not do any harm in this case unless they actually require resizing. That means, passing an expression with the wrong dimensions will result in a run-time error (in debug mode only) while passing expressions of the correct size will just work fine.
\b Note: In the above discussion the terms Matrix and Array and MatrixBase and ArrayBase can be exchanged and all arguments still hold.
\section TopicSummary Summary
- To summarize, the implementation of functions taking non-writable (const referenced) objects is not a big issue and does not lead to problematic situations in terms of compiling and running your program. However, such an implementation is likely to introduce unnecessary temporary objects in your code. In case you care about optimal run-time performance, use (const) references to MatrixBase or ArrayBase, i.e. write templated functions.
- Writing functions taking parameters that are writable (non-const) requires to implement the functions such that they take references to MatrixBase or ArrayBase. Furthermore, writable objects require to use Eigen's macro EIGEN_REF_TO_TEMPORARY and casting away constness within the function.
- Regarding resizing objects, in particular in functions that take as parameters MatrixBase (or ArrayBase), the actual resizing has to be performed on the derived class.
*/
}
|
