Let $\mathbb F$ be a field and $P \in \mathbb F [x_1,\ldots, x_n]$ be a homogeneous polynomial such that $|\mathbb F| > °(P)$ and $φ, ψ\colon \mathbb F^n \to \mathbb F^n$ be two maps such that $P(\mathbf{x} + λ\mathbf{y}) = P(φ(\mathbf{x}) + λψ(\mathbf{y}))$ for all $λ\in \mathbb F$ and $\mathbf{x}, \mathbf{y} \in \mathbb F^n.$ We provide the characterization of all such $φ$ and $ψ$ for all polynomials in the case if $\mathrm{char}(\mathbb F) = 0$ and for all polynomials satisfying certain condition in the case if $\mathrm{char}(\mathbb F) > 0$. This characterization generalizes the existing results regarding the linear maps on matrices preserving the determinant, the immanant and other homogeneous polynomial functions of matrix entries. To obtain the main result of this paper, we introduce the vector space $\mathcal L_{P} \subseteq {\mathbb F^n}^*$ spanned by the range of the gradient field of $P \in \mathbb F[x_1,\ldots, x_n]$. Being a linear invariant associated with $P,$ this space has several remarkable properties and may also be used for studying the linear maps preserving $P$. In addition, we demonstrate how the main result could be applied to the particular polynomial matrix invariants. Namely, we provide an explicit description of corresponding pairs of nonlinear maps $φ, ψ$ for the case where $P$ is equal to the Cullis' determinant of $n\times k$ rectangular matrix (with the assumption that $n \ge k + 2$ and $k \ge 3$).