This article introduces Hessian approximation algorithms to estimate the search direction of the quasi-Newton methods for solving optimization problems of continuous parameters. The proposed algorithms are quite different from other well-known quasi-Newton methods, such as symmetric rank-one, Davidon-Fletcher-Powell, and Broyden-Fletcher-Goldfarb-Shanno, in that the Hessian matrix is not calculated from the gradient information, rather directly from the function values. The proposed algorithms are designed for a class of hybrid algorithms that combine evolutionary search with the gradient-based methods of quasi-Newton type. The function values calculated for the evolutionary search are used for estimation of the Hessian matrix (or its inverse) as well as the gradient vector. Since the estimation process of the Hessian matrix is independent of that of the gradient vector, more reliable Hessian estimation with a small population is possible compared with the previous methods based upon the classical quasi-Newton methods. Numerical experiments show that the proposed algorithms are very competitive with state-of-the-art evolutionary algorithms for continuous optimization problems.