The classical Metropolis sampling method is a cornerstone of many statistical modeling applications that range from physics, chemistry, and biology to economics. This method is particularly suitable for sampling the thermal distributions of classical systems. The challenge of extending this method to the simulation of arbitrary quantum systems is that, in general, eigenstates of quantum Hamiltonians cannot be obtained efficiently with a classical computer. However, this challenge can be overcome by quantum computers.Here, we present a quantum algorithm which fully generalizes the classical Metropolis algorithm to the quantum domain. The meaning of quantum generalization is twofold: The proposed algorithm is not only applicable to both classical and quantum systems, but also offers a quantum speedup relative to the classical counterpart. Furthermore, unlike the classical method of quantum Monte Carlo, this quantum algorithm does not suffer from the negative-sign problem associated with fermionic systems. Applications of this algorithm include the study of low-temperature properties of quantum systems, such as the Hubbard model, and preparing the thermal states of sizable molecules to simulate, for example, chemical reactions at an arbitrary temperature.Metropolis method | statistical physics | quantum computing | quantum simulation | simulated annealing I nteracting many-body problems cover a wide range of applications in many areas in science and technology. These are best exemplified by the Ising spin model, which was originally developed for explaining ferromagnetism and its associated phase transitions (1). The Ising model was later found to be related to many practical applications, for example, error-correcting codes, image restoration, associative memory, and optimization problems (see, e.g., ref. 2). However, there is no efficient method for finding solutions to many-body problems in general. For example, the problem of finding the ground state of the Ising model with arbitrary couplings is known to be a nondeterministic polynomial (NP)-complete problem, meaning that, if an efficient algorithm for solving the Ising model exists, then all of the problems in the class of NP, for example the graph isomorphism problem and some variants of the traveling salesman problem, can be solved efficiently as well (see, e.g., ref.3).The most fundamental challenge for solving many-body problems is that they generally require an exponentially large amount of spatial and temporal computing resources to find the solutions as the system size increases (4). Although no universal method for solving general many-body problems has been found, powerful classical methods such as Markov-chain Monte Carlo (MCMC) or quantum Monte Carlo (QMC) have been invented and proven to be successful in many applications (see, e.g., ref. 5). These methods, however, have certain limitations. For example, the running time of MCMC scales as Oð1∕δÞ (6), where δ is the gap of the transition matrix. For problems such as spin glasses (7) where δ is vani...