Graph filters are a fundamental tool in the field of graph signal processing. This paper focuses on the design of autoregressive moving average (ARMA) graph filters. In the proposed algorithm, the denominator part of an ARMA graph filter is decomposed as a cascade of a few second-order factors (SOFs) and a higher order factor (HOF), whose coefficients are updated sequentially in each iteration. In the proposed design algorithm, stability constraints are only imposed on the roots of SOFs and coefficients of the HOF are left unconstrained to enhance the design accuracy. Moreover, the number of SOFs can be automatically determined, which is convenient in practical applications. Simulation results demonstrate that the proposed algorithm can achieve higher computational efficiency and also approximation accuracy, compared to the state-of-the-arts of graph filters.