Deceptive linear algebra problem

Question

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.)...

Answer 1

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.

Answer 2

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.

Answer 3

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...).

Answer 4

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 .