<p>Polynomial systems arising from the practice are often highly sparse, that is, the number of isolated solutions of a polynomial system is generally far less than its Bézout number. Therefore, the full exploration of the sparsity is an important topic in the field of homotopy method for solving polynomial systems. In this paper, the authors exploit the product structure of each polynomial to characterize the sparsity and further present a numerical method based on polynomial decomposition, in which the homotopy is the combination of the random product homotopy and the coefficient-parameter homotopy and the method is the combination of the symbolic methods and the numerical methods, to solve polynomial systems. Numerical results show that the proposed polynomial decomposition algorithm is more efficient than the existing homotopy methods in some cases, especially when the system has both sparse and dense polynomials.</p>