code: def _construct_X_M(self, omega, **kwargs):
          
          X = self._construct_X(omega, weighted=True, **kwargs)
          M = np.dot(X.T, X)
  
          if getattr(self, 'regularization', None) is not None:
              diag = M.ravel(order='K')[::M.shape[0] + 1]
              if self.regularize_by_trace:
                  diag += diag.sum() * np.asarray(self.regularization)
              else:
                  diag += np.asarray(self.regularization)
  
          return X, M
docs: Construct the weighted normal matrix of the problem