eigen.rst revision 12391
12810SN/AEigen 210764Sandreas.hansson@arm.com##### 39347SAndreas.Sandberg@arm.com 49347SAndreas.Sandberg@arm.com`Eigen <http://eigen.tuxfamily.org>`_ is C++ header-based library for dense and 59347SAndreas.Sandberg@arm.comsparse linear algebra. Due to its popularity and widespread adoption, pybind11 69347SAndreas.Sandberg@arm.comprovides transparent conversion and limited mapping support between Eigen and 79347SAndreas.Sandberg@arm.comScientific Python linear algebra data types. 89347SAndreas.Sandberg@arm.com 99347SAndreas.Sandberg@arm.comTo enable the built-in Eigen support you must include the optional header file 109347SAndreas.Sandberg@arm.com:file:`pybind11/eigen.h`. 119347SAndreas.Sandberg@arm.com 129347SAndreas.Sandberg@arm.comPass-by-value 139347SAndreas.Sandberg@arm.com============= 142810SN/A 152810SN/AWhen binding a function with ordinary Eigen dense object arguments (for 162810SN/Aexample, ``Eigen::MatrixXd``), pybind11 will accept any input value that is 172810SN/Aalready (or convertible to) a ``numpy.ndarray`` with dimensions compatible with 182810SN/Athe Eigen type, copy its values into a temporary Eigen variable of the 192810SN/Aappropriate type, then call the function with this temporary variable. 202810SN/A 212810SN/ASparse matrices are similarly copied to or from 222810SN/A``scipy.sparse.csr_matrix``/``scipy.sparse.csc_matrix`` objects. 232810SN/A 242810SN/APass-by-reference 252810SN/A================= 262810SN/A 272810SN/AOne major limitation of the above is that every data conversion implicitly 282810SN/Ainvolves a copy, which can be both expensive (for large matrices) and disallows 292810SN/Abinding functions that change their (Matrix) arguments. Pybind11 allows you to 302810SN/Awork around this by using Eigen's ``Eigen::Ref<MatrixType>`` class much as you 312810SN/Awould when writing a function taking a generic type in Eigen itself (subject to 322810SN/Asome limitations discussed below). 332810SN/A 342810SN/AWhen calling a bound function accepting a ``Eigen::Ref<const MatrixType>`` 352810SN/Atype, pybind11 will attempt to avoid copying by using an ``Eigen::Map`` object 362810SN/Athat maps into the source ``numpy.ndarray`` data: this requires both that the 372810SN/Adata types are the same (e.g. ``dtype='float64'`` and ``MatrixType::Scalar`` is 382810SN/A``double``); and that the storage is layout compatible. The latter limitation 392810SN/Ais discussed in detail in the section below, and requires careful 402810SN/Aconsideration: by default, numpy matrices and eigen matrices are *not* storage 419347SAndreas.Sandberg@arm.comcompatible. 422810SN/A 432810SN/AIf the numpy matrix cannot be used as is (either because its types differ, e.g. 442810SN/Apassing an array of integers to an Eigen paramater requiring doubles, or 454626SN/Abecause the storage is incompatible), pybind11 makes a temporary copy and 462810SN/Apasses the copy instead. 472810SN/A 4810509SAli.Saidi@ARM.comWhen a bound function parameter is instead ``Eigen::Ref<MatrixType>`` (note the 495338Sstever@gmail.comlack of ``const``), pybind11 will only allow the function to be called if it 5010509SAli.Saidi@ARM.comcan be mapped *and* if the numpy array is writeable (that is 512810SN/A``a.flags.writeable`` is true). Any access (including modification) made to 522810SN/Athe passed variable will be transparently carried out directly on the 532810SN/A``numpy.ndarray``. 545314SN/A 5510622Smitch.hayenga@arm.comThis means you can can write code such as the following and have it work as 5610622Smitch.hayenga@arm.comexpected: 579725Sandreas.hansson@arm.com 5810622Smitch.hayenga@arm.com.. code-block:: cpp 5910622Smitch.hayenga@arm.com 6010622Smitch.hayenga@arm.com void scale_by_2(Eigen::Ref<Eigen::VectorXd> v) { 612810SN/A v *= 2; 624626SN/A } 634626SN/A 642810SN/ANote, however, that you will likely run into limitations due to numpy and 652810SN/AEigen's difference default storage order for data; see the below section on 662810SN/A:ref:`storage_orders` for details on how to bind code that won't run into such 672810SN/Alimitations. 684626SN/A 6910764Sandreas.hansson@arm.com.. note:: 702810SN/A 7110766Sandreas.hansson@arm.com Passing by reference is not supported for sparse types. 7210768Sandreas.hansson@arm.com 7310768Sandreas.hansson@arm.comReturning values to Python 7410768Sandreas.hansson@arm.com========================== 7510768Sandreas.hansson@arm.com 7610768Sandreas.hansson@arm.comWhen returning an ordinary dense Eigen matrix type to numpy (e.g. 7710768Sandreas.hansson@arm.com``Eigen::MatrixXd`` or ``Eigen::RowVectorXf``) pybind11 keeps the matrix and 7810768Sandreas.hansson@arm.comreturns a numpy array that directly references the Eigen matrix: no copy of the 792810SN/Adata is performed. The numpy array will have ``array.flags.owndata`` set to 802810SN/A``False`` to indicate that it does not own the data, and the lifetime of the 812810SN/Astored Eigen matrix will be tied to the returned ``array``. 822810SN/A 832810SN/AIf you bind a function with a non-reference, ``const`` return type (e.g. 842810SN/A``const Eigen::MatrixXd``), the same thing happens except that pybind11 also 852810SN/Asets the numpy array's ``writeable`` flag to false. 8610764Sandreas.hansson@arm.com 8710764Sandreas.hansson@arm.comIf you return an lvalue reference or pointer, the usual pybind11 rules apply, 882810SN/Aas dictated by the binding function's return value policy (see the 892810SN/Adocumentation on :ref:`return_value_policies` for full details). That means, 902810SN/Awithout an explicit return value policy, lvalue references will be copied and 912810SN/Apointers will be managed by pybind11. In order to avoid copying, you should 9210766Sandreas.hansson@arm.comexplictly specify an appropriate return value policy, as in the following 9310768Sandreas.hansson@arm.comexample: 9410768Sandreas.hansson@arm.com 952810SN/A.. code-block:: cpp 962810SN/A 972810SN/A class MyClass { 982810SN/A Eigen::MatrixXd big_mat = Eigen::MatrixXd::Zero(10000, 10000); 992810SN/A public: 1004920SN/A Eigen::MatrixXd &getMatrix() { return big_mat; } 1012810SN/A const Eigen::MatrixXd &viewMatrix() { return big_mat; } 1024920SN/A }; 1034920SN/A 1044920SN/A // Later, in binding code: 1054920SN/A py::class_<MyClass>(m, "MyClass") 1065314SN/A .def(py::init<>()) 10710766Sandreas.hansson@arm.com .def("copy_matrix", &MyClass::getMatrix) // Makes a copy! 10810764Sandreas.hansson@arm.com .def("get_matrix", &MyClass::getMatrix, py::return_value_policy::reference_internal) 1095314SN/A .def("view_matrix", &MyClass::viewMatrix, py::return_value_policy::reference_internal) 1104920SN/A ; 1114920SN/A 1124920SN/A.. code-block:: python 1135314SN/A 1144920SN/A a = MyClass() 1152810SN/A m = a.get_matrix() # flags.writeable = True, flags.owndata = False 1162810SN/A v = a.view_matrix() # flags.writeable = False, flags.owndata = False 1174920SN/A c = a.copy_matrix() # flags.writeable = True, flags.owndata = True 1184626SN/A # m[5,6] and v[5,6] refer to the same element, c[5,6] does not. 11910764Sandreas.hansson@arm.com 1202810SN/ANote in this example that ``py::return_value_policy::reference_internal`` is 12110766Sandreas.hansson@arm.comused to tie the life of the MyClass object to the life of the returned arrays. 12210764Sandreas.hansson@arm.com 12310764Sandreas.hansson@arm.comYou may also return an ``Eigen::Ref``, ``Eigen::Map`` or other map-like Eigen 1242810SN/Aobject (for example, the return value of ``matrix.block()`` and related 1252810SN/Amethods) that map into a dense Eigen type. When doing so, the default 1262810SN/Abehaviour of pybind11 is to simply reference the returned data: you must take 1272810SN/Acare to ensure that this data remains valid! You may ask pybind11 to 1282810SN/Aexplicitly *copy* such a return value by using the 1294666SN/A``py::return_value_policy::copy`` policy when binding the function. You may 1304666SN/Aalso use ``py::return_value_policy::reference_internal`` or a 1314666SN/A``py::keep_alive`` to ensure the data stays valid as long as the returned numpy 1324666SN/Aarray does. 1334871SN/A 1344666SN/AWhen returning such a reference of map, pybind11 additionally respects the 1354666SN/Areadonly-status of the returned value, marking the numpy array as non-writeable 1364666SN/Aif the reference or map was itself read-only. 13710766Sandreas.hansson@arm.com 1384871SN/A.. note:: 1394666SN/A 1404666SN/A Sparse types are always copied when returned. 1414666SN/A 1424666SN/A.. _storage_orders: 14310766Sandreas.hansson@arm.com 1444666SN/AStorage orders 1454666SN/A============== 1464666SN/A 1474626SN/APassing arguments via ``Eigen::Ref`` has some limitations that you must be 14810764Sandreas.hansson@arm.comaware of in order to effectively pass matrices by reference. First and 14910764Sandreas.hansson@arm.comforemost is that the default ``Eigen::Ref<MatrixType>`` class requires 1502810SN/Acontiguous storage along columns (for column-major types, the default in Eigen) 1513149SN/Aor rows if ``MatrixType`` is specifically an ``Eigen::RowMajor`` storage type. 1522810SN/AThe former, Eigen's default, is incompatible with ``numpy``'s default row-major 1532810SN/Astorage, and so you will not be able to pass numpy arrays to Eigen by reference 1542810SN/Awithout making one of two changes. 1552810SN/A 15610764Sandreas.hansson@arm.com(Note that this does not apply to vectors (or column or row matrices): for such 1572810SN/Atypes the "row-major" and "column-major" distinction is meaningless). 1584666SN/A 1592810SN/AThe first approach is to change the use of ``Eigen::Ref<MatrixType>`` to the 1602810SN/Amore general ``Eigen::Ref<MatrixType, 0, Eigen::Stride<Eigen::Dynamic, 1612810SN/AEigen::Dynamic>>`` (or similar type with a fully dynamic stride type in the 1622810SN/Athird template argument). Since this is a rather cumbersome type, pybind11 1632810SN/Aprovides a ``py::EigenDRef<MatrixType>`` type alias for your convenience (along 1642810SN/Awith EigenDMap for the equivalent Map, and EigenDStride for just the stride 1652810SN/Atype). 1664626SN/A 1672810SN/AThis type allows Eigen to map into any arbitrary storage order. This is not 1682810SN/Athe default in Eigen for performance reasons: contiguous storage allows 1692810SN/Avectorization that cannot be done when storage is not known to be contiguous at 1702810SN/Acompile time. The default ``Eigen::Ref`` stride type allows non-contiguous 1712810SN/Astorage along the outer dimension (that is, the rows of a column-major matrix 1724626SN/Aor columns of a row-major matrix), but not along the inner dimension. 1732810SN/A 1742810SN/AThis type, however, has the added benefit of also being able to map numpy array 1752810SN/Aslices. For example, the following (contrived) example uses Eigen with a numpy 1762810SN/Aslice to multiply by 2 all coefficients that are both on even rows (0, 2, 4, 1772810SN/A...) and in columns 2, 5, or 8: 1784626SN/A 1792810SN/A.. code-block:: cpp 1804666SN/A 1812810SN/A m.def("scale", [](py::EigenDRef<Eigen::MatrixXd> m, double c) { m *= c; }); 1822810SN/A 1839347SAndreas.Sandberg@arm.com.. code-block:: python 1849347SAndreas.Sandberg@arm.com 1859347SAndreas.Sandberg@arm.com # a = np.array(...) 18610509SAli.Saidi@ARM.com scale_by_2(myarray[0::2, 2:9:3]) 1879347SAndreas.Sandberg@arm.com 1889347SAndreas.Sandberg@arm.comThe second approach to avoid copying is more intrusive: rearranging the 1899347SAndreas.Sandberg@arm.comunderlying data types to not run into the non-contiguous storage problem in the 1909347SAndreas.Sandberg@arm.comfirst place. In particular, that means using matrices with ``Eigen::RowMajor`` 1912810SN/Astorage, where appropriate, such as: 1922810SN/A 1932810SN/A.. code-block:: cpp 1942810SN/A 1952810SN/A using RowMatrixXd = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>; 1962810SN/A // Use RowMatrixXd instead of MatrixXd 1972810SN/A 1982810SN/ANow bound functions accepting ``Eigen::Ref<RowMatrixXd>`` arguments will be 1994666SN/Acallable with numpy's (default) arrays without involving a copying. 2004666SN/A 2012810SN/AYou can, alternatively, change the storage order that numpy arrays use by 2022810SN/Aadding the ``order='F'`` option when creating an array: 2032810SN/A 2042810SN/A.. code-block:: python 20510679Sandreas.hansson@arm.com 2062810SN/A myarray = np.array(source, order='F') 20710679Sandreas.hansson@arm.com 2082810SN/ASuch an object will be passable to a bound function accepting an 2094908SN/A``Eigen::Ref<MatrixXd>`` (or similar column-major Eigen type). 2104908SN/A 2114908SN/AOne major caveat with this approach, however, is that it is not entirely as 2122810SN/Aeasy as simply flipping all Eigen or numpy usage from one to the other: some 2132810SN/Aoperations may alter the storage order of a numpy array. For example, ``a2 = 2142810SN/Aarray.transpose()`` results in ``a2`` being a view of ``array`` that references 2152810SN/Athe same data, but in the opposite storage order! 2164626SN/A 2172810SN/AWhile this approach allows fully optimized vectorized calculations in Eigen, it 2184666SN/Acannot be used with array slices, unlike the first approach. 2192810SN/A 2204626SN/AWhen *returning* a matrix to Python (either a regular matrix, a reference via 2212810SN/A``Eigen::Ref<>``, or a map/block into a matrix), no special storage 2222810SN/Aconsideration is required: the created numpy array will have the required 2232810SN/Astride that allows numpy to properly interpret the array, whatever its storage 2242810SN/Aorder. 2254666SN/A 2262810SN/AFailing rather than copying 2272810SN/A=========================== 22810192Smitch.hayenga@arm.com 22910192Smitch.hayenga@arm.comThe default behaviour when binding ``Eigen::Ref<const MatrixType>`` eigen 23010192Smitch.hayenga@arm.comreferences is to copy matrix values when passed a numpy array that does not 23110192Smitch.hayenga@arm.comconform to the element type of ``MatrixType`` or does not have a compatible 23210192Smitch.hayenga@arm.comstride layout. If you want to explicitly avoid copying in such a case, you 23310192Smitch.hayenga@arm.comshould bind arguments using the ``py::arg().noconvert()`` annotation (as 23410192Smitch.hayenga@arm.comdescribed in the :ref:`nonconverting_arguments` documentation). 23510192Smitch.hayenga@arm.com 23610192Smitch.hayenga@arm.comThe following example shows an example of arguments that don't allow data 23710192Smitch.hayenga@arm.comcopying to take place: 23810192Smitch.hayenga@arm.com 23910192Smitch.hayenga@arm.com.. code-block:: cpp 24010192Smitch.hayenga@arm.com 24110192Smitch.hayenga@arm.com // The method and function to be bound: 24210192Smitch.hayenga@arm.com class MyClass { 24310192Smitch.hayenga@arm.com // ... 2442810SN/A double some_method(const Eigen::Ref<const MatrixXd> &matrix) { /* ... */ } 2452813SN/A }; 2462810SN/A float some_function(const Eigen::Ref<const MatrixXf> &big, 24710766Sandreas.hansson@arm.com const Eigen::Ref<const MatrixXf> &small) { 2482810SN/A // ... 2492813SN/A } 2502810SN/A 2512810SN/A // The associated binding code: 2525715Shsul@eecs.umich.edu using namespace pybind11::literals; // for "arg"_a 2532810SN/A py::class_<MyClass>(m, "MyClass") 2542810SN/A // ... other class definitions 2559725Sandreas.hansson@arm.com .def("some_method", &MyClass::some_method, py::arg().noconvert()); 2562810SN/A 2572810SN/A m.def("some_function", &some_function, 2582810SN/A "big"_a.noconvert(), // <- Don't allow copying for this arg 2592810SN/A "small"_a // <- This one can be copied if needed 2602810SN/A ); 2612810SN/A 2622810SN/AWith the above binding code, attempting to call the the ``some_method(m)`` 2632810SN/Amethod on a ``MyClass`` object, or attempting to call ``some_function(m, m2)`` 2642810SN/Awill raise a ``RuntimeError`` rather than making a temporary copy of the array. 2652810SN/AIt will, however, allow the ``m2`` argument to be copied into a temporary if 2662810SN/Anecessary. 2679347SAndreas.Sandberg@arm.com 2689347SAndreas.Sandberg@arm.comNote that explicitly specifying ``.noconvert()`` is not required for *mutable* 2699347SAndreas.Sandberg@arm.comEigen references (e.g. ``Eigen::Ref<MatrixXd>`` without ``const`` on the 2709347SAndreas.Sandberg@arm.com``MatrixXd``): mutable references will never be called with a temporary copy. 2719347SAndreas.Sandberg@arm.com 2729347SAndreas.Sandberg@arm.comVectors versus column/row matrices 2739347SAndreas.Sandberg@arm.com================================== 2749347SAndreas.Sandberg@arm.com 2759347SAndreas.Sandberg@arm.comEigen and numpy have fundamentally different notions of a vector. In Eigen, a 2769347SAndreas.Sandberg@arm.comvector is simply a matrix with the number of columns or rows set to 1 at 2779347SAndreas.Sandberg@arm.comcompile time (for a column vector or row vector, respectively). Numpy, in 2789347SAndreas.Sandberg@arm.comcontast, has comparable 2-dimensional 1xN and Nx1 arrays, but *also* has 2799347SAndreas.Sandberg@arm.com1-dimensional arrays of size N. 280 281When passing a 2-dimensional 1xN or Nx1 array to Eigen, the Eigen type must 282have matching dimensions: That is, you cannot pass a 2-dimensional Nx1 numpy 283array to an Eigen value expecting a row vector, or a 1xN numpy array as a 284column vector argument. 285 286On the other hand, pybind11 allows you to pass 1-dimensional arrays of length N 287as Eigen parameters. If the Eigen type can hold a column vector of length N it 288will be passed as such a column vector. If not, but the Eigen type constraints 289will accept a row vector, it will be passed as a row vector. (The column 290vector takes precendence when both are supported, for example, when passing a 2911D numpy array to a MatrixXd argument). Note that the type need not be 292expicitly a vector: it is permitted to pass a 1D numpy array of size 5 to an 293Eigen ``Matrix<double, Dynamic, 5>``: you would end up with a 1x5 Eigen matrix. 294Passing the same to an ``Eigen::MatrixXd`` would result in a 5x1 Eigen matrix. 295 296When returning an eigen vector to numpy, the conversion is ambiguous: a row 297vector of length 4 could be returned as either a 1D array of length 4, or as a 2982D array of size 1x4. When encoutering such a situation, pybind11 compromises 299by considering the returned Eigen type: if it is a compile-time vector--that 300is, the type has either the number of rows or columns set to 1 at compile 301time--pybind11 converts to a 1D numpy array when returning the value. For 302instances that are a vector only at run-time (e.g. ``MatrixXd``, 303``Matrix<float, Dynamic, 4>``), pybind11 returns the vector as a 2D array to 304numpy. If this isn't want you want, you can use ``array.reshape(...)`` to get 305a view of the same data in the desired dimensions. 306 307.. seealso:: 308 309 The file :file:`tests/test_eigen.cpp` contains a complete example that 310 shows how to pass Eigen sparse and dense data types in more detail. 311