diff --git a/include/eigenpy/angle-axis.hpp b/include/eigenpy/angle-axis.hpp
index 7957e480f9c12501f8a16004ac8dde03f8deb06d..ea78f4dff615ad9b0ab3fad268642d80630b2f73 100644
--- a/include/eigenpy/angle-axis.hpp
+++ b/include/eigenpy/angle-axis.hpp
@@ -61,13 +61,15 @@ namespace eigenpy
       .def(bp::init<Matrix3>
            ((bp::arg("self"),bp::arg("rotation matrix")),
             "Initialize from a rotation matrix"))
-      .def(bp::init<Quaternion>((bp::arg("self"),bp::arg("quaternion")),"Initialize from a quaternion."))
-      .def(bp::init<AngleAxis>((bp::arg("self"),bp::arg("copy")),"Copy constructor."))
+      .def(bp::init<Quaternion>((bp::arg("self"),bp::arg("quaternion")),
+                                "Initialize from a quaternion."))
+      .def(bp::init<AngleAxis>((bp::arg("self"),bp::arg("copy")),
+                               "Copy constructor."))
       
       /* --- Properties --- */
       .add_property("axis",
-                    bp::make_function((const Vector3 & (AngleAxis::*)()const)&AngleAxis::axis,
-                                      bp::return_value_policy<bp::copy_const_reference>()),
+                    bp::make_function((Vector3 & (AngleAxis::*)())&AngleAxis::axis,
+                                      bp::return_internal_reference<>()),
                     &AngleAxisVisitor::setAxis,"The rotation axis.")
       .add_property("angle",
                     (Scalar (AngleAxis::*)()const)&AngleAxis::angle,
@@ -138,7 +140,7 @@ namespace eigenpy
                             bp::no_init)
       .def(AngleAxisVisitor<AngleAxis>());
       
-      // Cast to Eigen::RotationBase and vice-versa
+      // Cast to Eigen::RotationBase
       bp::implicitly_convertible<AngleAxis,RotationBase>();
     }
 
diff --git a/include/eigenpy/decompositions/EigenSolver.hpp b/include/eigenpy/decompositions/EigenSolver.hpp
index 36eb02949345d68c69e5e407501ca3e7df84506e..3c53acfaec7a3093e08e46761ed0c01bce13c3ab 100644
--- a/include/eigenpy/decompositions/EigenSolver.hpp
+++ b/include/eigenpy/decompositions/EigenSolver.hpp
@@ -6,6 +6,7 @@
 #define __eigenpy_decomposition_eigen_solver_hpp__
 
 #include "eigenpy/eigenpy.hpp"
+#include "eigenpy/eigen-to-python.hpp"
 
 #include <Eigen/Core>
 #include <Eigen/Eigenvalues>
@@ -37,31 +38,30 @@ namespace eigenpy
        
       .def("eigenvalues",&Solver::eigenvalues,bp::arg("self"),
            "Returns the eigenvalues of given matrix.",
-           bp::return_value_policy<bp::return_by_value>())
+           bp::return_internal_reference<>())
       .def("eigenvectors",&Solver::eigenvectors,bp::arg("self"),
-           "Returns the eigenvectors of given matrix.",
-           bp::return_value_policy<bp::return_by_value>())
+           "Returns the eigenvectors of given matrix.")
       
       .def("compute",&EigenSolverVisitor::compute_proxy<MatrixType>,
            bp::args("self","matrix"),
            "Computes the eigendecomposition of given matrix.",
-           bp::return_value_policy<bp::reference_existing_object>())
+           bp::return_self<>())
       .def("compute",(Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> & matrix, bool))&Solver::compute,
            bp::args("self","matrix","compute_eigen_vectors"),
            "Computes the eigendecomposition of given matrix.",
-           bp::return_value_policy<bp::reference_existing_object>())
+           bp::return_self<>())
       
       .def("getMaxIterations",&Solver::getMaxIterations,bp::arg("self"),
            "Returns the maximum number of iterations.")
       .def("setMaxIterations",&Solver::setMaxIterations,bp::args("self","max_iter"),
            "Sets the maximum number of iterations allowed.",
-           bp::return_value_policy<bp::reference_existing_object>())
+           bp::return_self<>())
       
       .def("pseudoEigenvalueMatrix",&Solver::pseudoEigenvalueMatrix,bp::arg("self"),
            "Returns the block-diagonal matrix in the pseudo-eigendecomposition.")
       .def("pseudoEigenvectors",&Solver::pseudoEigenvectors ,bp::arg("self"),
            "Returns the pseudo-eigenvectors of given matrix.",
-           bp::return_value_policy<bp::return_by_value>())
+           bp::return_internal_reference<>())
       
       .def("info",&Solver::info,bp::arg("self"),
            "NumericalIssue if the input contains INF or NaN values or overflow occured. Returns Success otherwise.")
