docs: def _orthogonalize(X):
      """ Orthogonalize every column of design `X` w.r.t preceding columns
  
      Parameters
      ----------
      X: array of shape(n, p)
         the data to be orthogonalized
  
      Returns
      -------
      X: array of shape(n, p)
         the data after orthogonalization
  
      Notes
      -----
      X is changed in place. The columns are not normalized
      """
      if X.size == X.shape[0]:
          return X
      from scipy.linalg import pinv, norm
      for i in range(1, X.shape[1]):
          X[:, i] -= np.dot(np.dot(X[:, i], X[:, :i]), pinv(X[:, :i]))
          # X[:, i] /= norm(X[:, i])
      return X