How to efficiently control an FIR's magnitude response by altering its phase spectrum 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 &lt; i &lt; \frac{n}{2}$) and $P$ is conjugate symmetric (i.e. $p_i = - p_{n-i}$ for $0 &lt; i &lt; \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.
There are a few things that can make this easier


If you don't care about the phase response you can typically get the best result with calculating a minimum phase for your amplitude spectrum and using that as the complex target
If you want to control the in-between and the behavior at the band edges you need to specify targets points there as well. 
If you have more target points then FIR coefficients, you have an over determined system but you can still solve it with a least square error approach. this will not hit all your target points exactly but it's often a good compromise between filter fidelity and computational complexity. 
A least square error approach has the added advantage that you can apply weights to steer the solution towards a trade off that fits your requirements. For example, you typically need to control the out of  band behavior: you don't care about the exact behavior other than it shouldn't have too much gain. To achieve this, you can set the target simply to zero and play around with the weight until the gain in that band feels acceptable