diff --git a/include/eigenpy/decompositions/LDLT.hpp b/include/eigenpy/decompositions/LDLT.hpp
index 6e6b3e3d9b6812fb8a317beb91268f18ccc86066..6c30ae645e34009b0f001566138a631ad844c53e 100644
--- a/include/eigenpy/decompositions/LDLT.hpp
+++ b/include/eigenpy/decompositions/LDLT.hpp
@@ -53,22 +53,22 @@ namespace eigenpy
       
       .def("matrixLDLT",&Solver::matrixLDLT,bp::arg("self"),
            "Returns the LDLT decomposition matrix.",
-           bp::return_value_policy<bp::return_by_value>())
+           bp::return_internal_reference<>())
       
       .def("rankUpdate",(Solver & (Solver::*)(const Eigen::MatrixBase<VectorType> &, const RealScalar &))&Solver::template rankUpdate<VectorType>,
            bp::args("self","vector","sigma"),
-           bp::return_value_policy<bp::reference_existing_object>() )
+           bp::return_self<>())
     
 #if EIGEN_VERSION_AT_LEAST(3,3,0)
       .def("adjoint",&Solver::adjoint,bp::arg("self"),
            "Returns the adjoint, that is, a reference to the decomposition itself as if the underlying matrix is self-adjoint.",
-           bp::return_value_policy<bp::reference_existing_object>())
+           bp::return_self<>())
 #endif
       
       .def("compute",(Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> & matrix))&Solver::compute,
            bp::args("self","matrix"),
            "Computes the LDLT of given matrix.",
-           bp::return_value_policy<bp::reference_existing_object>())
+           bp::return_self<>())
       
       .def("info",&Solver::info,bp::arg("self"),
            "NumericalIssue if the input contains INF or NaN values or overflow occured. Returns Success otherwise.")
diff --git a/include/eigenpy/decompositions/LLT.hpp b/include/eigenpy/decompositions/LLT.hpp
index f854aba368401fc869bb6845adde8c4d2f743ea3..2dfce7ce603f54e5396407bf7e55a2f402e61ddf 100644
--- a/include/eigenpy/decompositions/LLT.hpp
+++ b/include/eigenpy/decompositions/LLT.hpp
@@ -43,7 +43,7 @@ namespace eigenpy
            "Returns the upper triangular matrix U.")
       .def("matrixLLT",&Solver::matrixLLT,bp::arg("self"),
            "Returns the LLT decomposition matrix.",
-           bp::return_value_policy<bp::return_by_value>())
+           bp::return_internal_reference<>())
       
       .def("rankUpdate",(Solver (Solver::*)(const VectorType &, const RealScalar &))&Solver::template rankUpdate<VectorType>,
            bp::args("self","vector","sigma"))
@@ -51,13 +51,13 @@ namespace eigenpy
 #if EIGEN_VERSION_AT_LEAST(3,3,0)
       .def("adjoint",&Solver::adjoint,bp::arg("self"),
            "Returns the adjoint, that is, a reference to the decomposition itself as if the underlying matrix is self-adjoint.",
-           bp::return_value_policy<bp::reference_existing_object>())
+           bp::return_self<>())
 #endif
       
       .def("compute",(Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> & matrix))&Solver::compute,
            bp::args("self","matrix"),
            "Computes the LLT of given matrix.",
-           bp::return_value_policy<bp::reference_existing_object>())
+           bp::return_self<>())
       
       .def("info",&Solver::info,bp::arg("self"),
            "NumericalIssue if the input contains INF or NaN values or overflow occured. Returns Success otherwise.")
diff --git a/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp b/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
index 562ba0c5bca53c72f0b1eb79a8493b024ec01061..ee129c7dbda3b7efbbf9ff1701f33e7c3f3438aa 100644
--- a/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
+++ b/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
@@ -6,6 +6,7 @@
 #define __eigenpy_decomposition_self_adjoint_eigen_solver_hpp__
 
 #include "eigenpy/eigenpy.hpp"
+#include "eigenpy/eigen-to-python.hpp"
 
 #include <Eigen/Core>
 #include <Eigen/Eigenvalues>
