Update HOGP docstring (#2709)

Summary:
Came across this error in the context of #2695, added a bit more explanation as well.

Pull Request resolved: #2709

Reviewed By: saitcakmak

Differential Revision: D68807233

Pulled By: Balandat

fbshipit-source-id: f3c3bd6373592465171d2f938faa8dcc30e32a3f

botorch/models/higher_order_gp.py

@@ -163,11 +163,14 @@ class HigherOrderGP(BatchedMultiOutputGPyTorchModel, ExactGP, FantasizeMixin):
     NOTE: This model requires the use of specialized Kronecker solves in
     linear operator, which are disabled by default in BoTorch. These are enabled
     by default in the `HigherOrderGP.posterior` call. However, they need to be
-    manually enabled by the user during model fitting.
+    manually enabled by the user during model fitting. Note also that we're using
+    `fit_gpytorch_mll_torch()` here instead of `fit_gpytorch_mll()` since the
+    approximate computations result in a non-smooth MLL that the default
+    L-BFGS-B optimizer invoked by `fit_gpytorch_mll()` does not handle well.
 
     Example:
         >>> from linear_operator.settings import _fast_solves
-        >>> model = SingleTaskGP(train_X, train_Y)
+        >>> model = HigherOrderGP(train_X, train_Y)
         >>> mll = ExactMarginalLogLikelihood(model.likelihood, model)
         >>> with _fast_solves(True):
         >>>     fit_gpytorch_mll_torch(mll)