Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
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 '13 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 –  AimForClarity Dec 21 '13 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 –  AimForClarity Dec 21 '13 at 5:59

1 Answer 1

up vote 4 down vote accepted

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} *)
share|improve this answer
    
Could you also say what the LinearAlgebra`VandermondeInverse[Range@d]; function does, it is undocumented. because Range@d is not a matrix but a 1D list. –  AimForClarity Dec 27 '13 at 5:25
    
@AimForClarity Why don't you try evaluating it? :) Also try LinearAlgebra`VandermondeMatrix@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 '13 at 6:50
    
thanks, i'll just need to pick up a nice book on it –  AimForClarity Dec 28 '13 at 17:36

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.