We present an interpolation method for efficient approximation of parametrized functions. The method recognizes and exploits the low-dimensional manifold structure of the parametrized functions to provide good approximation. Basic ingredients include a specific problem-dependent basis set defining a low-dimensional representation of the parametrized functions, and a set of ?best interpolation points? capturing the spatial-parameter variation of the parametrized functions. The best interpolation points are defined as solution of a least-squares minimization problem which can be solved efficiently using standard optimization algorithms. The approximation is then determined from the basis set and the best interpolation points through an inexpensive and stable interpolation procedure. In addition, an a posteriori error estimator is introduced to quantify the approximation error and requires little additional cost. Numerical results are presented to demonstrate the accuracy and efficiency of the method.