Deceptive linear algebra problem Does anyone recognize this problem? There are $2p$ equations and $2p$ unknowns, and it feels like a classic, but I've never encountered it before:
Given $d_1, d_2, \ldots d_p$, find $a_1, a_2, \ldots a_p$ and $n_1, n_2, \ldots n_p$:
$$ \left[ \begin{array}{ccc}
d_1\\
d_2\\
d_3\\
\vdots\\
d_{2p}
\end{array} \right]
= \left[ \begin{array}{ccc}
a_1, ~a_2, \cdots ~a_p\\
a_1^2, ~a_2^2, \cdots ~a_p^2\\
a_1^3, ~a_2^3, \cdots ~a_p^3\\
\vdots\\
a_1^{2p}, ~a_2^{2p}, \cdots ~a_p^{2p}
\end{array} \right] \cdot
\left[ \begin{array}{c} 
n_1\\
n_2\\
\vdots\\
n_p
\end{array} \right] $$
(and where $a \succeq 0$ and $n \succeq 0$ and only real values are considered, so wherever there can be positive and negative roots for the values of $a_i$, only the positive roots need be considered). 
It seems straightforward to solve with elimination: From the first row, let $a_1 = \frac{d_1-\sum_{i>1} a_i n_i}{n_1}$. Then from the second row, solve for $a_2$ (which will now involve a quadratic, etc.). Because the polynomials get larger, solving it with elimination in this manner can be quite slow (I crashed sage doing this on a fairly small problem).
Has this family of problems been characterized well in the literature and can it be solved efficiently in practice? Even a numerical solution would suffice, as long as it will converge certainly (no multi-start Levenberg-Marquardt, etc.)...
 A: Actually, the problem precisely as you've stated it can't be solved: you can't write $(1,0,0,0)$ as a sum of multiples of two vectors $(a,a^2,a^3,a^4)$ and $(b,b^2,b^3,b^4)$, but this is easily fixed by starting with the zeroth power instead of the first. (Even more degenerately, $(1,0)$ isn't a multiple of $(a,a^2)$ for any $a$.)
Dongryul Kim's answer is basically correct, but let me expand a little bit on it: the span of the vectors $(1,a,\dots, a^{n-1})$ for $a\in A$ for some finite subset $A$ in any field has dimension equal to the size of $A$ if $\# A< n$, and spans the whole space otherwise (note how important it is that we start at the zeroth power to make this work).  If $\# A< n$, then the orthogonal to this space is spanned by vectors of the form $(0,\dots,0,e_{\# A},\dots, e_1,1,0\dots,0)$ where $\prod_{a\in A}(X-a)=X^{\# A}+e_1X^{\# A-1}+\cdots+e_{\# A}$ (these are obviously orthogonal to the vectors, and have the correct dimensional span).  Thus, if we want to write $(d_1,\dots, d_n)$ as a sum of these vectors, we have to find numbers $1,e_1,\dots,e_q$ such that $(d_1,\dots, d_n)$ is orthogonal to all of these vectors.  That is, we want to consider the equation 
\begin{equation} \begin{bmatrix} d_{q+1} & d_{q} & \cdots & d_1 \\ d_{q+2} & d_{q+1} & \cdots & d_2 \\ \vdots & \vdots & \ddots & \vdots \\ d_{n} & d_{n-1} & \cdots & d_{n-q} \end{bmatrix} \begin{bmatrix} 1 \\ e_1 \\ \vdots \\ e_q \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix}. \end{equation}
For $q=1$, this will only have a non-trivial solution if all the $d_i$'s are 0.  Otherwise, there will be a unique smallest integer $q>1$ such that this equation has a non-trivial solution.  In that case, we have that $(d_1,\dots, d_n)$ is in the span of the vectors $(1,a,\dots, a^{n-1})$, for $a$ ranging over the roots of $X^{q}+e_1X^{q-1}+\cdots+e_q$, and one can easily solve for the $n_i$'s; by construction, we've seen that these roots are the unique minimal set such that $(d_1,\dots, d_n)$ is in the span.  By basic linear algebra, we must have  $q\leq \lceil n/2\rceil$, since at that point, the matrix has more columns than rows, and must have a non-trivial solution. Obviously, generically, we'll have equality there.
A: Again not an answer, but maybe a useful attack line. Your system is a more general version of the ones to be solved to determine nodes and weights of the Gaussian quadrature. Indeed, suppose that $d_i = \int_{a}^b x^i \omega(x) \mathrm{d}x$ for each $i$, and for a certain weight function $\omega(x)$; then, you are looking for a set of nodes $a_i$ and weights $n_i$ such that the quadrature formula $\sum_{i=1}^p n_i f(a_i) \approx \int_{a}^b \omega(x) f(x) \mathrm{d}x$ is eaxct for all polynomials of degree at most $2p-1$.
Look for the derivation of the Gaussian weights and nodes; the theorems there might be generalized to your setting. If you only have one of these systems to solve, maybe you are lucky and it corresponds to a known classical weighting function (Hermite, Laguerre...).
A: Let $$f(X) = (X - a_1) \cdots (X - a_p) = X^p + e_1 X^{p-1} + \dots + e_p.$$ Using $f(a_i) = 0$ for each $i$, we get 
\begin{equation} \begin{bmatrix} d_{p+1} & d_p & \cdots & d_1 \\ d_{p+2} & d_{p+1} & \cdots & d_2 \\ \vdots & \vdots & \ddots & \vdots \\ d_{2p} & d_{2p-1} & \cdots & d_p \end{bmatrix} \begin{bmatrix} 1 \\ e_1 \\ \vdots \\ e_p \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \end{equation}
Note that the dimension of the first matrix is $p \times (p+1)$. In the general case, it will be possible to determine $e_1, \dots, e_p$ uniquely. Then calculate the roots of $f(X)$ to obtain $a_i$, and the remaining $n_i$ part will be a simple elimination problem. 
If the dimension of the null space of the '$d_i$ matrix' is greater than $1$, it probably means that $a_i$ are not pairwise distinct. Then you might want to check whether system is also solvable with a smaller number of variables. 
A: Not an answer, but a suggestion. Eliminate the $n_i$'s. So let $A$ be the top $p$-by-$p$ part of your matrix and let $B$ be the bottom $p$-by-$p$ part, also let $b$ be the top $p$ entries of the $d$-vector and let $c$ be the bottom $p$ entries. Then you have $b=An$ and $c=Bn$, so $c=BA^{-1}b$. The next thing to observe is that
$$
  B = A\begin{pmatrix} a_1^p & 0 & \ldots & 0\\
        0 & a_2^p & \ldots & 0 \\
       \vdots &&\cdots &\vdots \\
      0 & 0 &\ldots &a_p^p \end{pmatrix}.
$$
So you get $p$ equations for the $a_1,\ldots,a_p$ that look like
$$
  c = A\begin{pmatrix} a_1^p & 0 & \ldots & 0\\
        0 & a_2^p & \ldots & 0 \\
       \vdots &&\cdots &\vdots \\
      0 & 0 &\ldots &a_p^p \end{pmatrix} A^{-1} b.
$$
The matrix $A$ is a vanderMonde matrix, and there are formulas (of sorts) for its inverse, for example in Macon, N.; A. Spitzbart. "Inverses of Vandermonde Matrices". The American Mathematical Monthly 65 (2): 95–100. https://www.jstor.org/stable/2308881 .
