A fourth order polynomial yield criterion (Poly4) and a non-quadratic yield function are coupled under associated flow rule to describe the evolution of anisotropic yielding behavior analytically. The Poly4 part takes a great role in describing the evolution of these anisotropic behaviors throughout deformation history. The non-quadratic part is isotropic under uniaxial tension loading, which is supposed to control the curvature of yield surface. The parameters in the proposed yield criterion can be determined by analytical expressions using the four hardening curves along 0, 45, 90° and biaxial directions and also r-values, which doesn't need an optimization process and an interpolation method for describing the evolution of yield loci. Beyond the r-values along 0, 45, and 90°, additional two r-values can be used to analytically determine the nine parameters based on the availability of the data. To guarantee the convexity of the proposed yield criterion, an analytical method is developed by using the method of order principal minor determinant. The accuracy and effectiveness of the proposed yield criterion are verified with six different materials and compared with the Yld2000-2d yield criterion. The results show that the proposed yield function maintains great accuracy for anisotropic hardening with increasing equivalent plastic strain (EPS) beyond the initial yield loci. This yield criterion not only can accurately predict the yield stresses along four monotonic loading conditions, but also significantly reduces the error along other monotonic loading conditions. The proposed yield criterion is implemented into the implicit FEM code AutoForm R8 using user yield function in the development version and it accurately predicted the eight ears of AA3104-H19 and AA5042-H2 after a circular cup drawing.