Re: Re: Problems with NSolve

*To*: mathgroup at smc.vnet.net*Subject*: [mg88240] Re: [mg88214] Re: Problems with NSolve*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Tue, 29 Apr 2008 06:49:56 -0400 (EDT)*References*: <fuumfl$8dr$1@smc.vnet.net> <48131E36.1090106@gmail.com> <200804280841.EAA06082@smc.vnet.net>

Kshitij Wagh wrote: > Hi > Sorry for the "fuzziness" :) > chpoly was found by using Det[A - x*IdentityMatrix[dim]] (where A is a > matrix function of u with \alpha etc. being some random parameters). \Alpha > etc are random real numbers between -1 and 1 with working precision 30. Also > please find attached a notebook of a working sample. > > [contact the sender to get this attachment - moderator] > > So I think the problem seems here, that instead of getting polynomials, I > am getting rational functions. The thread Carl sent mentions some > alternatives - I will try them out. > Thanks, > Kshitij > > On Sat, Apr 26, 2008 at 8:21 AM, Jean-Marc Gulliet < > jeanmarc.gulliet at gmail.com> wrote: > > >>kgwagh at gmail.com wrote: >> >> >>>Hi all >>> >>>I am trying to find the number of eigenvalue crossings for a matrix as >>>a function of the parameter 'u', on which the elements of the >>>(symmetric) matrix depend on linearly. The matrix elements also >>>involve randomly chosen constants. The plan is to find the >>>distribution of the crossings of these type of matrices as I scan over >>>the random numbers. >>> >>>So far I have been using the following : >>> >>>NSolve[{chpoly[u, dim, \[Alpha], \[Gamma], \[Epsilon], x] == 0, >>> D[chpoly[u, dim, \[Alpha], \[Gamma], \[Epsilon], x], x] == 0}, {x, >>> u}, WorkingPrecision -> prec] >>> >>>where chpoly is the characteristic polynomial of the matrix (with the >>>eigenvalue variable being x) and \alpha, \gamma and \epsilon are >>>constant parameter arrays of random numbers. For a double (or higher) >>>degeneracy of the eigenvalues both the characteristic equation and its >>>derivative should be zero. This approach has worked successfully only >>>upto 12*12 matrices (where one such computation takes 40 secs on my >>>laptop). For 13*13 my laptop takes 4000 sec. This seems to be somewhat >>>surprising, because these polynomials are of the order 'n' (where n >>>is dimension of the matrix) in both u and x - and n=13 does not sound >>>very computationally unreasonable. So I was wondering if there was a >>>faster approach I could take. >>> >>>Also, the problem essentially entails me to know the number (not the >>>values) of the real solutions to this system of polynomials. >>>CountRoots seemed ideal but it does not work for more than one >>>equation. So is there any alternative along this route? >>> >>>Any other alternatives are also welcome. >>>Thanks, >>>Kshitij. >>> >>> >> >>Could you post, please, a *fully working* example of what you are doing? It >>is really hard to optimize fuzzy functions (fuzzy not as in fuzzy logic but >>as in fuzzy definitions or concepts :-) >> >>Regards, >>-- Jean-Marc >> >> Without a specific matrix it is still fuzzy. I will show an example that should get you going, but it might not behave quite the same with actual parameters in place (due to numerical issues). First some code to generate a matrix of random polynomials, linear in some variable 'u'. I'm working over the integers but I'd guess this will all be fine if I use random reals of precision 30. randomLinearPolynomial[n_, u_] := RandomInteger[{-n,n}.{1,u}] randomMatrix[dim_,max_] := Table[ randomLinearPolynomial[max,u], {dim,dim}] Now to address what I believe is the main bottleneck. You are computing bivariate determinants, and Det is not using any method adapted to the specifications of your particular family of matrices. For dimension higher than 11 or so, it will do some thing that might well create denominators (at least, if inputs are not exact), and moreover I expect would be slow. Given that matrix elements are guaranteed to be linear polynomials in a particular variable, you can get a characteristic polynomial much faster via interpolation. One does this by finding the characteristic polynomial for n+1 specific values of that variable (where n is the dimension), using some other variable which I call 't' in the code below. To recover the correct coefficients in t, each of which is a polynomial of degree up to n in u, we'll use standard polynomial interpolation. myCharacteristicPolynomial[mat_, u_, x_] := Module[{n = Length[mat], polys, coeffs}, polys = Map[CharacteristicPolynomial[mat /. u -> #, x] &, Range[1, n + 1]]; coeffs = Transpose[Map[ PadRight[CoefficientList[#, x], n + 1] &, polys]]; coeffs = Map[Expand[ InterpolatingPolynomial[#, u]] &, coeffs]; Expand[coeffs.x^Range[0, n]] ] Here is an example: In[59]:= myMatrix = randomMatrix[5, 10] Out[59]= {{8 - 8*u, -2 + 10*u, 1 + 2*u, 10 + 5*u, -5*u}, {-3 + 7*u, -9 + 4*u, -1 - 8*u, -4 + 5*u, -6 - 4*u}, {10 - 5*u, -4*u, 6 - 5*u, -9, 9 - 7*u}, {7 + 5*u, 7 + 7*u, 5 - 6*u, -5 + 7*u, -6 - 5*u}, {-2 + u, 2 - 7*u, 1 - 7*u, 7 + 2*u, -6 + 5*u}} In[60]:= Timing[c1 = CharacteristicPolynomial[myMatrix, t]] Out[60]= {0.00399999999999778, -82800 + 14894*t + 2266*t^2 + 71*t^3 - 6*t^4 - t^5 + 260131*u - 14113*t*u - 3792*t^2*u - 129*t^3*u + 3*t^4*u - 184171*u^2 - 33715*t*u^2 + 987*t^2*u^2 + 299*t^3*u^2 + 195079*u^3 + 28987*t*u^3 + 133*t^2*u^3 - 113971*u^4 - 17800*t*u^4 + 15563*u^5} In[61]:= Timing[c2 = myCharacteristicPolynomial[myMatrix, u, t]] Out[61]= {0.00400000000000237, -82800 + 14894*t + 2266*t^2 + 71*t^3 - 6*t^4 - t^5 + 260131*u - 14113*t*u - 3792*t^2*u - 129*t^3*u + 3*t^4*u - 184171*u^2 - 33715*t*u^2 + 987*t^2*u^2 + 299*t^3*u^2 + 195079*u^3 + 28987*t*u^3 + 133*t^2*u^3 - 113971*u^4 - 17800*t*u^4 + 15563*u^5} In[62]:= c1 == c2 Out[62]= True No difference in timings, for dimension of 5. Now we redo at higher dimension. In[63]:= myBiggerMatrix = randomMatrix[11, 10]; In[65]:= Timing[c1 = CharacteristicPolynomial[myBiggerMatrix, t];] Out[65]= {4.38027, Null} In[66]:= Timing[ c2 = myCharacteristicPolynomial[myBiggerMatrix, u, t];] Out[66]= {0.064004, Null} In[67]:= c1 == c2 Out[67]= True So we have made substantial improvement in regard to this bottleneck. Next is the issue of finding, or at least counting, the real solutions to the system that results when we look for multiple eigenvalues. One approach, as you indicate, is to solve it and count the real solutions. findcrossings1[poly_, u_, t_] := Cases[{u, t} /. NSolve[{poly, D[poly, t]}, {u, t}], {_Real, _Real}] In[70]:= Timing[solns = findcrossings1[c1, u, t]] Out[70]= {21.0053, {{-0.301473, 778.797}, {-8.11089, -140.741}, {3.86324, 75.9679}, {7.27872, 80.1776}, {3.77224, 65.1796}, {4.55024, 55.6788}, {-0.740926, 18.6957}, {0.853927, 20.6467}, {-6.79623, 8.78344}, {-0.871391, 13.772}, {0.266494, -14.1587}, {0.135089, -13.0663}, {-0.260045, -13.2098}, {-0.860066, 11.1762}, {0.655937, 9.11692}, {0.967805, -4.66542}, {1.3373, 7.43368}, {0.00336349, -5.3958}, {-0.0986887, -4.9421}, {0.0149776, -3.17222}, {0.18333, -1.30277}, {-0.0141373, 0.155747}}} In[71]:= Length[solns] Out[71]= 22 This is a bit slow, and while I expect it to not deteriorate as fast as CharacteristicPolynomial, it will all the same become problematic as dimension rises. Note that we have two polynomials in two unknowns. A faster, and basically reliable, alternative, is to find the resultant of the these polynomials, and use real root counting on that. In[72]:= findcrossings2[poly_, u_, t_] := CountRoots[Resultant[poly, D[poly, t], t], u] In[73]:= Timing[findcrossings2[c1, u, t]] Out[73]= {0.94806, 22} There is a very small caveat to this. In a nongeneric selection of your parameters \[Alpha] et al, you could have a case where a real solution in the linear matrix variable u corresponds to crossings of two (nonreal) complex conjugate pairs of eigenvalues. In this nongeneric case findcrossings2 would not recognize that the real solution in u corresponded to complex crossings in t, hence would overcount. For your purposes of assessing statistics, this probability zero situation matters not at all. Daniel Lichtblau Wolfram Research

**Follow-Ups**:**Re: Re: Re: Problems with NSolve***From:*Daniel Lichtblau <danl@wolfram.com>

**References**:**Re: Problems with NSolve***From:*"Kshitij Wagh" <kgwagh@gmail.com>

**Re: Does Mathematica really need more printed, introductory documentation?**

**Complex Plot**

**Re: Problems with NSolve**

**Re: Re: Re: Problems with NSolve**