This work develops a computationally efficient sector-model subspace iteration to obtain accurate eigensolutions for the vibration of multi-stage systems where each stage is cyclically symmetric. The method fully captures the inter-stage coupling. The component stages are cyclically symmetric but the full multi-stage model is generally asymmetric because the number of sectors can differ among the stages. The sector-model subspace iteration is equivalent to the full-model subspace iteration in its accuracy, but it is more computationally efficient. Different from the full-model subspace iteration that relies on a full multi-stage system model, the sector-model subspace iteration only uses matrices of single-stage sector models. Parallel computation applies to most steps in each iteration and many matrices do not need to be recomputed for each iteration. The improved computational efficiency of the sector-model subspace iteration is more significant when the full model has a large number of stages, a large number of sectors within each stage, and/or many degrees of freedom in each single-stage sector. The iteration converges rapidly and can start with any initial trial vectors. Faster convergence is achieved using high-quality trial vectors from reduced-order modeling methods, such as single-stage analyses, multi-stage cyclic symmetry reduction, and component mode synthesis. In contrast to reduced-order modeling methods, the iteration provides eigensolutions with controllable accuracy based on user-selected convergence criteria. The application and advantages of the sector-model subspace iteration are demonstrated using two example multi-stage bladed disks.