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

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{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}$$ 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.

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{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}.$$ 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.

• Thanks a lot for the description, I wouldn't have understood Dongryul Kim's answer without it.
– user
Aug 6, 2015 at 2:57
• Hi Ben, do you know if this same approach would work if the columns of the $A$ matrix were of the form $(1,a,a^2,a^4,a^8...)$ instead (but still with unique powers $0,1,2,4,8...$)? When I use $a_1 = 0.7, a_2 = 0.2$ and $n_1=2,n_2=1$, I can recover the $a_i$ values from $d_0,d_1,d_2,d_3$ ($e=(1,-0.15,0.005)$), but cannot seem to recover them from $d_0,d_1,d_2,d_4$ ($e=(1,0.665,-0.056125)$), even when I try inserting 0s into $e$ to correspond to the terms $X^4, X^2, X^1, X^0$. Does the approach require an evenly spaced sequence of powers or does it generalize to unique powers? Vielen Dank!
– user
Aug 21, 2015 at 11:59
• @user Seems unlikely; the important property evenly spaces powers have is that their orthogonal is easily described. Aug 24, 2015 at 19:07

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

• This is actually very similar to the motivation I was using to approximate a more complex problem. Thanks for pointing me in the direction regarding quadrature.
– user
Aug 6, 2015 at 3:00

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 .

• Thank you, I tried to find the name of this type of matrix.
– user
Aug 6, 2015 at 3:00