@@ -37,10 +38,10 @@ namespace eigenpy
        
       .def("eigenvalues",&Solver::eigenvalues,bp::arg("self"),
            "Returns the eigenvalues of given matrix.",
-           bp::return_value_policy<bp::return_by_value>())
+           bp::return_internal_reference<>())
       .def("eigenvectors",&Solver::eigenvectors,bp::arg("self"),
            "Returns the eigenvectors of given matrix.",
-           bp::return_value_policy<bp::return_by_value>())
+           bp::return_internal_reference<>())
       
       .def("compute",&SelfAdjointEigenSolverVisitor::compute_proxy<MatrixType>,
            bp::args("self","matrix"),
@@ -49,23 +50,21 @@ namespace eigenpy
       .def("compute",(Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> & matrix, int options))&Solver::compute,
            bp::args("self","matrix","options"),
            "Computes the eigendecomposition of given matrix.",
-           bp::return_value_policy<bp::reference_existing_object>())
+           bp::return_self<>())
       
       .def("computeDirect",&SelfAdjointEigenSolverVisitor::computeDirect_proxy,
            bp::args("self","matrix"),
            "Computes eigendecomposition of given matrix using a closed-form algorithm.",
-           bp::return_value_policy<bp::reference_existing_object>())
+           bp::return_self<>())
       .def("computeDirect",(Solver & (Solver::*)(const MatrixType & matrix, int options))&Solver::computeDirect,
            bp::args("self","matrix","options"),
            "Computes eigendecomposition of given matrix using a closed-form algorithm.",
-           bp::return_value_policy<bp::reference_existing_object>())
+           bp::return_self<>())
        
       .def("operatorInverseSqrt",&Solver::operatorInverseSqrt,bp::arg("self"),
-           "Computes the inverse square root of the matrix.",
-           bp::return_value_policy<bp::return_by_value>())
+           "Computes the inverse square root of the matrix.")
       .def("operatorSqrt",&Solver::operatorSqrt,bp::arg("self"),
-           "Computes the inverse square root of the matrix.",
-           bp::return_value_policy<bp::return_by_value>())
+           "Computes the inverse square root of the matrix.")
       
       .def("info",&Solver::info,bp::arg("self"),
            "NumericalIssue if the input contains INF or NaN values or overflow occured. Returns Success otherwise.")
diff --git a/include/eigenpy/quaternion.hpp b/include/eigenpy/quaternion.hpp
index d1462f2e8c1c1c06ecd29ae01b57a19198899e08..6d5f1ffc8c83bf32e4d248a7ac95937f3346dbdf 100644
--- a/include/eigenpy/quaternion.hpp
+++ b/include/eigenpy/quaternion.hpp
@@ -155,7 +155,7 @@ namespace eigenpy
       .def("coeffs",(const Vector4 & (Quaternion::*)()const)&Quaternion::coeffs,
            bp::arg("self"),
            "Returns a vector of the coefficients (x,y,z,w)",
-           bp::return_value_policy<bp::copy_const_reference>())
+           bp::return_internal_reference<>())
       .def("matrix",&Quaternion::matrix,
            bp::arg("self"),
            "Returns an equivalent 3x3 rotation matrix. Similar to toRotationMatrix.")
@@ -216,9 +216,13 @@ namespace eigenpy
       .def("__setitem__",&QuaternionVisitor::__setitem__)
       .def("__getitem__",&QuaternionVisitor::__getitem__)
       .def("assign",&assign<Quaternion>,
-           bp::args("self","quat"),"Set *this from an quaternion quat and returns a reference to *this.",bp::return_self<>())
+           bp::args("self","quat"),
+           "Set *this from an quaternion quat and returns a reference to *this.",
+           bp::return_self<>())
       .def("assign",(Quaternion & (Quaternion::*)(const AngleAxis &))&Quaternion::operator=,
-           bp::args("self","aa"),"Set *this from an angle-axis aa and returns a reference to *this.",bp::return_self<>())
+           bp::args("self","aa"),
+           "Set *this from an angle-axis aa and returns a reference to *this.",
+           bp::return_self<>())
       .def("__str__",&print)
       .def("__repr__",&print)
       
@@ -251,9 +255,9 @@ namespace eigenpy
     { return self = quat; }
 
     static Quaternion* FromTwoVectors(const Vector3& u, const Vector3& v)
-    { 
+    {
       Quaternion* q(new Quaternion); q->setFromTwoVectors(u,v);
-      return q; 
+      return q;
     }
   
     static bool __eq__(const Quaternion & u, const Quaternion & v)