Question:
Extensively searching the space of all possible vectors of length $n$ to satisfy a (non-overdetermined) requirement is possible in principle. Hence, there is a way to calculate a complex FIR kernel of length $n$ whose magnitude response goes through $2n$ given points, when its phase response is allowed to be completely arbitrary. (Corollary, a real FIR kernel with $n$ given magnitude points and $n$ mirrored points can be calculated.) The extensive search method is hopelessly inefficient, though. Is there a practical way to calculate an FIR with such requirements?
Context:
It is common DSP practice to use the DFT to convert an array of real magnitude points $M = \{m_0, m_1, \dots, m_{n-2}, m_{n-1}\}$, $m_i \in \mathbb{R}$ along with an array of real phase points $P = \{p_0, p_1, \dots, p_{n-2}, p_{n-1}\}$, $p_i \in \mathbb{R}$ to an FIR kernel $K = \{k_0, k_1, \dots, k_{n-2}, k_{n-1}\}$. $K$'s elements are generally complex, but if $M$ is symmetric (i.e. $m_i = m_{n-i}$ for $0 < i < \frac{n}{2}$) and $P$ is conjugate symmetric (i.e. $p_i = - p_{n-i}$ for $0 < i < \frac{n}{2}$), then $P$'s elements are completely real.
The magnitude response (absolute value of the DFT) of $K$ will always match $M$ at the individual points $m_i$. Between these fixed points, however, it can attain very unintuitive and undesirable values. This undesirable behavior can easily be visualized by appending a large number of zeros to $K$ and then calculating the absolute value of the DFT. Changing any $p_i$ changes the behavior of the magnitude response of $K$ in between the $m_i$. Thus, two FIRs derived from the same magnitude spectrum $M$, but with different phase spectra $P$, can have very dissimilar magnitude responses. This is illustrated by the figure below:
The upper plot shows three randomly chosen examples of $P$, called $P_0$, $P_1$ and $P_2$. The lower plot shows a randomly chosen example of $M$ and the resulting magnitude responses when an FIR kernel is formed using $M$ as the magnitude spectrum and one of the $P_i$ as the corresponding phase spectrum. $n$ was chosen to be $16$ here. The three magnitude responses all go through the fixed set of points defined by $M$, but vary considerably everywhere else, due to their different phase spectra.
There are many audio applications where we do not really care about the phase response of a filter, but we absolutely do care about the magnitude response. Given a magnitude spectrum $M$, we can hence regard all the $p_i$ as degrees of freedom to shape our magnitude response. More specifically, this should allow us to define $n$ additional magnitude points that the magnitude response of $K$ is required to go through. These additional requirements should be able to help to substantially reduce any undesirable magnitude response behavior of $K$ without the need to increase its length $n$. Calculating $P$ to satisfy such additional requirements seems nontrivial, though.
Are my assumptions correct? Am I missing something? Is there any practical (polynomial complexity) way to do this?
Edit:
This seems to come down to the question of wether a complex polynomial can feasibly be interpolated by a number of absolute value requirements. The coefficients of the FIR kernel $K$ define a polynomial in the complex plane and $K$'s frequency response lies on the unit circle. Any magnitude requirements are positions on the unit circle with a value which the absolute value of that polynomial is supposed to go trough. Consequently, I asked a more specific question on Mathematics Stack Exchange. If I gain any further insight from there, I will update or answer my own question here.