We present a polynomial hybrid Monte Carlo (PHMC) algorithm for lattice QCD with odd numbers of flavors of O(a)-improved Wilson quark action. The algorithm makes use of the non-Hermitian Chebyshev polynomial to approximate the inverse square root of the fermion matrix required for an odd number of flavors. The systematic error from the polynomial approximation is removed by a noisy Metropolis test for which a new method is developed. Investigating the property of our PHMC algorithm in the Nf=2 QCD case, we find that it is as efficient as the conventional HMC algorithm for a moderately large lattice size (163×48) with intermediate quark masses (mPS/mV∼0.7–0.8). We test our odd-flavor algorithm through extensive simulations of two-flavor QCD treated as an Nf=1+1 system, and comparing the results with those of the established algorithms for Nf=2 QCD. These tests establish that our PHMC algorithm works on a moderately large lattice size with intermediate quark masses (163×48,mPS/mV∼0.7–0.8). Finally we experiment with the (2+1)-flavor QCD simulation on small lattices (43×8 and 83×16), and confirm the agreement of our results with those obtained with the R algorithm and extrapolated to a zero molecular dynamics step size.