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