code: public static void twoLoopHp(Vec x_grad, List<Double> rho, List<Vec> s, List<Vec> y, Vec q, double[] alphas)
      {
          //q ← ∇ fk;
          x_grad.copyTo(q);
          if(s.isEmpty())
              return;//identity, we are done
          //for i = k−1,k−2,...,k−m
          for(int i = 0; i < s.size(); i++)
          {
              Vec s_i = s.get(i);
              Vec y_i = y.get(i);
              double alpha_i = alphas[i] = rho.get(i)*s_i.dot(q);
              q.mutableSubtract(alpha_i, y_i);
          }
          
          //r ← Hk0q;  and see eq (7.20), done in place in q
          q.mutableMultiply(s.get(0).dot(y.get(0))/y.get(0).dot(y.get(0)));
          //for i = k−m,k−m+1,...,k−1
          for(int i = s.size()-1; i >= 0; i--)
          {
              //β ← ρ_i y_i^T r ;
              double beta = rho.get(i)*y.get(i).dot(q);
              //r ← r + si (αi − β)
              q.mutableAdd(alphas[i]-beta, s.get(i));
          }
      }
docs: See Algorithm 7.4 (L-BFGS two-loop recursion).
  @param x_grad the initial value &nabla; f<sub>k</sub>
  @param rho
  @param s
  @param y
  @param q the location to store the value of H<sub>k</sub> &nabla; f<sub>k</sub>
  @param alphas temp space to do work, should be as large as the number of history vectors