diff --git a/src/shogun/mathematics/linalg/LinalgBackendBase.h b/src/shogun/mathematics/linalg/LinalgBackendBase.h
index cf2dbec009b..b139fc06d03 100644
--- a/src/shogun/mathematics/linalg/LinalgBackendBase.h
+++ b/src/shogun/mathematics/linalg/LinalgBackendBase.h
@@ -389,8 +389,8 @@ namespace shogun
  */
 #define BACKEND_GENERIC_IN_PLACE_MATRIX_PROD(Type, Container)                  \
 	virtual void matrix_prod(                                                  \
-	    SGMatrix<Type>& a, Container<Type>& b, Container<Type>& result,        \
-	    bool transpose_A, bool transpose_B) const                              \
+	    const SGMatrix<Type>& a, const Container<Type>& b,                     \
+	    Container<Type>& result, bool transpose_A, bool transpose_B) const     \
 	{                                                                          \
 		SG_SNOTIMPLEMENTED;                                                    \
 	}
@@ -529,7 +529,7 @@ namespace shogun
  */
 #define BACKEND_GENERIC_IN_PLACE_SCALE(Type, Container)                        \
 	virtual void scale(                                                        \
-	    Container<Type>& a, Type alpha, Container<Type>& result) const         \
+	    const Container<Type>& a, Type alpha, Container<Type>& result) const   \
 	{                                                                          \
 		SG_SNOTIMPLEMENTED;                                                    \
 	}
diff --git a/src/shogun/mathematics/linalg/LinalgBackendEigen.h b/src/shogun/mathematics/linalg/LinalgBackendEigen.h
index 639ea84caa0..905b0890e38 100644
--- a/src/shogun/mathematics/linalg/LinalgBackendEigen.h
+++ b/src/shogun/mathematics/linalg/LinalgBackendEigen.h
@@ -199,8 +199,8 @@ namespace shogun
 /** Implementation of @see LinalgBackendBase::matrix_prod */
 #define BACKEND_GENERIC_IN_PLACE_MATRIX_PROD(Type, Container)                  \
 	virtual void matrix_prod(                                                  \
-	    SGMatrix<Type>& a, Container<Type>& b, Container<Type>& result,        \
-	    bool transpose_A, bool transpose_B) const;
+	    const SGMatrix<Type>& a, const Container<Type>& b,                     \
+	    Container<Type>& result, bool transpose_A, bool transpose_B) const;
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_MATRIX_PROD, SGVector)
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_MATRIX_PROD, SGMatrix)
 #undef BACKEND_GENERIC_IN_PLACE_MATRIX_PROD
@@ -267,7 +267,7 @@ namespace shogun
 /** Implementation of @see linalg::scale */
 #define BACKEND_GENERIC_IN_PLACE_SCALE(Type, Container)                        \
 	virtual void scale(                                                        \
-	    Container<Type>& a, Type alpha, Container<Type>& result) const;
+	    const Container<Type>& a, Type alpha, Container<Type>& result) const;
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_SCALE, SGVector)
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_SCALE, SGMatrix)
 #undef BACKEND_GENERIC_IN_PLACE_SCALE
@@ -531,13 +531,13 @@ namespace shogun
 		/** Eigen3 matrix * vector in-place product method */
 		template <typename T>
 		void matrix_prod_impl(
-		    SGMatrix<T>& a, SGVector<T>& b, SGVector<T>& result, bool transpose,
-		    bool transpose_B = false) const;
+		    const SGMatrix<T>& a, const SGVector<T>& b, SGVector<T>& result,
+		    bool transpose, bool transpose_B = false) const;
 
 		/** Eigen3 matrix in-place product method */
 		template <typename T>
 		void matrix_prod_impl(
-		    SGMatrix<T>& a, SGMatrix<T>& b, SGMatrix<T>& result,
+		    const SGMatrix<T>& a, const SGMatrix<T>& b, SGMatrix<T>& result,
 		    bool transpose_A, bool transpose_B) const;
 
 		/** Return the largest element in the vector with Eigen3 library */
@@ -594,11 +594,13 @@ namespace shogun
 
 		/** Eigen3 vector inplace scale method: result = alpha * A */
 		template <typename T>
-		void scale_impl(SGVector<T>& a, T alpha, SGVector<T>& result) const;
+		void
+		scale_impl(const SGVector<T>& a, T alpha, SGVector<T>& result) const;
 
 		/** Eigen3 matrix inplace scale method: result = alpha * A */
 		template <typename T>
-		void scale_impl(SGMatrix<T>& a, T alpha, SGMatrix<T>& result) const;
+		void
+		scale_impl(const SGMatrix<T>& a, T alpha, SGMatrix<T>& result) const;
 
 		/** Eigen3 set const method */
 		template <typename T, template <typename> class Container>
