       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

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, cons}, vars, Method ->
>
> gives:
>
> {0.366453, {x -> 1., x -> -1., x -> 0.654654, x ->
> -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, ..., x, for n = 4.

Table[x[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,

N[x,30]

-0.999273257811212882702294983064549830406931695104

and of x, both for n = 100.

N[x,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

<<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
0.9992732410210722

and

X
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,

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