MathGroup Archive 2006

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Vandermonde Matrix/Optimization question

  • To: mathgroup at smc.vnet.net
  • Subject: [mg71967] Re: Vandermonde Matrix/Optimization question
  • From: Paul Abbott <paul at physics.uwa.edu.au>
  • Date: Wed, 6 Dec 2006 06:03:59 -0500 (EST)
  • Organization: The University of Western Australia
  • References: <ejevma$1vc$1@smc.vnet.net>

In article <ejevma$1vc$1 at smc.vnet.net>, scottmacd123 at gmail.com wrote:

> I am trying to find the numerical values of x[i] variables that
> maximize the determinant of a Vandermonde matrix. 

What is the application? 

> For small N x N matrices where N <= 7, Mathematica has no problem doing 
> this calculuation. Indeed:
> 
> VandermondMatrix[n_]:=Table[x[i]^j,{i,1,n},{j,1,n}]

This is not the standard definition of the Vandermonde matrix. See

  http://mathworld.wolfram.com/VandermondeMatrix.html

Since you are multiplying the i-th row of the Vandermonde matrix by 
x[i], your determinant is

  Product[x[i], {i, 1, n}]

times that of the standard Vandermonde matrix.

> vars[n_] := Table[x[i], {i, 1, n}]
> cons[n_] := Table[-1 <= x[i] <= 1, {i, 1, n}]
> NMaximize[{Det[VandermondMatrix[4], cons[4]}, vars[4], Method ->
> {"NelderMead","RandomSeed"->100}]
> 
> gives:
> 
> {0.366453, {x[1] -> 1., x[2] -> -1., x[3] -> 0.654654, x[4] ->
> -0.654654}}
> 
> I want to do this same calculation for larger N, say N=100. Clearly the
> code I have above will fail because I am asking Mathematica to
> calculate the determinant of a 100x100 matrix symbolically, which is
> not feasible. Obviously the determinant must be evaluated numerically.

Not true. The determinant of the standard Vandermonde matrix is 
especially simple -- it is just the product of the differences between 
the x[i]'s. Hence, the determinant of your matrix is (up to a sign) just

  Product[x[i] Product[x[i] - x[j], {j, i + 1, n}], {i, 1, n}]

So, your optimization problem is to maximize the value of the product of 
n values, x[i], multiplied by their differences x[i] - x[j], with each 
x[i] in [-1, 1].

> While the above does not work, I hope it is clear what I want to do:
> have NMaximize fill in numerical values for the x[i] variables before
> calling Det[].

Even though computing the Det is simple, I expect that NMaximize still 
will struggle to determine the optimal x[i] for very large n. However, I 
have found a closed-form solution to your problem for _even_ n using 
"experimental mathematics". Analytic and numerical experiments showed 
that +1 and -1 always arise as optimal values. For even n, one observes 
that the optimal values arise as pairs of opposite sign,

  x[k] == - x[n - k + 1]

for k = 2, 3, ..., n/2. Using these results it is straightforward to 
compute accurate values of the x[i] for n = 2, 4, ..., 14. Using

  <<NumberTheory`Recognize`

one can compute the polynomials whose roots are the optimal x[i]. The 
optimal x[i] turn out to be the roots of an associated Legendre 
polynomial,

 x[n_?EvenQ][i_] := 
   Root[Function[x, Sqrt[1 - x^2] LegendreP[n, 1, x]/x], i]

or, equivalently (but slightly simpler and better, since more is known 
about orthogonal polynomials), the roots of a Gegenbauer polynomial:

 x[n_?EvenQ][i_] := 
   Root[Function[x, GegenbauerC[n + 1, -1/2, x]/x], i]

I expect that this is a known result, though I have not managed to find 
it in the literature. A search for "Gegenbauer" and "Vandermonde" turned 
up some relevant results, but I have not yet found an explicit statement 
of this result. Another connection is that 

  Integrate[LegendreP[n, x], x] / x

generates these polynomials for even n.

For example, here are the exact x[1], ..., x[4], for n = 4.

  Table[x[4][i], {i, 4}]

  {-1, -Sqrt[3/7], Sqrt[3/7], 1}

Clearly these agree with your values above.

Here is a high-precision value of x[2],

  N[x[100][2],30]

  -0.999273257811212882702294983064549830406931695104

and of x[10], both for n = 100.

  N[x[100][10],30]

  -0.958521719914778750791263985570009960595780493123

For large n, there are asymptotic formulas for the zeros of Gegenbauer 
polynomials in terms of the zeros of Bessel functions. See Chapter 22 of 
Abramowitz and Stegun at

  http://www.convertit.com/Go/ConvertIt/Reference/AMS55.ASP

After loading

  <<NumericalMath`

we obtain the required zeros.
  
  zeros = BesselJZeros[1, 100];

Indexing these as
  
  j[i_/; 1 <= i <= 100]:= zeros[[i]]

we implement the approximate formula AS 22.16.2

  X[n_][i_] := 1 - (j[i-1]^2 (1 + 1/(n + 1)))/(2 (n + 1)^2)

Then we obtain (apart from a sign, but recall the sign-symmetry of the 
roots for even n),

  X[100][2]
  0.9992732410210722

and

  X[100][10]
  0.9582358644272112

There is a large literature on the topic of zeros of orthogonal 
polynomials and I expect that you can find better asymptotic 
expressions. However, Mathematica can compute the roots of very large 
polynomials exactly, and to arbitrary precision so, at least for even n, 
your problem is "solved".

Cheers,
Paul

_______________________________________________________________________
Paul Abbott                                      Phone:  61 8 6488 2734
School of Physics, M013                            Fax: +61 8 6488 1014
The University of Western Australia         (CRICOS Provider No 00126G)    
AUSTRALIA                               http://physics.uwa.edu.au/~paul


  • Prev by Date: Re: Re: Re: Reduction of Radicals
  • Next by Date: RE: Re: RE: RE: Functional decomposition (solving f[f[x]] = g[x] for given g)
  • Previous by thread: Re: Re: are there any methods of figuring out how "large" a piece of typeset textual data will be?
  • Next by thread: need help defining a MakeBoxes rule for Piecewise