Vectorization time-varying recursive linear function

I try to vectorize this simple recursive relation (all quantities are scalars)

x_{0} = 0;
x_{n} = x_{n-1}*a_{n} + b_{n} for n=1,2,...,N

In MATLAB code it can be carried out by for loop

% test inputs
b=rand(1,10);
a=0.9+zeros(size(b));

xk=0;
x=zeros(size(b));
for k=1:length(x)
    xk = a(k)*xk+b(k);
    x(k) = xk;
end

For a(:) constant this can be vectorized by IIR filter

ac = unique(a);
if length(ac)==1
    x = filter(1, [1 -ac], b);
end

I would though it could have some time-varying IIR filter that I can use to vectorize the case where a is time-dependent.

But I couldn't find anywhere such stock function. anyone have an idea?

Answer 1

[David Goodmanson]

Hi Bruno,

a = rand(1,50);
b = rand(1,50);

% method 1
xk = 0;
x = zeros(1,50);
for k = 1:50
    xk = a(k)*xk + b(k);
    x(k) = xk;
end

% method 2
cpa = cumprod([1 a(2:end)])
x1 = filter(1,[1 -1],b./cpa).*cpa;

max(abs(x1-x))
  ans =  4.4409e-16

[Bruno Luong]

Thanks David, very clever workaround.

The problem is that I might have some zeros in A, in that case the filter returns NaN onwards.

[David Goodmanson]

Hi Bruno,

Also, if one of the a's is nonzero but very small, there are probably going to be numerical accuracy issues. It's unfortunate that Matlab apparently does not have a built-in function for this type of iteration.