How to justify the solution of this problem? Assume $\mathbf{x} \in \mathbb R_+^N$ with support $P=\{p_1,p_2,\cdots,p_K\}$ ($P$ is unknown).
We already know that $$f_1(\mathbf{x}) = f_2(\mathbf{x}) = \cdots = f_{N-1}(\mathbf{x})$$
where 
$$f_l(\mathbf{x}) =\frac{ |DFT[x]|}{\sum_{i=0}^{N-1} x_i} =\frac{ |\sum_{i=0}^{N-1} x_i w^{il}|}{\sum_{i=0}^{N-1} x_i}, \qquad w:=e^{-j2\pi/N}$$
This problem is comprised of $2K$ unknowns ($P$ and non-zero elements of $\mathbf{x}$) and $N-2$ equations, therefore this is an over-determined system of equations, and generally it has Least Square solution.
After evaluation of LS solution for many $(N,K)$ pairs, I reached these results (facts):


*

*For almost all cases, the solution for $P$ was $P$(s) whose Difference Multiset has minimum variance. Where Difference Multiset of a point-set is defined : $D = \{a_1,a_2,\cdots,a_{N-1}\}$ and $a_d$ is the number of occurrences of $d = p_i-p_j \mod N, \quad i \ne j$

*The  $\mathbf{x}$ solution for all cases, was a solution with close non-zero elements, i.e. ($x_i \approx x_j, i,j\in P$)
This figure shows the LS solution for the $(N,K) = (20,7)$

Someone asked how I solved this problem. For finding the best $P$ ($P$ leading to least error in LS solution) I used combinatorial search and after finding $P$, the remaining problem is a continuous least square that is solvable using many optimization methods, such as Gauss-Newton, Levenberg-Marquardt, ... .
Here is how $f_l(x^*)$ is distributed for $l=0,1,\cdots,N-1$:

My question is, how to prove the first fact analytically or even intuitively?
I asked it in MO too.
 A: Here might be some intuition:  
If you write out your conditions, note that for
$\|x\|_1=1$,
$f_l(x)^2 = \sum_{i,\ j=0}^{N-1} x_i x_j \omega^{l(i-j)}
 = \sum_{k=0}^{N-1} \omega^{lk} \left(\sum_{i=0}^{N-1} x_i x_{(i+k\ \mathrm{mod}\ N)} \right)$
Let $c_k = \sum_{i=0}^{N-1} x_i x_{(i+k\ \mathrm{mod}\ N)}$. Then $f_l$ is constant for $l=1,\ldots,N-1$ if and only if the DFT of $(c_0,\ldots,c_{N-1})$
is constant (except for $l=0$). 
Noting that $\sum_{k=1}^{N-1} w^{lk} = -1 + \sum_{k=0}^{N-1} w^{lk} = -1$, we can compute that the DFT of $(b,c,\ldots,c)$ is $(b+(N-1)c,b-c,\ldots,b-c)$.
As the DFT is invertible, this implies that $f_l$ is constant in $l$ if and only if $c_k$ is constant for $k>0$. 
I think how I've arranged so far is more directly related to multisets. So now we're back to the same question, but now we can study $c_k$ rather than the original DFT magnitudes.

Some more thoughts:
Relaxing the problem means looking at full vectors (non-sparse) $x$. We can also ask what kind of vectors will the $c_k$ be constant for. It seems like there are $N$ unknowns and $N-2$ equations, but let us throw in $\sum_{i} x_i = 1$ (which is required for the $c_k$ formula that I wrote), and also $\sum_{i} x_i^2 = C^2$, which is just the missing term $c_0$, and vary $C$ to see what solutions look like. It's not even clear to me what this picture looks like, or if it is relevant. For $C^2 = 1/N$, then the constant solution satisfies the equation, and should be the only one. Otherwise, who knows? Once this is understood, then adding on the sparse constraint...  My guess as to what happens is that interesting behavior will occur around $C^2 \sim 1/k$ for each sparsity level $k$?
