Looking for a very fast root finding method in Matlab

Question

Behnam 2020-10-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/603937-looking-for-a-very-fast-root-finding-method-in-matlab

编辑： Bruno Luong 2020-10-6

采纳的回答： Bruno Luong

在 MATLAB Online 中打开

I am trying to solve an equation like below in Matlab

ai, bi, and ci are vectors of length 10 and these vectors change through the code. this code has to be executed more than 1 million times. So, the method for finding x must be very fast. I tried the following codes:

method 1:

s=tf('s');%transfer function
sys=0;
for i=1:10
    sys =sys+a(i)/(s*b(i)+c(i));
end
[roots1] = zpkdata(sys,'v');

In this case, the for loop takes 0.1 second to be completed on my PC and [roots] = zpkdata(sys,'v');executes in less than 0.005 seconds. So, preparing the equation sys before being solved by zpkdata takes a long time for a million times run in the real case. Seemingly vectorized operation does not work for 'tf' argument type.

Next, I checked the following method: method 2:

syms x
sys=sum(a./(x*b+c));
[roots1]=vpasolve(sys,x)

This symbolic method was again slow and took 0.13 seconds to get executed. Do you know any fast method suitable for my case? Or, do you know a way to prepare sys more quickly in method 1?

Many tnx.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bruno Luong 2020-10-2

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/603937-looking-for-a-very-fast-root-finding-method-in-matlab#answer_504133

编辑：Bruno Luong 2020-10-2

在 MATLAB Online 中打开

% randomm coefficients for test
a=rand(1,10);
b=rand(1,10);
c=rand(1,10);
% assume a ~= 0
tic
n = length(a);
bn = b./a;
cn = c./a;
P = 1;
Q = [bn(1) cn(1)];
for k=2:n
    P = [P*bn(k), 0] + [0, P*cn(k)] + Q;
    Q = [Q*bn(k), 0] + [0, Q*cn(k)];
end
x = roots(P);
toc % Elapsed time is 0.000385 seconds.
% Check
x = x(:).'
fx = sum(a(:)./(b(:).*x+c(:)),1)

7 个评论
显示 5更早的评论隐藏 5更早的评论

Paul 2020-10-4

编辑：Paul 2020-10-4

在 MATLAB Online 中打开

Bruno,

How concerned should the OP be about the accuracy of your approach? Especially if n gets large and/or the coefficients of the final polynomial have a large dynamic range? I'm asking because I was inspired by the OP's control theoretic approach and came up with a method that on average isn't quite as fast as yours, but does seem to offer more accuracy, at least for the problem domain we are using as an example. Here is the code:

>> type useroots
function [x,fx] = useroots(a,b,c)
% assume a ~= 0
bn = b./a;
cn = c./a;
P = 1;
Q = [bn(1) cn(1)];
for k=2:numel(a)
    P = [P*bn(k), 0] + [0, P*cn(k)] + Q;
    Q = [Q*bn(k), 0] + [0, Q*cn(k)];
