diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h
index c6bda1115..95842c160 100644
--- a/Eigen/src/QR/SelfAdjointEigenSolver.h
+++ b/Eigen/src/QR/SelfAdjointEigenSolver.h
@@ -105,6 +105,25 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     /** \returns the computed eigen values */
     RealVectorType eigenvalues(void) const { return m_eivalues; }
 
+    /** \returns the positive square root of the matrix
+      *
+      * \note the matrix itself must be positive in order for this to make sense.
+      */
+    MatrixType operatorSqrt() const
+    {
+      return m_eivec * m_eivalues.cwise().sqrt().asDiagonal() * m_eivec.adjoint();
+    }
+
+    /** \returns the positive inverse square root of the matrix
+      *
+      * \note the matrix itself must be positive definite in order for this to make sense.
+      */
+    MatrixType operatorInverseSqrt() const
+    {
+      return m_eivec * m_eivalues.cwise().inverse().cwise().sqrt().asDiagonal() * m_eivec.adjoint();
+    }
+
+
   protected:
     MatrixType m_eivec;
     RealVectorType m_eivalues;
diff --git a/test/eigensolver.cpp b/test/eigensolver.cpp
index c093e9474..653edb7ea 100644
--- a/test/eigensolver.cpp
+++ b/test/eigensolver.cpp
@@ -108,6 +108,9 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox(
           symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal().eval()), largerEps));
 
+  MatrixType sqrtSymmA = eiSymm.operatorSqrt();
+  VERIFY(symmA.isApprox(sqrtSymmA*sqrtSymmA, ei_sqrt(test_precision<RealScalar>())));
+  VERIFY(sqrtSymmA.isApprox(symmA*eiSymm.operatorInverseSqrt(), ei_sqrt(test_precision<RealScalar>())));
 }
 
 template<typename MatrixType> void eigensolver(const MatrixType& m)