code: def optimize_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 = get_weight(x, y, perc).astype(bool)
          if issparse(weights): weights = weights.A
      else:
          weights = None
  
      x, y = x.astype(np.float64), y.astype(np.float64)
  
      n_vars = x.shape[1]
      offset, gamma = np.zeros(n_vars), np.zeros(n_vars)
  
      for i in range(n_vars):
          xi = x[:, i] if weights is None else x[:, i][weights[:, i]]
          yi = y[:, i] if weights is None else y[:, i][weights[:, i]]
  
          if fit_offset:
              offset[i], gamma[i] = minimize(lambda m: np.sum((-yi + xi * m[1] + m[0])**2), method="L-BFGS-B",
                                             x0=(0, 0.1), bounds=[(0, None), (None, None)]).x
          else:
              gamma[i] = minimize(lambda m: np.sum((-yi + xi * m) ** 2), x0=0.1, method="L-BFGS-B").x
      offset[np.isnan(offset)], gamma[np.isnan(gamma)] = 0, 0
      return offset, gamma
docs: Just to compare with closed-form solution