code: def leastsq_NxN(x, y, fit_offset=False, perc=None):
      
      if perc is not None:
          if not fit_offset and isinstance(perc, (list, tuple)): perc = perc[1]
          weights = csr_matrix(get_weight(x, y, perc)).astype(bool)
          x, y = weights.multiply(x).tocsr(), weights.multiply(y).tocsr()
      else:
          weights = None
  
      with warnings.catch_warnings():
          warnings.simplefilter("ignore")
          xx_ = prod_sum_obs(x, x)
          xy_ = prod_sum_obs(x, y)
  
          if fit_offset:
              n_obs = x.shape[0] if weights is None else sum_obs(weights)
              x_ = sum_obs(x) / n_obs
              y_ = sum_obs(y) / n_obs
              gamma = (xy_ / n_obs - x_ * y_) / (xx_ / n_obs - x_ ** 2)
              offset = y_ - gamma * x_
  
              # fix negative offsets:
              idx = offset < 0
              gamma[idx] = xy_[idx] / xx_[idx]
              offset = np.clip(offset, 0, None)
          else:
              gamma = xy_ / xx_
              offset = np.zeros(x.shape[1])
      offset[np.isnan(offset)], gamma[np.isnan(gamma)] = 0, 0
      return offset, gamma
docs: Solution to least squares: gamma = cov(X,Y) / var(X)