docs: @Override
    public <M, E, T extends E> SufficientStatistics train(GradientOracle<M, E> oracle,
        SufficientStatistics initialParameters, Iterable<T> trainingData) {
      int circularBufferSize = numVectorsInApproximation;
      SufficientStatistics currentParameters = initialParameters;
      SufficientStatistics nextParameters = currentParameters.duplicate();
      SufficientStatistics previousParameters = currentParameters.duplicate();
      SufficientStatistics previousGradient = currentParameters.duplicate();
      SufficientStatistics direction = currentParameters.duplicate();
      List<SufficientStatistics> pointDeltas = Lists.newArrayList(Collections.nCopies(circularBufferSize, null));
      List<SufficientStatistics> gradientDeltas = Lists.newArrayList(Collections.nCopies(circularBufferSize, null));
      List<Double> scalings = Lists.newArrayList(Collections.nCopies(circularBufferSize, null));
  
      MapReduceExecutor executor = MapReduceConfiguration.getMapReduceExecutor();
      List<T> dataList = Lists.newArrayList(trainingData);
      GradientEvaluation gradientEvaluation = null;
      for (int i = 0; i < maxIterations || maxIterations < 0; i++) {
        log.notifyIterationStart(i);
        log.logParameters(i, initialParameters);
  
        if (gradientEvaluation == null) {
          gradientEvaluation = evaluateGradient(currentParameters, dataList,
              oracle, executor, log);
        }
        SufficientStatistics gradient = gradientEvaluation.getGradient();
        
        double gradientL2Norm = gradient.getL2Norm();
        if (gradientL2Norm < gradientConvergenceThreshold) {
          return currentParameters;
        }
  
        // We haven't converged yet. Figure out which direction to move in.
        log.startTimer("lbfgs/compute_search_direction");
        // Store the requisite data for approximating the inverse
        // Hessian.
        if (i > 0) {
          int index = (i - 1) % circularBufferSize;
          SufficientStatistics pointDelta = pointDeltas.get(index);
          
          if (pointDelta == null) {
            pointDelta = currentParameters.duplicate();
          } else {
            pointDelta.zeroOut();
            pointDelta.increment(currentParameters, 1.0);
          }
  
          pointDelta.increment(previousParameters, -1.0);
          pointDeltas.set(index, pointDelta);
  
          // Note that gradient and previousGradient are actually the
          // *negative* gradient (i.e., a descent direction).
          SufficientStatistics gradientDelta = gradientDeltas.get(index);
          
          if (gradientDelta == null) {
            gradientDelta = previousGradient.duplicate();
          } else {
            gradientDelta.zeroOut();
            gradientDelta.increment(previousGradient, 1.0);
          }
          gradientDelta.increment(gradient, -1.0);
          gradientDeltas.set(index, gradientDelta);
  
          scalings.set(index, 1.0 / (pointDelta.innerProduct(gradientDelta)));
        }
  
        previousParameters.zeroOut();
        previousParameters.increment(currentParameters, 1.0);
        previousGradient.zeroOut();
        previousGradient.increment(gradient, 1.0);
  
        // Compute this iteration's search direction.
        // TODO: deal with the memory here. 
        int hessianVectorCount = (int) Math.min(numVectorsInApproximation, i);
        direction.zeroOut();
        direction.increment(gradient, 1.0);
        direction.multiply(-1.0);
        double[] weights = new double[hessianVectorCount];
        for (int j = 0; j < hessianVectorCount; j++) {
          int index = (i - (j + 1)) % circularBufferSize;
          double weight = scalings.get(index) * (pointDeltas.get(index).innerProduct(direction));
          direction.increment(gradientDeltas.get(index), -1.0 * weight);
          weights[hessianVectorCount - (j + 1)] = weight;
        }
  
        // The assumption here is that the initial Hessian estimate is
        // the identity. Multiply direction by the Hessian estimate here
        // to pick another value.
        for (int j = 0; j < hessianVectorCount; j++) {
          int index = (i + j - hessianVectorCount) % circularBufferSize;
          double weight = scalings.get(index) * (gradientDeltas.get(index).innerProduct(direction));
          direction.increment(pointDeltas.get(index), weights[j] - weight);
        }
        log.stopTimer("lbfgs/compute_search_direction");
  
        log.logStatistic(i, "parameter l2 norm", previousParameters.getL2Norm());
        log.logStatistic(i, "gradient l2 norm", gradientL2Norm);
        log.logStatistic(i, "direction l2 norm", direction.getL2Norm());
        log.logStatistic(i, "search errors", gradientEvaluation.getSearchErrors());
        log.logStatistic(i, "objective value", gradientEvaluation.getObjectiveValue());
  
        log.startTimer("lbfgs/compute_step_size");
        // Perform a backtracking line search to find a step size.
        double stepSize = 1.0 / LINE_SEARCH_CONSTANT;
        double currentObjectiveValue = gradientEvaluation.getObjectiveValue();
        double nextObjectiveValue, curInnerProd, cond1Rhs, nextInnerProd, cond2Rhs;
  
        do {
          stepSize = stepSize * LINE_SEARCH_CONSTANT;
          // if (nextParameters == null) {
          nextParameters = currentParameters.duplicate();
          // } else {
            // nextParameters.zeroOut();
            // nextParameters.increment(currentParameters, 1.0);
          // }
          nextParameters.increment(direction, -1.0 * stepSize);
          gradientEvaluation = evaluateGradient(nextParameters, dataList, oracle, executor, log);
  
          // Check the Wolfe conditions to ensure sufficient descent.
          nextObjectiveValue = gradientEvaluation.getObjectiveValue();
          curInnerProd = gradient.innerProduct(direction);
          SufficientStatistics nextGradient = gradientEvaluation.getGradient();
          nextInnerProd = nextGradient.innerProduct(direction);
  
          cond1Rhs = currentObjectiveValue - (WOLFE_CONDITION_C1 * stepSize * curInnerProd);
          cond2Rhs = -1.0 * WOLFE_CONDITION_C2 * curInnerProd;
        } while ((nextObjectiveValue <= cond1Rhs || Double.isNaN(cond1Rhs) || Double.isNaN(nextObjectiveValue))
            && stepSize > minStepSize);
  
        log.logStatistic(i, "step size", stepSize);
        log.stopTimer("lbfgs/compute_step_size");
  
        if (stepSize <= minStepSize) {
          throw new LbfgsConvergenceError("L-BFGS could not find a suitable step size.",
              currentParameters, direction, stepSize, i);
        }
        currentParameters = nextParameters;
  
        log.notifyIterationEnd(i);
      }
      return currentParameters;
    }