Decomposing a diagonal positive real matrix

I would like to 'decompose' a diagonal positive real matrix $E$ of rank $D$ onto $\sum_{i=1}^{D}c(i)N^i$: $$E = \left( \begin{array}{ccc} 0 & & & \\ & a & & \\ & & b & \\ & & & \ldots \\ \end{array} \right) \Rightarrow \sum_{i=1}^{D}c(i)N^i$$ where $$N = \left( \begin{array}{ccc} 0 & & & \\ & 1 & & \\ & & 2 & \\ & & & \ldots \\ \end{array} \right)$$

and $\bigl\{c(0), c(1), c(2), \ldots \bigr\}$ is the coeffients we want to find given $E$. Just a note about $E$: it is well behaved and $0<a<b<\ldots$. We can also assume $E$ is such that $c(i) < c(i+1)$.

First I just re-express the problem as a set of linear equations, where $E[i]$ is the $i$th eigenvalue of the rank $D$ matrix and $D_{max}$ is the truncation of the expansion (I assume I have truncate because of numerical difficulty in solving the equations). $D$ should be of $\mathcal{O}(10^3)$:

$$E[p] = \sum_{i=1}^{D_\max}c(i)p^i : p\in(1,\ldots,D)$$

This sounds more complicated than the rather simple problem I think it should be. Now in terms of Mathematica my toy attempt is as follows:

 d = 10;
Dmax = 5;
(* Generate matrix E *)
evals = Table[1.25* x - 1.21*10^-5 x^2 + 10^-6 x^3 + 0*10^-8 x^7, {x, 1, d}]
(* Decompose *)
eqns = Table[ evals[[p]] == Sum[c[i] p^i, {i, 1, Dmax}], {p, 1, d}];
solns = NSolve[eqns, Table[c[i], {i, 0, Dmax}], Reals]  ;
Chop@solns


So as long as this is soluble, this works and gives:

 (* {{c[1] -> 1.25, c[2] -> -0.0000121, c[3] -> 1.*10^-6, c[4] -> 0, c[5] -> 0}} *)


That is as long as 0 multiples $10^{-8} x^7$. Otherwise this is not a problem to solve but to minimize (which is the case of interest). (The NSolve here is just to demonstrate the idea). Using NMinimize with the constrains: $c[1] > c[2] > c[3] > c[4] > c[5]$:

 eqns = Table[Abs[evals[[p]] - Sum[c[i] p^i, {i, 1, Dmax}] ],          {p, 1, d}];
solns = NMinimize[{Total[eqns],  Greater @@ Table[c[i], {i, 1, Dmax}]}, Table[c[i], {i, 1, Dmax}]]


 (* {995.29, {c[1] -> 0.91, c[2] -> 0.468, c[3] -> 0.247, c[4] -> 0.245, c[5] -> -0.030}} *)


This is quite far from the result above, and it is much worse if the $10^{-8}x^7$ term is included.

Does anyone have any ideas on how to do this more correctly or efficiently?

• You end up with a system of equations involving a Vandermonde matrix, which is fixed for a given dimension $D$. You can perhaps precompute the inverse of this and just multiply evals by it. I get similar results with this approach.
– rm -rf
Dec 21, 2013 at 5:09
• This sounds very interesting! Perhaps you can put up some code, because I am not sure what exactly you mean at the moment. Also, It would be nice for the code to work for an approximate decomposition, in some sense a projection of from the {N^1,N^2,..N^D} onto {N^1,N^2,..N^Dmax}. as I tried to simulate with 10^−8 * x^7 Dec 21, 2013 at 5:42
• Thinking about it, do you mean: LeastSquares[ Table[ p^i , {p, 1, d}, {i, 1, Dmax}], evals]. That seems to work, except when I make D = 100 (It must be the large number 100^dmax). Any suggestions on how to get around this Dec 21, 2013 at 5:59

The set of equations given by

$$e_p = \sum_{i=1}^D c_i\ p^i,\quad p\in\{1,\ldots,D\}$$

where the $e_p$ are the diagonal elements (or eigenvalues) of $\mathbf{E}$ can be written as

$$\mathbf{e} = \mathbf{c}\ \left( \begin{array}{ccc} 1 & 1 & \ldots & 1^{D-1} \\ 1 & 2 & \ldots & 2^{D-1} \\ \vdots & & \ddots & \vdots \\ 1 & D & \ldots & D^{D-1} \\ \end{array} \right) \left( \begin{array}{ccc} 1 & & & \\ & 2 & & \\ & & \ddots & \\ & & & D \\ \end{array} \right)$$

where $\mathbf{e}=\{e_1,\ldots,e_D\}$ and $\mathbf{c}=\{c_1,\ldots,c_D\}$. Observe that the first matrix is a Vandermonde matrix that is explicitly invertible and that the two matrices are dependent only on the rank $D$. So you can precompute the inverse matrix and obtain $\mathbf{c}$ for different $\mathbf{e}$.

d = 10;
e = Table[1.25*x - 1.21*10^-5 x^2 + 10^-6 x^3 + 0*10^-8 x^7, {x, 1, d}];
Vinv = DiagonalMatrix[1/Range@d].LinearAlgebraVandermondeInverse[Range@d];
Chop[e.Vinv]
(* {1.25, -0.0000121, 1.*10^-6, 0, 0, 0, 0, 0, 0, 0} *)


where Vinv is the inverse of the product of the two matrices above.

This also works when you include the $10^{-8}x^7$ term:

e = Table[1.25*x - 1.21*10^-5 x^2 + 10^-6 x^3 + 10^-8 x^7, {x, 1, d}];
Chop[e.Vinv]
(* {1.25, -0.0000121, 1.*10^-6, 0, 0, 0, 1.*10^-8, 0, 0, 0} *)

• Could you also say what the LinearAlgebraVandermondeInverse[Range@d]; function does, it is undocumented. because Range@d is not a matrix but a 1D list. Dec 27, 2013 at 5:25
• @AimForClarity Why don't you try evaluating it? :) Also try LinearAlgebraVandermondeMatrix@Range@d (for an appropriate value of d`). Note that a square Vandermonde matrix is fully specified by a 1D list, and a unique inverse exists if all the elements of the 1D list are distinct.
– rm -rf
Dec 27, 2013 at 6:50
• thanks, i'll just need to pick up a nice book on it Dec 28, 2013 at 17:36