  We present a novel specialization of the variational Monte Carlo linear
method for the optimization of the recently introduced cluster Jastrow
antisymmetric geminal power ansatz, achieving a lower-order polynomial cost
scaling than would be possible with a naive application of the linear method
and greatly improving optimization performance relative to the previously
employed quasi-Newton approach. We test the methodology on highly
multi-reference triple-bond stretches, achieving accuracies superior to
traditional coupled cluster theory and multi-reference perturbation theory in
both the typical example of N$_2$ and the transition-metal-oxide example of
[ScO]$^+$.
