Hello,
i think i have figured it out:
I have set up the optimizer with history_size = 3 and max_iter = 1. After each optimizer.step() call you can print the optimizer state with print(optimizer.state[optimizer._params[0]]) and the length of the old directories which are taken into account in each iteration with print(len(optimizer.state[optimizer._params[0]]['old_dirs'])). The output of the latter is:
0
1
2
3
3
3
…
So the two loop recursion
              # iteration in L-BFGS loop collapsed to use just one buffer
                q = flat_grad.neg()
                for i in range(num_old - 1, -1, -1):
                    al[i] = old_stps[i].dot(q) * ro[i]
                    q.add_(old_dirs[i], alpha=-al[i])

                # multiply by initial Hessian
                # r/d is the final direction
                d = r = torch.mul(q, H_diag)
                for i in range(num_old):
                    be_i = old_dirs[i].dot(r) * ro[i]
                    r.add_(old_stps[i], alpha=al[i] - be_i)

is correct (see e.g. [Jorge Nocedal, Stephen Wright: Numerical Optimization, Algorithm 7.4]). This means that max_iter is really just the maximal number of iterations per optimization step and therefore max_iter >1 doesn’t have an effect. Moreover max_eval is also not really reasonable. Btw. after  [Jorge Nocedal, Stephen Wright: Numerical Optimization] it is often sufficient to set history_size between 5-20 and the line search is necessary such that the BFGS update is stable.