  When using Wannier functions to study the electronic structure of
multi-parameter Hamiltonians $H^{(\boldsymbol k,\bf \lambda)}$ carrying a
dependence on crystal momentum $\boldsymbol k$ and an additional periodic
parameter $\bf\lambda$, one usually constructs several sets of Wannier
functions for a set of values of $\bf\lambda$. We present the concept of higher
dimensional Wannier functions (HDWFs), which provide a minimal and accurate
description of the electronic structure of multi-parameter Hamiltonians based
on a single set of HDWFs. The obstacle of non-orthogonality of Bloch functions
at different $\bf\lambda$ is overcome by introducing an auxiliary real space,
which is reciprocal to the parameter $\bf\lambda$. We derive a generalized
interpolation scheme and emphasize the essential conceptual and computational
simplifications in using the formalism, for instance, in the evaluation of
linear response coefficients. We further implement the necessary machinery to
construct HDWFs from ab initio within the full-potential linearized augmented
plane-wave method (FLAPW). We apply our implementation to accurately
interpolate the Hamiltonian of a one-dimensional magnetic chain of Mn atoms in
two important cases of $\bf\lambda$: (i) the spin-spiral vector $\boldsymbol
q$, and (ii) the direction of the ferromagnetic magnetization $\hat{\boldsymbol
m}$. Using the generalized interpolation of the energy, we extract the
corresponding values of magneto-crystalline anisotropy energy, Heisenberg
exchange constants, and spin stiffness, which compare very well with the values
obtained from direct first principles calculations. For toy models we
demonstrate that the method of HDWFs can also be used in applications such as
the virtual crystal approximation, ferroelectric polarization and spin torques.
