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