Ising machines are expected to solve combinatorial optimization problems efficiently by representing them as Ising models or equivalent quadratic unconstrained binary optimization (QUBO) models . However, upper bound exists on the computable problem size due to the hardware limitations of Ising machines. This paper propose a new hybrid annealing method based on partial QUBO extraction, called subQUBO model extraction, with multiple solution instances. For a given QUBO model, the proposed method obtains N I quasi-optimal solutions (quasi-ground-state solutions) in some way using a classical computer. The solutions giving these quasi-optimal solutions are called solution instances. We extract a size-limited subQUBO model as follows based on a strong theoretical background: we randomly select N S (N S < N I ) solution instances among them and focus on a particular binary variable x i in the N S solution instances. If x i value is much varied over N S solution instances, it is included in the subQUBO model; otherwise, it is not. We find a (quasi-)ground-state solution of the extracted subQUBO model using an Ising machine and add it as a new solution instance. By repeating this process, we can finally obtain a (quasi-)ground-state solution of the original QUBO model. Experimental evaluations confirm that the proposed method can obtain better quasi-ground-state solution than existing methods for large-sized QUBO models.