The linear system in Quasi Newton method

Question

I have implemented a Quasi Newton method for my problem, where I use the Hessian matrix approximation based approach. Hence, there is a linear system solve in every iteration. I solve the linear system using Conjugate Gradients and use the compact representation given in Representations of quasi-Newton matrices and their use in limited memory methods to get matrix vector products. I am noticing that the linear system is solved in very less number of iterations, 20 to 30 usually. This happens for even very large problems (one million variables). Is this observed in general when working with Quasi Newton? I can think of one possible reason for this fast convergence. Since, I ensure the curvature condition on the Quasi Newton vector pair, the Hessian approximation is guaranteed to be strictly positive definite and the Quadratic form in Conjugate Gradients has a "nice" shape. I am not sure whether this sufficiently explains what is happening.

Edit: Based on the accepted answer below, I have this question. The situation in the Quasi Newton approach, where the linear system can be provably solved in very less number of iterations is very interesting. Are there many works where this has been favorably exploited?

Answer 1

It is a well-known theorem that conjugate gradients is guaranteed to converge (in exact arithmetic) within $r+1$ steps for an identity-plus-rank-$r$ matrix, regardless of its size. In finite precision, you would find that in the worst case, CG would nearly stagnate for the first $r+1$ iterations, then rapidly converge towards the solution for every iteration after that.

The limited memory Hessian approximation is exactly an identity-plus-rank-$r$ matrix, where $r$ is (possibly a small scalar multiple of) the number of vectors you are currently carrying. Hence CG is guaranteed to solve it in $O(r)$ iterations. Assuming that the matrix-vector products with the Hessian are efficiently performed, this is a total time complexity of $O(nr^2)$.

Actually, CG is not the most efficient way to implement quasi-Newton. Instead of forming the limited memory approximation of the Hessian and then inverting it afterwards, you can form the limited memory approximation of the Hessian-inverse directly, using steps no more expensive than the ones you have been performing. Then the matrix solve at each iteration is replaced by a single matrix-vector multiplication, which can be performed in $O(nr)$ complexity.

Proofs for all of these statements can be found in any standard textbook on nonlinear optimization, such as Nocedal & Wright.