How to find the roots of large number of polynomials arranged in matrix array

14 views (last 30 days)
%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

Accepted Answer

Stephan
Stephan on 7 Nov 2018
Edited: Stephan on 7 Nov 2018
Hi,
you could use cellfun to apply roots to your Matrix:
% 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
Bruno Luong on 6 Nov 2018
A vectorized code can be derived from this formula for a~=0.
For a = 0, use Cardano formula
  5 Comments
Bruno Luong
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
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.

Sign in to comment.


Steven Lord
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
Oussama GASSAB
Oussama GASSAB on 7 Nov 2018
thank you so much for your answers, I really appreciate them. I think the problem is solved by using all your important comments. I would like to calculate the roots because I want to find some roots that satisfy some conditions (i need only a few numbers of roots). I have reduced the interval linspace(0,1,10) to find the required regions for each xn and then I repeat the simulation in that regions until I get the required results. it takes only 2s simulation. thank you, everyone, I really felt that I m not alone

Sign in to comment.

Categories

Find more on Polynomials in Help Center and File Exchange

Products


Release

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!