ground/gcs/src/libs/eigen/doc/TutorialMatrixArithmetic.dox

   1 namespace Eigen {
   2
   3 /** \eigenManualPage TutorialMatrixArithmetic Matrix and vector arithmetic
   4
   5 This page aims to provide an overview and some details on how to perform arithmetic
   6 between matrices, vectors and scalars with Eigen.
   7
   8 \eigenAutoToc
   9
  10 \section TutorialArithmeticIntroduction Introduction
  11
  12 Eigen offers matrix/vector arithmetic operations either through overloads of common C++ arithmetic operators such as +, -, *,
  13 or through special methods such as dot(), cross(), etc.
  14 For the Matrix class (matrices and vectors), operators are only overloaded to support
  15 linear-algebraic operations. For example, \c matrix1 \c * \c matrix2 means matrix-matrix product,
  16 and \c vector \c + \c scalar is just not allowed. If you want to perform all kinds of array operations,
  17 not linear algebra, see the \ref TutorialArrayClass "next page".
  18
  19 \section TutorialArithmeticAddSub Addition and subtraction
  20
  21 The left hand side and right hand side must, of course, have the same numbers of rows and of columns. They must
  22 also have the same \c Scalar type, as Eigen doesn't do automatic type promotion. The operators at hand here are:
  23 \li binary operator + as in \c a+b
  24 \li binary operator - as in \c a-b
  25 \li unary operator - as in \c -a
  26 \li compound operator += as in \c a+=b
  27 \li compound operator -= as in \c a-=b
  28
  29 <table class="example">
  30 <tr><th>Example:</th><th>Output:</th></tr>
  31 <tr><td>
  32 \include tut_arithmetic_add_sub.cpp
  33 </td>
  34 <td>
  35 \verbinclude tut_arithmetic_add_sub.out
  36 </td></tr></table>
  37
  38 \section TutorialArithmeticScalarMulDiv Scalar multiplication and division
  39
  40 Multiplication and division by a scalar is very simple too. The operators at hand here are:
  41 \li binary operator * as in \c matrix*scalar
  42 \li binary operator * as in \c scalar*matrix
  43 \li binary operator / as in \c matrix/scalar
  44 \li compound operator *= as in \c matrix*=scalar
  45 \li compound operator /= as in \c matrix/=scalar
  46
  47 <table class="example">
  48 <tr><th>Example:</th><th>Output:</th></tr>
  49 <tr><td>
  50 \include tut_arithmetic_scalar_mul_div.cpp
  51 </td>
  52 <td>
  53 \verbinclude tut_arithmetic_scalar_mul_div.out
  54 </td></tr></table>
  55
  56
  57 \section TutorialArithmeticMentionXprTemplates A note about expression templates
  58
  59 This is an advanced topic that we explain on \ref TopicEigenExpressionTemplates "this page",
  60 but it is useful to just mention it now. In Eigen, arithmetic operators such as \c operator+ don't
  61 perform any computation by themselves, they just return an "expression object" describing the computation to be
  62 performed. The actual computation happens later, when the whole expression is evaluated, typically in \c operator=.
  63 While this might sound heavy, any modern optimizing compiler is able to optimize away that abstraction and
  64 the result is perfectly optimized code. For example, when you do:
  65 \code
  66 VectorXf a(50), b(50), c(50), d(50);
  67 ...
  68 a = 3*b + 4*c + 5*d;
  69 \endcode
  70 Eigen compiles it to just one for loop, so that the arrays are traversed only once. Simplifying (e.g. ignoring
  71 SIMD optimizations), this loop looks like this:
  72 \code
  73 for(int i = 0; i < 50; ++i)
  74   a[i] = 3*b[i] + 4*c[i] + 5*d[i];
  75 \endcode
  76 Thus, you should not be afraid of using relatively large arithmetic expressions with Eigen: it only gives Eigen
  77 more opportunities for optimization.
  78
  79 \section TutorialArithmeticTranspose Transposition and conjugation
  80
  81 The transpose \f$ a^T \f$, conjugate \f$ \bar{a} \f$, and adjoint (i.e., conjugate transpose) \f$ a^* \f$ of a matrix or vector \f$ a \f$ are obtained by the member functions \link DenseBase::transpose() transpose()\endlink, \link MatrixBase::conjugate() conjugate()\endlink, and \link MatrixBase::adjoint() adjoint()\endlink, respectively.
  82
  83 <table class="example">
  84 <tr><th>Example:</th><th>Output:</th></tr>
  85 <tr><td>
  86 \include tut_arithmetic_transpose_conjugate.cpp
  87 </td>
  88 <td>
  89 \verbinclude tut_arithmetic_transpose_conjugate.out
  90 </td></tr></table>
  91
  92 For real matrices, \c conjugate() is a no-operation, and so \c adjoint() is equivalent to \c transpose().
  93
  94 As for basic arithmetic operators, \c transpose() and \c adjoint() simply return a proxy object without doing the actual transposition. If you do <tt>b = a.transpose()</tt>, then the transpose is evaluated at the same time as the result is written into \c b. However, there is a complication here. If you do <tt>a = a.transpose()</tt>, then Eigen starts writing the result into \c a before the evaluation of the transpose is finished. Therefore, the instruction <tt>a = a.transpose()</tt> does not replace \c a with its transpose, as one would expect:
  95 <table class="example">
  96 <tr><th>Example:</th><th>Output:</th></tr>
  97 <tr><td>
  98 \include tut_arithmetic_transpose_aliasing.cpp
  99 </td>
 100 <td>
 101 \verbinclude tut_arithmetic_transpose_aliasing.out
 102 </td></tr></table>
 103 This is the so-called \ref TopicAliasing "aliasing issue". In "debug mode", i.e., when \ref TopicAssertions "assertions" have not been disabled, such common pitfalls are automatically detected.
 104
 105 For \em in-place transposition, as for instance in <tt>a = a.transpose()</tt>, simply use the \link DenseBase::transposeInPlace() transposeInPlace()\endlink  function:
 106 <table class="example">
 107 <tr><th>Example:</th><th>Output:</th></tr>
 108 <tr><td>
 109 \include tut_arithmetic_transpose_inplace.cpp
 110 </td>
 111 <td>
 112 \verbinclude tut_arithmetic_transpose_inplace.out
 113 </td></tr></table>
 114 There is also the \link MatrixBase::adjointInPlace() adjointInPlace()\endlink function for complex matrices.
 115
 116 \section TutorialArithmeticMatrixMul Matrix-matrix and matrix-vector multiplication
 117
 118 Matrix-matrix multiplication is again done with \c operator*. Since vectors are a special
 119 case of matrices, they are implicitly handled there too, so matrix-vector product is really just a special
 120 case of matrix-matrix product, and so is vector-vector outer product. Thus, all these cases are handled by just
 121 two operators:
 122 \li binary operator * as in \c a*b
 123 \li compound operator *= as in \c a*=b (this multiplies on the right: \c a*=b is equivalent to <tt>a = a*b</tt>)
 124
 125 <table class="example">
 126 <tr><th>Example:</th><th>Output:</th></tr>
 127 <tr><td>
 128 \include tut_arithmetic_matrix_mul.cpp
 129 </td>
 130 <td>
 131 \verbinclude tut_arithmetic_matrix_mul.out
 132 </td></tr></table>
 133
 134 Note: if you read the above paragraph on expression templates and are worried that doing \c m=m*m might cause
 135 aliasing issues, be reassured for now: Eigen treats matrix multiplication as a special case and takes care of
 136 introducing a temporary here, so it will compile \c m=m*m as:
 137 \code
 138 tmp = m*m;
 139 m = tmp;
 140 \endcode
 141 If you know your matrix product can be safely evaluated into the destination matrix without aliasing issue, then you can use the \link MatrixBase::noalias() noalias()\endlink function to avoid the temporary, e.g.:
 142 \code
 143 c.noalias() += a * b;
 144 \endcode
 145 For more details on this topic, see the page on \ref TopicAliasing "aliasing".
 146
 147 \b Note: for BLAS users worried about performance, expressions such as <tt>c.noalias() -= 2 * a.adjoint() * b;</tt> are fully optimized and trigger a single gemm-like function call.
 148
 149 \section TutorialArithmeticDotAndCross Dot product and cross product
 150
 151 For dot product and cross product, you need the \link MatrixBase::dot() dot()\endlink and \link MatrixBase::cross() cross()\endlink methods. Of course, the dot product can also be obtained as a 1x1 matrix as u.adjoint()*v.
 152 <table class="example">
 153 <tr><th>Example:</th><th>Output:</th></tr>
 154 <tr><td>
 155 \include tut_arithmetic_dot_cross.cpp
 156 </td>
 157 <td>
 158 \verbinclude tut_arithmetic_dot_cross.out
 159 </td></tr></table>
 160
 161 Remember that cross product is only for vectors of size 3. Dot product is for vectors of any sizes.
 162 When using complex numbers, Eigen's dot product is conjugate-linear in the first variable and linear in the
 163 second variable.
 164
 165 \section TutorialArithmeticRedux Basic arithmetic reduction operations
 166 Eigen also provides some reduction operations to reduce a given matrix or vector to a single value such as the sum (computed by \link DenseBase::sum() sum()\endlink), product (\link DenseBase::prod() prod()\endlink), or the maximum (\link DenseBase::maxCoeff() maxCoeff()\endlink) and minimum (\link DenseBase::minCoeff() minCoeff()\endlink) of all its coefficients.
 167
 168 <table class="example">
 169 <tr><th>Example:</th><th>Output:</th></tr>
 170 <tr><td>
 171 \include tut_arithmetic_redux_basic.cpp
 172 </td>
 173 <td>
 174 \verbinclude tut_arithmetic_redux_basic.out
 175 </td></tr></table>
 176
 177 The \em trace of a matrix, as returned by the function \link MatrixBase::trace() trace()\endlink, is the sum of the diagonal coefficients and can also be computed as efficiently using <tt>a.diagonal().sum()</tt>, as we will see later on.
 178
 179 There also exist variants of the \c minCoeff and \c maxCoeff functions returning the coordinates of the respective coefficient via the arguments:
 180
 181 <table class="example">
 182 <tr><th>Example:</th><th>Output:</th></tr>
 183 <tr><td>
 184 \include tut_arithmetic_redux_minmax.cpp
 185 </td>
 186 <td>
 187 \verbinclude tut_arithmetic_redux_minmax.out
 188 </td></tr></table>
 189
 190
 191 \section TutorialArithmeticValidity Validity of operations
 192 Eigen checks the validity of the operations that you perform. When possible,
 193 it checks them at compile time, producing compilation errors. These error messages can be long and ugly,
 194 but Eigen writes the important message in UPPERCASE_LETTERS_SO_IT_STANDS_OUT. For example:
 195 \code
 196   Matrix3f m;
 197   Vector4f v;
 198   v = m*v;      // Compile-time error: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES
 199 \endcode
 200
 201 Of course, in many cases, for example when checking dynamic sizes, the check cannot be performed at compile time.
 202 Eigen then uses runtime assertions. This means that the program will abort with an error message when executing an illegal operation if it is run in "debug mode", and it will probably crash if assertions are turned off.
 203
 204 \code
 205   MatrixXf m(3,3);
 206   VectorXf v(4);
 207   v = m * v; // Run-time assertion failure here: "invalid matrix product"
 208 \endcode
 209
 210 For more details on this topic, see \ref TopicAssertions "this page".
 211
 212 */
 213
 214 }