Let (a) over right arrow := (a(1),..., a(n)) is an element of [1, infinity)(n), p is an element of (0, 1), and alpha := 1/p - 1. For any x is an element of R-n and t is an element of [0, infinity), let Phi(p)(x, t) := {t/1 + (t vertical bar x vertical bar(v)(a))(1-p) if nu alpha is not an element of N, t/1 + (t[x](a)(v))(1-p) [log (e + vertical bar x vertical bar a())]p if nu alpha is not an element of N, where [center dot]((a) over right arrow) := 1+vertical bar center dot vertical bar((a) over right arrow), vertical bar center dot vertical bar((a) over right arrow) denotes the anisotropic quasi-homogeneous norm with respect to (a) over right arrow, and v := a(1)+ ... + a(n). Let H-d(p) (R-n), L-alpha(a)(R-n), and H-d(Phi p)(R-n) be, respectively, the anisotropic Hardy space, the anisotropic Campanato space, and the anisotropic Musielak-Orlicz Hardy space associated with Phi(p) on R-n. In this article, via first establishing the wavelet characterization of anisotropic Campanato spaces, we prove that for any f is an element of H-a(p)(R-n) and g is an element of L-alpha(a)(R-n), the product of f and g can be decomposed into S(f, g) + T(f, g) in the sense of tempered distributions, where S is a bilinear operator bounded from H-(a) over right arrow(p)(R-n) x L-alpha((a) over right arrow) (R-n) to L-1(R-n) and T is a bilinear operator bounded from H-(a) over right arrow(p) (R-n) x L-alpha((a) over right arrow) (R-n) to H-(a) over right arrow(Phi p) (R-n). Moreover, this bilinear decomposition is sharp in the dual sense that any Y. H-(d) over right arrow(Phi p) (R-n) that fits into the above bilinear decomposition should satisfy (L-1(R-n) + Y)* = (L-1(R-n) + H-(d) over right arrow(Phi p) (R-n))*. As applications, for any non-constant b is an element of L-alpha((a) over right arrow)(R-n) and any sublinear operator T satisfying some mild bounded assumptions, we find the largest subspace of H-(a) over right arrow(p) (R-n), denoted by H-(a) over right arrow ,b(p) (R-n), such that the commutator [b, T] is bounded from H-(a) over right arrow ,b(p)(R-n) to L-1(R-n). In addition, when T is an anisotropic Calderon-Zygmund operator, the boundedness of [b, T] from H-(a) over right arrow ,b(p)(R-n) to L-1(R-n) (or to H-(a) over right arrow(1) (R-n)) is also presented. The key of their proofs is the wavelet characterization of function spaces under consideration.