diff --git a/src/shogun/mathematics/linalg/LinalgNamespace.h b/src/shogun/mathematics/linalg/LinalgNamespace.h
index 75c76f69ab2..b93cf0f3d2a 100644
--- a/src/shogun/mathematics/linalg/LinalgNamespace.h
+++ b/src/shogun/mathematics/linalg/LinalgNamespace.h
@@ -818,7 +818,8 @@ namespace shogun
 		 */
 		template <typename T>
 		void element_prod(
-		    Block<SGMatrix<T>>& a, Block<SGMatrix<T>>& b, SGMatrix<T>& result)
+		    const Block<SGMatrix<T>>& a, const Block<SGMatrix<T>>& b,
+		    SGMatrix<T>& result)
 		{
 			REQUIRE(
 			    a.m_row_size == b.m_row_size && a.m_col_size == b.m_col_size,
@@ -850,7 +851,8 @@ namespace shogun
 		 * @return The result of the operation
 		 */
 		template <typename T>
-		SGMatrix<T> element_prod(Block<SGMatrix<T>>& a, Block<SGMatrix<T>>& b)
+		SGMatrix<T>
+		element_prod(const Block<SGMatrix<T>>& a, const Block<SGMatrix<T>>& b)
 		{
 			REQUIRE(
 			    a.m_row_size == b.m_row_size && a.m_col_size == b.m_col_size,
@@ -877,7 +879,8 @@ namespace shogun
 		 * @param result Result matrix
 		 */
 		template <typename T>
-		void element_prod(SGMatrix<T>& a, SGMatrix<T>& b, SGMatrix<T>& result)
+		void element_prod(
+		    const SGMatrix<T>& a, const SGMatrix<T>& b, SGMatrix<T>& result)
 		{
 			REQUIRE(
 			    a.num_rows == b.num_rows && a.num_cols == b.num_cols,
@@ -912,7 +915,7 @@ namespace shogun
 		 * @return The result of the operation
 		 */
 		template <typename T>
-		SGMatrix<T> element_prod(SGMatrix<T>& a, SGMatrix<T>& b)
+		SGMatrix<T> element_prod(const SGMatrix<T>& a, const SGMatrix<T>& b)
 		{
 			REQUIRE(
 			    a.num_rows == b.num_rows && a.num_cols == b.num_cols,
@@ -1015,7 +1018,7 @@ namespace shogun
 		 */
 		template <typename T>
 		void matrix_prod(
-		    SGMatrix<T>& A, SGVector<T>& b, SGVector<T>& result,
+		    const SGMatrix<T>& A, const SGVector<T>& b, SGVector<T>& result,
 		    bool transpose = false)
 		{
 			if (transpose)
@@ -1069,8 +1072,8 @@ namespace shogun
 		 * @return result Result vector
 		 */
 		template <typename T>
-		SGVector<T>
-		matrix_prod(SGMatrix<T>& A, SGVector<T>& b, bool transpose = false)
+		SGVector<T> matrix_prod(
+		    const SGMatrix<T>& A, const SGVector<T>& b, bool transpose = false)
 		{
 			SGVector<T> result;
 			if (transpose)
@@ -1113,7 +1116,7 @@ namespace shogun
 		 */
 		template <typename T>
 		void matrix_prod(
-		    SGMatrix<T>& A, SGMatrix<T>& B, SGMatrix<T>& result,
+		    const SGMatrix<T>& A, const SGMatrix<T>& B, SGMatrix<T>& result,
 		    bool transpose_A = false, bool transpose_B = false)
 		{
 			REQUIRE(
@@ -1214,8 +1217,8 @@ namespace shogun
 		 */
 		template <typename T>
 		SGMatrix<T> matrix_prod(
-		    SGMatrix<T>& A, SGMatrix<T>& B, bool transpose_A = false,
-		    bool transpose_B = false)
+		    const SGMatrix<T>& A, const SGMatrix<T>& B,
+		    bool transpose_A = false, bool transpose_B = false)
 		{
 			SGMatrix<T> result;
 
@@ -1278,8 +1281,8 @@ namespace shogun
 		 */
 		template <typename T>
 		void dgemv(
-		    T alpha, SGMatrix<T> a, bool transpose, SGVector<T> x, T beta,
-		    SGVector<T>& y)
+		    T alpha, const SGMatrix<T>& a, bool transpose, const SGVector<T>& x,
+		    T beta, SGVector<T>& y)
 		{
 			auto temp_vector = matrix_prod(a, x, transpose);
 			add(temp_vector, y, y, alpha, beta);
@@ -1302,8 +1305,8 @@ namespace shogun
 		 */
 		template <typename T>
 		void dgemm(
-		    T alpha, SGMatrix<T> a, SGMatrix<T> b, bool transpose_a,
-		    bool transpose_b, T beta, SGMatrix<T>& c)
+		    T alpha, const SGMatrix<T>& a, const SGMatrix<T>& b,
+		    bool transpose_a, bool transpose_b, T beta, SGMatrix<T>& c)
 		{
 			auto temp_matrix = matrix_prod(a, b, transpose_a, transpose_b);
 			add(temp_matrix, c, c, alpha, beta);
@@ -1408,7 +1411,7 @@ namespace shogun
 		 * @param result The vector of alpha * a
 		 */
 		template <typename T>
-		void scale(SGVector<T>& a, SGVector<T>& result, T alpha = 1)
+		void scale(const SGVector<T>& a, SGVector<T>& result, T alpha = 1)
 		{
 			REQUIRE(
 			    result.vlen == a.vlen,
@@ -1428,7 +1431,7 @@ namespace shogun
 		 * @param result The matrix of alpha * A
 		 */
 		template <typename T>
-		void scale(SGMatrix<T>& A, SGMatrix<T>& result, T alpha = 1)
+		void scale(const SGMatrix<T>& A, SGMatrix<T>& result, T alpha = 1)
 		{
 			REQUIRE(
 			    (A.num_rows == result.num_rows), "Number of rows of matrix A "
@@ -1452,7 +1455,7 @@ namespace shogun
 		 * @return Vector or matrix of alpha * A
 		 */
 		template <typename T, template <typename> class Container>
-		Container<T> scale(Container<T>& a, T alpha = 1)
+		Container<T> scale(const Container<T>& a, T alpha = 1)
 		{
 			auto result = a.clone();
 			scale(a, result, alpha);
diff --git a/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp b/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp
index 9d2f9a52941..36c7b5b0b64 100644
--- a/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp
+++ b/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp
@@ -129,8 +129,8 @@ DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_EXPONENT, SGMatrix)
 
 #define BACKEND_GENERIC_IN_PLACE_MATRIX_PROD(Type, Container)                  \
 	void LinalgBackendEigen::matrix_prod(                                      \
-	    SGMatrix<Type>& a, Container<Type>& b, Container<Type>& result,        \
-	    bool transpose_A, bool transpose_B) const                              \
+	    const SGMatrix<Type>& a, const Container<Type>& b,                     \
+	    Container<Type>& result, bool transpose_A, bool transpose_B) const     \
 	{                                                                          \
 		matrix_prod_impl(a, b, result, transpose_A, transpose_B);              \
 	}
@@ -140,7 +140,7 @@ DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_MATRIX_PROD, SGMatrix)
 
 #define BACKEND_GENERIC_IN_PLACE_SCALE(Type, Container)                        \
 	void LinalgBackendEigen::scale(                                            \
-	    Container<Type>& a, Type alpha, Container<Type>& result) const         \
+	    const Container<Type>& a, Type alpha, Container<Type>& result) const   \
 	{                                                                          \
 		scale_impl(a, alpha, result);                                          \
 	}
@@ -303,8 +303,8 @@ void LinalgBackendEigen::exponent_impl(
 
 template <typename T>
 void LinalgBackendEigen::matrix_prod_impl(
-    SGMatrix<T>& a, SGVector<T>& b, SGVector<T>& result, bool transpose,
-    bool transpose_B) const
+    const SGMatrix<T>& a, const SGVector<T>& b, SGVector<T>& result,
+    bool transpose, bool transpose_B) const
 {
 	typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
 	typename SGVector<T>::EigenVectorXtMap b_eig = b;
@@ -318,8 +318,8 @@ void LinalgBackendEigen::matrix_prod_impl(
 
 template <typename T>
 void LinalgBackendEigen::matrix_prod_impl(
-    SGMatrix<T>& a, SGMatrix<T>& b, SGMatrix<T>& result, bool transpose_A,
-    bool transpose_B) const
+    const SGMatrix<T>& a, const SGMatrix<T>& b, SGMatrix<T>& result,
+    bool transpose_A, bool transpose_B) const
 {
 	typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
 	typename SGMatrix<T>::EigenMatrixXtMap b_eig = b;
@@ -340,7 +340,7 @@ void LinalgBackendEigen::matrix_prod_impl(
 
 template <typename T>
 void LinalgBackendEigen::scale_impl(
-    SGVector<T>& a, T alpha, SGVector<T>& result) const
+    const SGVector<T>& a, T alpha, SGVector<T>& result) const
 {
 	typename SGVector<T>::EigenVectorXtMap a_eig = a;
 	typename SGVector<T>::EigenVectorXtMap result_eig = result;
@@ -350,7 +350,7 @@ void LinalgBackendEigen::scale_impl(
 
 template <typename T>
 void LinalgBackendEigen::scale_impl(
-    SGMatrix<T>& a, T alpha, SGMatrix<T>& result) const
+    const SGMatrix<T>& a, T alpha, SGMatrix<T>& result) const
 {
 	typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
 	typename SGMatrix<T>::EigenMatrixXtMap result_eig = result;