MathGroup Archive 2008

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

Search the Archive

Re: Re: Problems with NSolve

  • To: mathgroup at
  • Subject: [mg88240] Re: [mg88214] Re: Problems with NSolve
  • From: Daniel Lichtblau <danl at>
  • Date: Tue, 29 Apr 2008 06:49:56 -0400 (EDT)
  • References: <fuumfl$8dr$> <> <>

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> wrote:
>>kgwagh at 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.
>>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 :-)
>>-- 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 

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 

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

  • Prev by Date: Re: Does Mathematica really need more printed, introductory documentation?
  • Next by Date: Complex Plot
  • Previous by thread: Re: Problems with NSolve
  • Next by thread: Re: Re: Re: Problems with NSolve