LinAlgError with Sel_BW

[imadwalid]

Hello guys, I'am using MGWR and Sel_BW it throws an error but when when using a sample size > 2040 but works fine for < 2040

b_y = np.log(gdf["price"].head(2050).values.reshape((-1, 1)))
b_X = gdf["surface"].head(2050).values.reshape((-1, 1))
u = gdf['x'].head(2050)
v = gdf['y'].head(2050)
b_coords = list(zip(u, v))
b_X = (b_X - b_X.mean(axis = 0)) / b_X.std(axis = 0)
b_y = (b_y - b_y.mean(axis = 0)) / b_y.std(axis = 0)

mgwr_selector = Sel_BW(b_coords, b_y, b_X, multi=True)
print(mgwr_selector.fixed)
print(mgwr_selector.kernel)

mgwr_bw = mgwr_selector.search()
print(mgwr_bw)
mgwr_results = MGWR(b_coords, b_y, b_X, mgwr_selector).fit()

/usr/local/lib/python3.7/dist-packages/mgwr/kernels.py:60: RuntimeWarning: divide by zero encountered in true_divide
  self.kernel = self._kernel_funcs(self.dvec / self.bandwidth)
/usr/local/lib/python3.7/dist-packages/mgwr/kernels.py:60: RuntimeWarning: invalid value encountered in true_divide
  self.kernel = self._kernel_funcs(self.dvec / self.bandwidth)
---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-208-6365a372acef> in <module>()
----> 1 mgwr_bw = mgwr_selector.search(verbose=2)
      2 print(mgwr_bw)
      3 mgwr_results = MGWR(b_coords, b_y, b_X, mgwr_selector).fit()

12 frames
/usr/local/lib/python3.7/dist-packages/mgwr/sel_bw.py in search(self, search_method, criterion, bw_min, bw_max, interval, tol, max_iter, init_multi, tol_multi, rss_score, max_iter_multi, multi_bw_min, multi_bw_max, bws_same_times, pool, verbose)
    311 
    312         if self.multi:
--> 313             self._mbw()
    314             self.params = self.bw[3]  #params n by k
    315             self.sel_hist = self.bw[-2] #bw searching history

/usr/local/lib/python3.7/dist-packages/mgwr/sel_bw.py in _mbw(self)
    400                            self.max_iter_multi, self.rss_score, gwr_func,
    401                            bw_func, sel_func, multi_bw_min, multi_bw_max,
--> 402                            bws_same_times, verbose=self.verbose)
    403 
    404     def _init_section(self, X_glob, X_loc, coords, constant):

/usr/local/lib/python3.7/dist-packages/mgwr/search.py in multi_bw(init, y, X, n, k, family, tol, max_iter, rss_score, gwr_func, bw_func, sel_func, multi_bw_min, multi_bw_max, bws_same_times, verbose)
    221                 bw = bws[j]
    222             else:
--> 223                 bw = sel_func(bw_class, multi_bw_min[j], multi_bw_max[j])
    224                 gwr_sel_hist.append(deepcopy(bw_class.sel_hist))
    225 

/usr/local/lib/python3.7/dist-packages/mgwr/sel_bw.py in sel_func(bw_func, bw_min, bw_max)
    395                 search_method=search_method, criterion=criterion,
    396                 bw_min=bw_min, bw_max=bw_max, interval=interval, tol=tol,
--> 397                 max_iter=max_iter, pool=self.pool, verbose=False)
    398 
    399         self.bw = multi_bw(self.init_multi, y, X, n, k, family, self.tol_multi,

/usr/local/lib/python3.7/dist-packages/mgwr/sel_bw.py in search(self, search_method, criterion, bw_min, bw_max, interval, tol, max_iter, init_multi, tol_multi, rss_score, max_iter_multi, multi_bw_min, multi_bw_max, bws_same_times, pool, verbose)
    317                 -1]  #scalar, optimal bw from initial gwr model
    318         else:
--> 319             self._bw()
    320             self.sel_hist = self.bw[-1]
    321 

/usr/local/lib/python3.7/dist-packages/mgwr/sel_bw.py in _bw(self)
    337             self.bw = golden_section(a, c, delta, gwr_func, self.tol,
    338                                      self.max_iter, self.int_score,
--> 339                                      self.verbose)
    340         elif self.search_method == 'interval':
    341             self.bw = equal_interval(self.bw_min, self.bw_max, self.interval,

/usr/local/lib/python3.7/dist-packages/mgwr/search.py in golden_section(a, c, delta, function, tol, max_iter, int_score, verbose)
     60             score_b = dict[b]
     61         else:
---> 62             score_b = function(b)
     63             dict[b] = score_b
     64             if verbose:

/usr/local/lib/python3.7/dist-packages/mgwr/sel_bw.py in <lambda>(bw)
    327             self.coords, self.y, self.X_loc, bw, family=self.family, kernel=
    328             self.kernel, fixed=self.fixed, constant=self.constant, offset=self.
--> 329             offset, spherical=self.spherical).fit(lite=True, pool=self.pool))
    330 
    331         self._optimized_function = gwr_func

/usr/local/lib/python3.7/dist-packages/mgwr/gwr.py in fit(self, ini_params, tol, max_iter, solve, lite, pool)
    333                 rslt = map(self._local_fit, range(m))  #sequential
    334 
--> 335             rslt_list = list(zip(*rslt))
    336             influ = np.array(rslt_list[0]).reshape(-1, 1)
    337             resid = np.array(rslt_list[1]).reshape(-1, 1)

/usr/local/lib/python3.7/dist-packages/mgwr/gwr.py in _local_fit(self, i)
    249 
    250         if isinstance(self.family, Gaussian):
--> 251             betas, inv_xtx_xt = _compute_betas_gwr(self.y, self.X, wi)
    252             predy = np.dot(self.X[i], betas)[0]
    253             resid = self.y[i] - predy

/usr/local/lib/python3.7/dist-packages/spglm/iwls.py in _compute_betas_gwr(y, x, wi)
     35     xT = (x * wi).T
     36     xtx = np.dot(xT, x)
---> 37     xtx_inv_xt = linalg.solve(xtx, xT)
     38     betas = np.dot(xtx_inv_xt, y)
     39     return betas, xtx_inv_xt

/usr/local/lib/python3.7/dist-packages/scipy/linalg/basic.py in solve(a, b, sym_pos, lower, overwrite_a, overwrite_b, debug, check_finite, assume_a, transposed)
    214                                                (a1, b1))
    215         lu, ipvt, info = getrf(a1, overwrite_a=overwrite_a)
--> 216         _solve_check(n, info)
    217         x, info = getrs(lu, ipvt, b1,
    218                         trans=trans, overwrite_b=overwrite_b)

/usr/local/lib/python3.7/dist-packages/scipy/linalg/basic.py in _solve_check(n, info, lamch, rcond)
     29                          '.'.format(-info))
     30     elif 0 < info:
---> 31         raise LinAlgError('Matrix is singular.')
     32 
     33     if lamch is None:

LinAlgError: Matrix is singular.

train_paris.csv i tried changing kernel but nothing worked, i'am using only one feature the surface to predict the price, i'am adding the file so u can reproduce the error, thanks for your help.

[imadwalid]

i just tried changing the search criterion to BIC and worked, but didn't not with the other having the same problem !

[Ziqi-Li]

Hi @imadwalid There may be local clusters of similar values in the first 2050 records. A fix would be to set the lower bound of the golden-section search a bit higher. For example, in your dataset, setting to 50 works (default is like 40ish).

mgwr_bw = mgwr_selector.search(verbose=True,multi_bw_min=[50])