diff --git a/src/XrnetUtils.cpp b/src/XrnetUtils.cpp
index d9633f4..d07fe34 100644
--- a/src/XrnetUtils.cpp
+++ b/src/XrnetUtils.cpp
@@ -64,3 +64,11 @@ Eigen::MatrixXd computeResponseRcpp(SEXP X,
 double logit_inv(double x) {
     return 1 / (1 + std::exp(-x));
 }
+
+Eigen::Map<const Eigen::MatrixXd> createEigenMapFromRcppNumericMatrix(const Rcpp::NumericMatrix& mat) {
+    if (mat.rows() == 0 || mat.cols() == 0) {
+        Eigen::MatrixXd emptyMat;
+        return Eigen::Map<const Eigen::MatrixXd>(emptyMat.data(), emptyMat.rows(), emptyMat.cols());
+    }
+    return Eigen::Map<const Eigen::MatrixXd>(mat.begin(), mat.rows(), mat.cols());
+}
diff --git a/src/XrnetUtils.h b/src/XrnetUtils.h
index e25b48c..ba53b59 100644
--- a/src/XrnetUtils.h
+++ b/src/XrnetUtils.h
@@ -16,6 +16,8 @@ void compute_penalty(Eigen::Ref<Eigen::VectorXd> path,
 
 double logit_inv(double x);
 
+Eigen::Map<const Eigen::MatrixXd> createEigenMapFromRcppNumericMatrix(const Rcpp::NumericMatrix& mat);
+
 template <typename TX>
 Eigen::MatrixXd computeResponse(const TX & X,
                                 const Eigen::Ref<const Eigen::MatrixXd> & Fixed,
diff --git a/src/rcpp_wrappers_fit.cpp b/src/rcpp_wrappers_fit.cpp
index d82cec3..6acb681 100644
--- a/src/rcpp_wrappers_fit.cpp
+++ b/src/rcpp_wrappers_fit.cpp
@@ -6,6 +6,7 @@
 #include "GaussianSolver.h"
 #include "BinomialSolver.h"
 
+
 template <typename TX, typename TZ>
 Rcpp::List fitModel(const TX & x,
                     const bool & is_sparse_x,
@@ -198,7 +199,7 @@ Rcpp::List fitModelRcpp(SEXP x,
                 );
         else {
             Rcpp::NumericMatrix ext_mat(ext);
-            MapMat extmap((const double *) &ext_mat[0], ext_mat.rows(), ext_mat.cols());
+            MapMat extmap = createEigenMapFromRcppNumericMatrix(ext);
             return fitModel<MapMat, MapMat>(
                     xmap, is_sparse_x, y, extmap, fixed, weights_user,
                     intr, stnd, penalty_type, cmult, quantiles, num_penalty,
@@ -221,7 +222,7 @@ Rcpp::List fitModelRcpp(SEXP x,
         }
         else {
             Rcpp::NumericMatrix ext_mat(ext);
-            MapMat extmap((const double *) &ext_mat[0], ext_mat.rows(), ext_mat.cols());
+            MapMat extmap = createEigenMapFromRcppNumericMatrix(ext);
             return fitModel<MapMat, MapMat>(
                     xmap, is_sparse_x, y, extmap, fixed, weights_user, intr, stnd,
                     penalty_type, cmult, quantiles, num_penalty,
@@ -241,7 +242,7 @@ Rcpp::List fitModelRcpp(SEXP x,
             );
         else {
             Rcpp::NumericMatrix ext_mat(ext);
-            MapMat extmap((const double *) &ext_mat[0], ext_mat.rows(), ext_mat.cols());
+            MapMat extmap = createEigenMapFromRcppNumericMatrix(ext);
             return fitModel<MapSpMat, MapMat>(
                     Rcpp::as<MapSpMat>(x), is_sparse_x, y, extmap, fixed, weights_user,
                     intr, stnd, penalty_type, cmult, quantiles, num_penalty,