A procedure for energy optimizing an antisymmetrized geminal power state, which contains a description of correlation effects in an unbiased fashion, is presented. This procedure overcomes difficulties and shortcomings of past optimization procedures for antisymmetrized geminal power states. It is shown that the computational cost scales as r(5), where r is the size of the spin orbital atomic basis set, which is superior to the scaling costs of multiexcited complete active space self-consistent field calculations, but unlike such calculations their is no preselection of configurations or excitation level. The variational parameters are the geminal coefficients with respect to a fixed atomic orbital basis, which allows one to classify these variables into chemically significant and non significant ones according to their magnitudes and thus reduce size of calculations.