Decomposing a diagonal positive real matrix

Question

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}]]  

Which gives the bad result:

 (* {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?

Answer 1

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

Taking your example above:

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].LinearAlgebra`VandermondeInverse[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} *)