# Speed up calculation of recursively defined list

I have two lists $$a$$ and $$b$$ of length $$n$$ and $$n-1$$ respectively (typically I have $$n \approx 1000$$). I have to compute a list $$\theta$$ of length $$n$$ which is defined recursively by the following relation: $$\theta_i = a_i\theta_{i-1} - b_{i-1}^2\theta_{i-2}, \;\;\;\;\; i = 2,...,n$$ where $$\theta_0 = 1$$ and $$\theta_1 = a_1$$. I can easily code a simple loop that performs this operation:

theta[[1]] = a[[1]];
theta[[2]] = a[[2]]*theta[[1]] - b[[1]]^2;
Do[
theta[[i]] =
a[[i]]*theta[[i - 1]] - (b[[i - 1]]^2)*theta[[i - 2]];
, {i, 3, n}];


where I have computed explicitly $$\theta_2$$ for simplicity.

So far so good, but I know that the Do[] loops are usually not the best way to go in Mathematica, so I wonder if there is some way to speed up this code. Also consider that I absolutely need to insert this code into Compile[], where the lists $$a$$ and $$b$$ are the input of my function. This is because I want to give the attribute Listable to my function for performance reasons.

I know that sometimes there are clever ways to code recursive functions in general, but I really don't know much about this topic, so any help related to some basic concepts of recursive procedures is also appreciated. Thanks!

• Please post a[i] and b[i]. Nov 23, 2022 at 23:52
• a and b can be whatever, this is meant to be a small step inside a function where a and b are the input variables. Nov 24, 2022 at 0:28

A common approach to speed recursive calculations like this is memoization.
In Mathematica, this can be as simple as using the f[arg_]:=f[arg]=computation pattern:

listOfAs = RandomReal[{0, 1}, 1000];
listOfBs = RandomReal[{0, 1}, 1000];

Clear[theta]
theta[i_] := theta[i] = With[{a = listOfAs[[i]], b = listOfBs[[i - 1]]},
a theta[i - 1] - b^2 theta[i - 2] ]
theta[0] = 1; theta[1] = listOfAs[[1]];

AbsoluteTiming[theta/@Range[1000];]


This only takes a fraction of a second on my machine, so I'd be curious what application you had in mind and whether it indeed requires Compile.

• Thanks! Essentially I am comparing several approaches to compute the top-left element of the inverse matrix of a tridiagonal symmetric matrix. I need the fastest possible algoritm because I apply the function to thousands of matrices. So far I simply use the cofactor method, which is not listable,, but then I found that this can be done with such recursive methods and I was wondering if those could be compiled and made listable Nov 24, 2022 at 0:26
a[i_] = i;
b[i_] = -i;
RecurrenceTable[{θ[i] ==
a[i] θ[i - 1] - b[i - 1]^2 θ[i - 2], θ[0] == 1,
θ[1] == a[1]}, θ, {i, 1, 1000}]