Suppose I am given some non-negative constants $(C_p)_{p=1, ..., l}$ and I would like to find an integer $N$ and vector $v \in R^N$ such that $$ \sum_{i=1}^N (v_i)^p = C_p $$ for $p=1, ..., l$.
- Can I find $v$ constructively?
- How small can I take $N$?
It seems that for a generic right-hand side the $N = l$, but I cannot prove that.
I used to numerically solve a similar problem using the following algorithm I call continuous point insertion. The idea is to construct solution gradually, starting from the simpliest system $$ v_1 = c_1 $$ and adding one more equation at a time.
Assume that we've already solved the system with $k$ equations, so we know $\mathbf{v}^{(k)} \in \mathbb R^k$ satisfying system $$ \sum_{i=1}^k v_i = C_1\\ \sum_{i=1}^k v_i^2 = C_2\\ \vdots\\ \sum_{i=1}^k v_i^k = C_k\\ $$ Now, let's insert a new variable $v_{k+1}$ along with the $k+1$ equation parametrically $$\begin{align} \sum_{i=1}^k v_i + v_{k+1} &= C_1 + (1-\gamma) A_{k+1}\\ \sum_{i=1}^k v_i^2 + v_{k+1}^2 &= C_2 + (1-\gamma) A_{k+1}^2\\ \vdots\\ \sum_{i=1}^k v_i^k + v_{k+1}^k &= C_k + (1-\gamma) A_{k+1}^k\\ \sum_{i=1}^k v_i^{k+1} + v_{k+1}^{k+1} &= \gamma C_{k+1} + (1-\gamma) A_{k+1}^{k+1}+ (1 - \gamma)\sum_{i=1}^k (v_i^{(k)})^{k+1} \end{align} $$ The $A_{k+1}$ is some arbitrary parameter, an initial guess for $v_{k+1}$. When $\gamma = 0$ the system reduces to the previous system for which the solution is known $$ \mathbf{v}^{(k+1)}\big|_{\gamma = 0} = \left(v_1^{(k)}, v_2^{(k)}, \dots, v_k^{(k)}, A_{k+1}\right). $$ When $\gamma = 1$ is it the system we want to solve. Let $\mathbf{v}(\gamma)$ be the smooth (in $\gamma$) solution to the system above. Denoting the system as $\mathbf{F}(\mathbf{v}(\gamma)) = \mathbf{g}(\gamma)$ we can obtain a system of ODE for $\mathbf{v}(\gamma)$: $$ \mathbf v'(\gamma) = \left(\frac{\partial \mathbf F(\mathbf v(\gamma))}{\partial \mathbf v}\right)^{-1} \mathbf g'(\gamma) $$ with initial condition $$ \mathbf v(0) = \left(v_1^{(k)}, v_2^{(k)}, \dots, v_k^{(k)}, A_{k+1}\right) $$ This technique is known as numerical continuation.
Note that this ODE was extremely stiff and ill-conditioned for my problem, so I had to use double-double and quad-double precision to integrate it for sizes about $n = 20$.
The method of integration I used was a combination of explicit Euler method with Newton refinement at each step:
When we reach $\gamma = 1$ we solved the system with $k+1$ equation. Repeat the step until the system of $l$ equations is solved.
Note that both the Euler step and Newton refinement use the same Jacobi matrix $$ \frac{\partial F_p}{\partial v_i} = p v_i^{p-1} $$ which is not singular unless some $v_i$ are the same (it's a variation of the Wandermonde matrix). So we need to avoid crossing of the $v_i(\gamma)$ trajectories. Careful choice of $A_{k+1}$ may help here.