diff --git a/.gitignore b/.gitignore
index 174f8247..2bd819eb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -30,3 +30,5 @@ cmake-build-debug/
 #docs files
 *.puml
 .Rproj.user
+
+experiments/
\ No newline at end of file
diff --git a/.gitmodules b/.gitmodules
index ce546b35..f96e5a93 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,4 +1,4 @@
 [submodule "fdaPDE/core"]
 	path = fdaPDE/core
-	url = https://github.com/fdaPDE/fdaPDE-core
+	url = https://github.com/marcogalliani/fdaPDE-core
 	branch = stable
diff --git a/fdaPDE/core b/fdaPDE/core
index a2a9c88b..f1a90ebf 160000
--- a/fdaPDE/core
+++ b/fdaPDE/core
@@ -1 +1 @@
-Subproject commit a2a9c88b203b6f1b74bbc71d698f131f839f25ba
+Subproject commit f1a90ebffc41c9eefa8d4df7291254ab256414f4
diff --git a/fdaPDE/src/models/fpca.h b/fdaPDE/src/models/fpca.h
index 2036524b..f602d1d2 100644
--- a/fdaPDE/src/models/fpca.h
+++ b/fdaPDE/src/models/fpca.h
@@ -77,7 +77,8 @@ template <typename VariationalSolver> class fpca_power_iteration_impl {
             case OptimizeGCV: {
                 auto gcv_functor = [&](auto lambda) { return gcv_(X, lambda, V.col(i)); };
                 GridSearch<n_lambda> optimizer;
-                opt_lambda = optimizer.optimize(gcv_functor, lambda_grid);
+                auto opt_ = optimizer.optimize(gcv_functor, lambda_grid);
+                for (int i = 0; i < n_lambda; ++i) { opt_lambda[i] = opt_[i]; }
             } break;
             case OptimizeMSRE: {
             } break;
@@ -134,12 +135,12 @@ template <typename VariationalSolver> class fpca_power_iteration_impl {
     double gcv_(const matrix_t& X, const LambdaT lambda, const InitT& f0) {
         const auto& [f, s] = solve_(X, lambda, f0);
         // evaluate GCV index at convergence
-        std::array<double, n_lambda> lambda_vec;
-        std::copy(lambda.data(), lambda.data() + n_lambda, lambda_vec.begin());
-        if (edf_map_.find(lambda_vec) == edf_map_.end()) {   // cache Tr[S]
-            edf_map_[lambda_vec] = smoother_->edf();
+	std::array<double, n_lambda> lambda_;
+	for (int i = 0; i < n_lambda; ++i) { lambda_[i] = lambda[i]; }
+        if (edf_map_.find(lambda_) == edf_map_.end()) {   // cache Tr[S]
+            edf_map_[lambda_] = smoother_->edf();
         }
-        int dor = n_locs_ - edf_map_.at(lambda_vec);
+        int dor = n_locs_ - edf_map_.at(lambda_);
         return (n_locs_ / std::pow(dor, 2)) * ((smoother_->Psi() * f) - smoother_->response()).squaredNorm();
     }
     std::unordered_map<std::array<double, n_lambda>, double, internals::std_array_hash<double, n_lambda>> edf_map_;
@@ -199,7 +200,8 @@ template <typename VariationalSolver> class fpca_subspace_iteration_impl {
         case OptimizeGCV: {
             auto gcv_functor = [&](auto lambda) { return gcv_(X, rank, lambda, V); };
             GridSearch<n_lambda> optimizer;
-            opt_lambda = optimizer.optimize(gcv_functor, lambda_grid);
+            auto opt_ = optimizer.optimize(gcv_functor, lambda_grid);
+            for (int i = 0; i < n_lambda; ++i) { opt_lambda[i] = opt_[i]; }
         } break;
         case OptimizeMSRE: {
         } break;
@@ -313,7 +315,8 @@ template <typename VariationalSolver> class fpca_direct_impl {
         case OptimizeGCV: {
             auto gcv_functor = [&](auto lambda) { return gcv_(X, rank, lambda, flag); };
             GridSearch<n_lambda> optimizer;
-            opt_lambda = optimizer.optimize(gcv_functor, lambda_grid);
+            auto opt_ = optimizer.optimize(gcv_functor, lambda_grid);
+            for (int i = 0; i < n_lambda; ++i) { opt_lambda[i] = opt_[i]; }
         } break;
         case OptimizeMSRE: {
         } break;
@@ -371,7 +374,7 @@ template <typename VariationalSolver> class fpca_direct_impl {
     double gcv_(const matrix_t& X, int rank, const LambdaT lambda, int flag) {
         const auto& [F, S] = solve_(X, rank, lambda, flag);
         // evaluate GCV index at convergence (note that Tr[S] = \|D^(-1)\|_F^2)
-        int dor = n_locs_ - invD_.squaredNorm();
+        int dor = n_locs_ -  (invD_*smoother_->Psi().transpose()).squaredNorm();
         return (n_locs_ / std::pow(dor, 2)) * (X.transpose() * S - (smoother_->Psi() * F)).squaredNorm();
     }
     int n_locs_ = 0, n_units_ = 0, n_dofs_ = 0;