docs: def get_Q(self):
          """Get the model's estimate of Q = \mu P \mu^T
  
          We can then separately extract \mu subject to additional constraints,
          e.g. \mu P 1 = diag(O).
          """
          Z = self.Z.detach().clone().numpy()
          O = self.O.numpy()
          I_k = np.eye(self.k)
          return O @ Z @ np.linalg.inv(I_k + Z.T @ O @ Z) @ Z.T @ O