Accurate and efficient groupwise registration is important for population analysis. Current groupwise registration methods suffer from high computational cost, which hinders their application to large image datasets. To alleviate the computational burden while delivering accurate groupwise registration result, we propose to use a hierarchical graph set to model the complex image distribution with possibly large anatomical variations, and then turn the groupwise registration problem as a series of simple-to-solve graph shrinkage problems. Specifically, first, we divide the input images into a set of image clusters hierarchically, where images within each image cluster have similar anatomical appearances whereas images falling into different image clusters have varying anatomical appearances. After clustering, two types of graphs, i.e., intra-graph and inter-graph, are employed to hierarchically model the image distribution both within and across the image clusters. The constructed hierarchical graph set divides the registration problem of the whole image set into a series of simple-to-solve registration problems, where the entire registration process can be solved accurately and efficiently. The final deformation pathway of each image to the estimated population center can be obtained by composing each part of the deformation pathway along the hierarchical graph set. To evaluate our proposed method, we performed registration of a hundred of brain images with large anatomical variations. The results indicate that our method yields significant improvement in registration performance over state-of-the-art groupwise registration methods.