How to find the roots of large number of polynomials arranged in matrix array
14 views (last 30 days)
Show older comments
%suppose that
[x1,x2,x3,x4,x5]=meshgrid(linspace(0,1,101)) ; % i have computer with 500 GB RAM
%and we have the following polynomial
P=a*x^4+b*x^3+c*x^2+d*x+e ;
% where a,b,c,d,and e depend on x1,x2,x3,x4,and x5.
% i need to find the roots for all the polynomial when x1,x2,x3,x4,and x5 varied in the region 0 to 1
% the problem is that the "roots()" function does not work on the array matrices and if i use the "for loop" the simulation will take more than 7 days
% any help and suggestion will be approciated
0 Comments
Accepted Answer
Stephan
on 7 Nov 2018
Edited: Stephan
on 7 Nov 2018
Hi,
% Matrix containing 3 polynomials of order 4
A = [1 1 -2 2 -3; 2 3 -2 -2 1; 0 0 -4 2 -7]
% Convert to cell array of single polynomials
B =mat2cell(A,ones(1,size(A,1)),size(A,2))
% calculate roots for every polynomial and save in new cell array
C = cellfun(@roots,B,'UniformOutput',0)
% Access the roots of the first polynomial
C{1,:}
Best regards
Stephan
More Answers (2)
Bruno Luong
on 6 Nov 2018
5 Comments
Bruno Luong
on 7 Nov 2018
Edited: Bruno Luong
on 7 Nov 2018
I don't know the algorithm behind ROOTS, but I know it can computed as eigenvalues of the companion matrices.
All TMW needs is to have transform EIG function to vectorized form. Something that I want to do in my free time.
Walter Roberson
on 7 Nov 2018
I did some reading to see if I could figure out an efficient algorithm for your purpose.
What I found is that it has been shown that there are polynomials for which a small change in coefficients makes a big change in roots, like 10^29. Some of them are pretty simple polynomials to construct too. Small changes on the order of the last bit in single precision can drive to complex roots as well.
The implication is that as you marginally change the coefficients you can get very different roots or different numbers of real roots. So it is not nearly as simple as just tweaking the roots slightly with small coefficient changes.
Now if it were the case that you could fix your leading coefficients but varying the constant offset then there would be efficiencies available in calculation over multiple constant values. However your coefficients are functions of your parameter x so small changes in x could lead to changes to all of coefficients. I therefore think that the technique of decomposition would not be productive.
About all I have come up with is that you could start with roots() for one location, and then for adjacent locations fsolve the polynomial using the original roots as starting points. fsolve instead of fzero because fsolve is rated for finding complex roots but fzero is real only.
It appears there are formulas available that permit calculating bounds on the absolute values of roots so the search could be constrained in theory.
Steven Lord
on 7 Nov 2018
Assuming that we replace meshgrid in your code with ndgrid, each of x1 through x5 will be 5-dimensional arrays with 101^5 elements. Each one of those will require nearly 80 GB of contiguous memory to store.
Now you want to use these five 5-D arrays to generate polynomial coefficients. Assuming that you can't perform in-place manipulation (that blog post is old but still relevant) that will mean you need several more 80ish GB blocks of contiguous memory.
Now let's assume that constructing the coefficients and computing the roots of these polynomials take on average 1/1000 second. How long would it take to compute all the roots of your matrix serially?
>> days(seconds(101^5/1000))
ans =
121.6447
About 4 months. On the plus side, this calculation should be amenable to parallelization: give each of 101 workers responsibility for a specific value of x1 and have them compute the roots for each of the polynomials generated from that value of x1. So now you're looking at 29 hours on that cluster.
If you tell us what your ultimate goal is in trying to compute the roots of a large collection of polynomials, we may be able to suggest an alternative that doesn't require you to wait until the first week in March or create a 101 worker cluster.
4 Comments
See Also
Categories
Find more on Polynomials in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!