Looking for a very fast root finding method in Matlab

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.

Answer 1

[Bruno Luong]

% 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)

[John D'Errico]

Bruno was too fast for me. But this is the idea.

[Bruno Luong]

That's certainly because I'm usually careless to comment/explain as much as you John. ;-)

[Behnam]

Dear Bruno

I saved your code as a function named Bruno, not to forget your help

It works perfectly.

Thank you very very much.

[Paul]

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.

[Bruno Luong]

@Paul, yes good remark. The accuracy of this method is under caution. I wouldn't recommend this method for larger than 10 terms.

I can't run tzero since I don't own the toolbox.

I also suspect fx happens to be large when the zeros of f is close to the poles. The zeros is likely still accurate.

I would propose to run few Newton's step to improve the accuracy of the solutions, if they are not multiple roots.

a=rand(1,10);
b=rand(1,10);
c=rand(1,10);

% assume a ~= 0
nnewton = 3;
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);
% Newton steps to refine solutions
x = x(:).';
for n=1:nnewton
    d = 1./(bn(:).*x+cn(:));
    fx = sum(d,1);
    dfx = -sum(bn(:).*d.^2,1);
    x = x-fx./dfx;
end
toc 

% Check
x = x(:).'
fx = sum(a(:)./(b(:).*x+c(:)),1)

[Paul]

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]

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]

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

But again, this uses symbolic computations. Instead, we should build the coefficients directly.