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>
- Re: Re: Re: Problems with NSolve
- References:
- Re: Problems with NSolve
- From: "Kshitij Wagh" <kgwagh@gmail.com>
- Re: Problems with NSolve