end
x = roots(P);
% Check
fx = max(sum(a(:)./(b(:).*x(:).'+c(:)),1));
end
>> type usetzero
function [x,fx] = usetzero(a,b,c)
x = tzero(diag(-c./b),a(:)./b(:),ones(1,numel(a)),0);
% Check
fx = max(sum(a(:)./(b(:).*x(:).'+c(:)),1));
end

And here is the test code:

% random coefficients for test
n = 10;
ntrials = 1000;
a = rand(ntrials,n);
b = rand(ntrials,n);
c = rand(ntrials,n);
% a = [1 2];
% b = [3 4];
% c = [5 6];
[dtroots,dttzero,fxroots,fxtzero] = deal(zeros(ntrials,1));
for ii = 1:ntrials
    
    tic
    [xr,fxroots(ii)] = useroots(a(ii,:),b(ii,:),c(ii,:));
    dtroots(ii) = toc;
    
    tic
    [xz,fxtzero(ii)] = usetzero(a(ii,:),b(ii,:),c(ii,:));
    dttzero(ii) = toc;
    
end
figure(1);plot(1:ntrials,[dtroots dttzero],'-x'),grid
figure(2);plot(1:ntrials,abs([fxroots fxtzero]),'-x'),grid
mean([dtroots dttzero],1)
max([fxroots fxtzero],[],1)
max(rmoutliers([fxroots fxtzero]))
sum(fxtzero < fxroots)/ntrials

Here are the results:

ans =
   9.1427e-05   1.0471e-04
ans =
   5.5042e+02   4.2086e-06
ans =
   1.4760e-07   1.0820e-11
ans =
   9.9200e-01
 

Also, depending on the problem domain of the OP, either approach would have to deal with special cases of a_i = 0 (as you noted) or b_i = 0.

Paul 2020-10-5

在 MATLAB Online 中打开

I looked more closely at the worst case that I saw and I agree that the zeros seem to be pretty accurate. Is any zero close to a pole? I'm not sure. Here's the data that I see:

>> index = 334;
disp('the data')
[a(334,:);b(index,:);c(index,:)]
[xr,fxr]=useroots(a(index,:),b(index,:),c(index,:));
[xt,fxt]=usetzero(a(index,:),b(index,:),c(index,:));
xr = sort(xr); xt = sort(xt);
disp('roots');
[xr xt]
disp('fx for each root');
[sum(a(index,:)./(b(index,:).*xr + c(index,:)),2)  sum(a(index,:)./(b(index,:).*xt + c(index,:)),2)]
disp('evaluate the denominators at the fifth zero')
[((b(index,:).*xr(5) + c(index,:))) ; ((b(index,:).*xt(5) + c(index,:)))]
the data
ans =
  Columns 1 through 8
     5.403041674081475e-01     1.127955018909095e-01     2.091123542456026e-01     5.291713447535034e-01     8.596427432788267e-02     4.162109095747117e-01     6.974328024049009e-01     8.109584297027117e-02
     3.914807928171162e-01     1.146364750887373e-01     6.027741708170893e-01     5.186778923606967e-01     8.162768010754857e-01     4.668434760127809e-01     7.865056595354460e-01     6.845644168380302e-01
     3.077652099388081e-01     3.663647990555764e-01     3.056739234359613e-01     4.173917963553463e-01     6.568463572900645e-01     4.881625070399680e-01     7.140215217186643e-01     5.928806834423205e-01
  Columns 9 through 10
     1.577279766062935e-01     4.783351991950024e-01
     7.740945923505923e-01     7.475420835940132e-01
     4.959658147894286e-01     4.048304825339897e-01
roots
ans =
    -2.839538690133413e+00    -2.839538690133357e+00
    -1.007931214723476e+00    -1.007931215204822e+00
    -8.819362558135367e-01    -8.819362340251731e-01
    -8.599283416923192e-01    -8.599283771198578e-01
    -8.046893095146248e-01    -8.046892524985675e-01
    -7.960872390157562e-01    -7.960872818401908e-01
    -6.531383428012292e-01    -6.531383430055718e-01
    -5.815367632129854e-01    -5.815367630496052e-01
    -5.157640198572923e-01    -5.157640198874921e-01
fx for each root
ans =
     5.113687251423471e-13     3.552713678800501e-15
    -3.763499218933930e-07    -3.634870182622763e-13
     4.741643357686343e-05    -1.363353874239692e-13
    -1.484119453740718e-04    -2.256861364458018e-12
     5.504241178784788e+02    -4.750323030489056e-06
    -1.251901757180818e-03     2.383870878475136e-12
    -3.136275452675363e-07    -7.380762667708041e-13
     9.649716403714592e-08    -3.055333763768431e-13
    -1.706899546149998e-07     2.131628207280301e-12
evaluate the denominators at the fifth zero
ans =
  Columns 1 through 8
    -7.255198921435058e-03     2.741180530712299e-01    -1.793720078720927e-01     1.724129111635442e-05    -2.858140174866186e-06     1.124985526758360e-01     8.112882561774160e-02     4.201901553864418e-02
    -7.255176600743762e-03     2.741180596073497e-01    -1.793719735042861e-01     1.727086408481737e-05    -2.811599290053479e-06     1.124985792934105e-01     8.112887046119344e-02     4.201905456980826e-02
  Columns 9 through 10
    -1.269398282281745e-01    -1.967086405464007e-01
    -1.269397840923528e-01    -1.967085979244985e-01
    

The fifth zero is close to the fourth and fifth poles. At first glance it appears that the two solutions really aren't that different wrt to those pole locations, but maybe that (seemingly) small difference is significant.

Bruno Luong 2020-10-6

编辑：Bruno Luong 2020-10-6

f is a sum of the inverse distances to the poles (with phase). If f(x) is large then x must be "close" to one of the pole.

请先登录，再进行评论。

Answer 2

John D'Errico 2020-10-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/603937-looking-for-a-very-fast-root-finding-method-in-matlab#answer_504127

编辑：John D'Errico 2020-10-2

在 MATLAB Online 中打开

Symbolic methods (vpasolve) are NEVER efficient.

But all you have is a rational polynomial. The roots of a rational polynomial are the roots of the numerator polynomial. So compute that, or actually the coefficients of that polynomial. Then call roots. You won't be able to do it more quickly than that.

n = 10;
a = rand(1,n);
b = rand(1,n);
c = rand(1,n);
 
drop = @(V,ind) V(setdiff(1:length(V),ind));
 
tic
P = sym(0);
syms x
for i = 1:n;
   P = P + a(i) * prod(x*drop(b,i) + drop(c,i));
end
toc
 
roots(sym2poly(P))
ans =
      -2.8076
      -2.0032
      -1.3514
      -1.2919
      -1.0956
     -0.86728
     -0.40917
     -0.29335
    -0.029945