diff --git a/src/gwmodelpp/GWRMultiscale.cpp b/src/gwmodelpp/GWRMultiscale.cpp
index ed16fac..b19d0e4 100644
--- a/src/gwmodelpp/GWRMultiscale.cpp
+++ b/src/gwmodelpp/GWRMultiscale.cpp
@@ -426,17 +426,17 @@ vec GWRMultiscale::fitVarSerial(const vec &x, const vec &y, const uword var, mat
     if (mHasHatMatrix)
     {
         mat ci, si;
-        S = mat(mHasHatMatrix ? nDp : 1, nDp, fill::zeros);
+        S = mat(nDp, nDp, fill::zeros);
         for (uword i = 0; i < nDp  ; i++)
         {
             GWM_LOG_STOP_BREAK(mStatus);
             vec w = mSpatialWeights[var].weightVector(i);
             mat xtw = trans(x % w);
-            mat xtwx = xtw * x;
-            mat xtwy = xtw * y;
+            double xtwx = as_scalar(xtw * x);
+            double xtwy = as_scalar(xtw * y);
             try
             {
-                mat xtwx_inv = inv_sympd(xtwx);
+                double xtwx_inv = 1.0 / xtwx;
                 betas.col(i) = xtwx_inv * xtwy;
                 ci = xtwx_inv * xtw;
                 si = x(i) * ci;
@@ -570,12 +570,12 @@ double GWRMultiscale::bandwidthSizeCriterionVarCVSerial(BandwidthWeight *bandwid
         vec w = bandwidthWeight->weight(d);
         w(i) = 0.0;
         mat xtw = trans(mXi % w);
-        mat xtwx = xtw * mXi;
-        mat xtwy = xtw * mYi;
+        double xtwx = as_scalar(xtw * mXi);
+        double xtwy = as_scalar(xtw * mYi);
         try
         {
-            mat xtwx_inv = inv_sympd(xtwx);
-            vec beta = xtwx_inv * xtwy;
+            double xtwx_inv = 1.0 / xtwx;
+            double beta = xtwx_inv * xtwy;
             double res = mYi(i) - det(mXi(i) * beta);
             cv += res * res;
         }
@@ -607,11 +607,11 @@ double GWRMultiscale::bandwidthSizeCriterionVarAICSerial(BandwidthWeight *bandwi
         vec d = mSpatialWeights[var].distance()->distance(i);
         vec w = bandwidthWeight->weight(d);
         mat xtw = trans(mXi % w);
-        mat xtwx = xtw * mXi;
-        mat xtwy = xtw * mYi;
+        double xtwx = as_scalar(xtw * mXi);
+        double xtwy = as_scalar(xtw * mYi);
         try
         {
-            mat xtwx_inv = inv_sympd(xtwx);
+            double xtwx_inv = 1.0 / xtwx;
             betas.col(i) = xtwx_inv * xtwy;
             mat ci = xtwx_inv * xtw;
             mat si = mXi(i) * ci;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f5fe3bc..f808657 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -29,6 +29,7 @@ set(SAMPLE_DATA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/data)
 add_definitions(-DSAMPLE_DATA_DIR="${SAMPLE_DATA_DIR}")
 
 set(ADDON_SOURCES
+    data/londonhp.cpp
     data/londonhp100.cpp
     TerminateCheckTelegram.cpp
     FileTelegram.cpp
diff --git a/test/testGWRMultiscale.cpp b/test/testGWRMultiscale.cpp
index 5d70a99..5ceda76 100644
--- a/test/testGWRMultiscale.cpp
+++ b/test/testGWRMultiscale.cpp
@@ -10,6 +10,7 @@
 #include "gwmodelpp/spatialweight/CRSDistance.h"
 #include "gwmodelpp/spatialweight/BandwidthWeight.h"
 #include "gwmodelpp/spatialweight/SpatialWeight.h"
+#include "londonhp.h"
 #include "londonhp100.h"
 #include "TerminateCheckTelegram.h"
 #include "FileTelegram.h"
@@ -336,3 +337,47 @@ TEST_CASE("Multiscale GWR: cancel")
     }
 
 }
+
+TEST_CASE("Large Data Repeate", "[benchmark]")
+{
+    mat londonhp_coord, londonhp_data;
+    vector<string> londonhp_fields;
+    if (!read_londonhp(londonhp_coord, londonhp_data, londonhp_fields))
+    {
+        FAIL("Cannot load londonhp100 data.");
+    }
+
+    vec y = londonhp_data.col(0);
+    mat x = join_rows(ones(londonhp_data.n_rows), londonhp_data.cols(uvec({1, 3})));
+    uword nVar = 3;
+
+    BENCHMARK("specified bandwidth")
+    {
+        vector<SpatialWeight> spatials;
+        vector<bool> preditorCentered;
+        vector<GWRMultiscale::BandwidthInitilizeType> bandwidthInitialize;
+        vector<GWRMultiscale::BandwidthSelectionCriterionType> bandwidthSelectionApproach;
+        for (size_t i = 0; i < nVar; i++)
+        {
+            CRSDistance distance;
+            BandwidthWeight bandwidth(64, true, BandwidthWeight::Bisquare);
+            spatials.push_back(SpatialWeight(&bandwidth, &distance));
+            preditorCentered.push_back(i != 0);
+            bandwidthInitialize.push_back(GWRMultiscale::BandwidthInitilizeType::Initial);
+            bandwidthSelectionApproach.push_back(GWRMultiscale::BandwidthSelectionCriterionType::CV);
+        }
+
+        GWRMultiscale algorithm;
+        algorithm.setCoords(londonhp_coord);
+        algorithm.setDependentVariable(y);
+        algorithm.setIndependentVariables(x);
+        algorithm.setSpatialWeights(spatials);
+        algorithm.setHasHatMatrix(true);
+        algorithm.setPreditorCentered(preditorCentered);
+        algorithm.setBandwidthInitilize(bandwidthInitialize);
+        algorithm.setBandwidthSelectionApproach(bandwidthSelectionApproach);
+        algorithm.setBandwidthSelectThreshold(vector(3, 1e-5));
+        algorithm.fit();
+        return algorithm.betas();
+    };
+}